Self-diffusion in metals is believed to be understood rather well and their coefficients obey a linear temperathe Arrhenius equation D (T ) = D0exp(-Q/k BT ) with ture dependency on an Arrehenius plot (log(D) vs. 1/T ).
temperature-independent prefactor D0 and activation energy Q, which for mono-vacancy diffusion is the sum of vacancy formation enthalpy E f and migration enthalpy Em. The anharmonicity was initially suggested to manifest itself in the form of soft phonons due to the specific distribution of the d-electrons [14,15]. By now, the continuous advances in atomistic simulation methods and computational capabilities have provided new opportunities to gain insight into the microscopic origin of anomalous diffusion, and lead to a variety of explanations for non-Arrhenius diffusion. Proposed mechanisms include a vacancy-interstitial model in Ti [16], anharmonic effects in both /3-Ti [17] and in Mo [18], and concerted atomic motion [19,20].
The previous work has been performed either with classical molecular dynamics (MD) based on empirical potentials or with density-functional theory (DFT). On the classical MD side, in Ref. [16] diffusion coefficients in /3-Ti were determined from the mean-square displacement of the atoms. The results suggested that interstitials would play a non-negligible role at self-diffusion near melting temperature and the sum of Arrhenius-like vacancy and interstitial diffusions with different slopes would result in the observed anomalous diffusion. The spontaneous formation of Frenkel pairs was also found in Ref. [20], leading to temperature-dependent formation energies. While no diffusion coefficients were calculated in Ref. [20], it was argued that the temperature-dependent formation energy would give rise to anomalous diffusion. The fact that the high-temperature defect concentrations varied by several orders of magnitude between Refs. [16] However, a number of elements show a remarkable "curvature" and thus deviation from linearity. While such an anomaly is rare for non-bcc metals, it has been found to be especially strong in the group IVB metals Ti, Zr and Hf, whose bcc-phases are mechanically unstable at lower temperatures, while it is much weaker in other bcc metals like Mo and Nb [1] that are mechanically stable for all temperatures.
Focusing on bcc Mo (in the following "Mo") and bcc Ti (in the following "/3-Ti") as prototypes for mechanically stable and unstable bcc-metals, the underlying mechanism(s) for such a non-Arrhenius behavior has been under debate for many decades and remains unsettled. Given that self-diffusion in elementary metals is governed by vacancy jumps, initial explanations were mostly based on the contribution of "secondary" diffusion mechanisms and include diffusion via divacancies [2], next-nearest-neighbor (NNN) jumps [3,4], diffusion via interstitials [5] or diffusion enhancement by phase transformations [6,7]. While none of these mechanisms could be conclusively confirmed for either Mo or /3-Ti, it has been argued in the late 1980s that secondary mechanisms would not be consistent with experiments [8,9], and it was suggested that self-diffusion should be dominated by the traditional mono-vacancy jumps through nearestneighbor (NN) sites in both Mo and /3-Ti [9][10][11]. This and [20] as well as the fact that interstitial concentrations have been found exceedingly small in Mo even at its melting temperature casts doubt on this explanation. In another MD study [19], it was proposed that concerted motion of atoms significantly contributed to diffusion at temperatures > 83% of the melting temperature, in addition to the otherwise dominant vacancy diffusion which followed a straight Arrhenius behavior over the entire examined temperature range from 73% to 97% of the melting temperature. While the considerable variation in results of these studies may at least partially originate from the different interatomic potentials used [21][22][23], the fact that experiments show pronounced non-linearity already in the sub-1600 K range where all MD-studies find Arrhenius like behavior make further studies desirable.
Other than classical MD, density functional theory (DFT)-based methods do not rely on empirical interatomic potentials. Performing harmonic transition-state theory based DFT calculations of diffusion coefficients started to become standard in the late 1990 for semiconductors [24] and were adopted into the metallic systems a decade later [25][26][27][28]. Mattson et al. [18] proposed a large quadratic temperature dependence of the vacancy formation energy in bcc Mo based on ab intio molecular dynamics (AIMD) simulations, but did not study the temperature dependence of the formation entropy. Further, the two data points obtained at > 2400 K were insufficient to address the nonlinearity right above the transition temperature. Non-MD based first principles calculations on the high temperature bcc phases of group IVB metals like titanium, zirconium and hafnium have remained quite challenging due to their mechanical instability at zero-temperature [29,30], with a manifestation of decreasing energies and negative curvatures as a function of atomic displacement [30] or strain [29] in certain directions at 0 K. Since phonon frequencies are determined by square roots of energy versus displacement curvatures, this corresponds to vibrational modes with imaginary frequencies. Since within the common harmonic or quasiharmonic approximations, vibrational free energies in solids are calculated from phonon frequencies, the appearance of imaginary frequencies impedes free-energy calculations at finite temperatures for mechanically unstable crystals.
A number of approaches have been suggested to overcome the limitations from mechanical instabilities through molecular-dynamics [31] or phonon calculations on cells with self-consistently displaced atoms [32] to capture the temperature-induced anharmonicity, but most of them are computationally too expensive to perform calculations on supercells with point defects. Kadkhodaei et al. employed a combination of self-consistent ab initio lattice dynamics (SCAILD) and the temperature dependent effective potential (TDEP) method and examined in a thorough study the influence of phonon anharmonic-ity on the diffusion coefficient [17] in /3-Ti, but did not include thermal expansion and left the lattice parameters fixed. While the results could explain differences in magnitude between modeling and experiment that arise when straight harmonic transition state theory is used, it did not reproduce non-Arrhenius behavior, suggesting thermal expansion as a prime suspect to cause anomalous temperature dependence.
Herein, we perform for a comprehensive study of the thermal expansion and self-diffusion anomaly in the mechanically stable Mo and the unstable /3-Ti with an accurate and efficient ab initio approach. Our work utilizes the large displacement method (LDM), first proposed by Antolin et al. [30], from which we demonstrate that the thermal expansion coefficients are well captured in the temperature range where the quasi harmonic approximation (QHA) applies. For the minimum-volume at each temperature, we then calculate the free energy of vacancy formation and diffusion coefficients within LDM-based harmonic transition state theory, effectively forming a quasiharmonic version of it. Our predicted self-diffusion coefficients and their anomalous temperature dependence in both Mo and /3-Ti agree well with previous measurements, a strong indication that going beyond harmonic to quasiharmonic transition state theory explains anomalous temperature dependence of self-diffusion in a natural way that eliminates the need to invoke ad-hoc non traditional assumptions about the diffusion mechanism. Finally, we quantified the non-linearity in the self-diffusion coefficients and compared it with both experimental data and previous computational work. We show that the nonlinear Arrhenius self-diffusion arises mostly from the thermal expansion, and further pinpoint that it is not the migration enthalpy, but the formation entropy and, in the case of /3-Ti, attempt frequency that carries the majority of the temperature dependence. This work uncovers the nature of anomalous self-diffusion in Mo and /3-Ti, and validates that LDM is an effective way to perform such calculations in mechanically unstable metals.
In the mono-vacancy mechanism of self-diffusion in the bcc lattice, the atom-vacancy exchange jump happens between two nearest neighbor positions along ½ [111]. The bcc self-diffusion coefficient can be written as
where g is the geometrical factor (1 for cubic lattices), f the correlation factor (0.7272 for bcc [33,34]), a the lattice constant, Cv the vacancy concentration and r the atom-vacancy exchange jump frequency along ½ [111].
The atomic fraction of vacancy lattice sites Cv is given r-point calculations as described in previous work [37]. by There, the entropy is approximated by
where !'!G f is the free energy of vacancy formation within the bcc lattice. !'!S f and !'!H f are the vacancy formation entropy and enthalpy, respectively. The atomvacancy exchange jump frequency r can be calculated within transition state theory as [35]
where g [v(V )] is the phonon density of states [37] and the r-point phonon frequencies v(V ) are calculated for a supercell with volume V .
The attempt frequency v * can be calculated as The climbing image nudged elastic band method (CI-NEB) [38] with 3 images was employed to determine the transition state (saddle point) structures for Mo and Ta (see Sec. II 4 why Ta), with the saddle point configuration of Ti determined from those as described in the following section due to the relaxation problems in mechanically unstable high-temperature structures. The initial and final structures, i.e. stable states, were fully relaxed first. Then during the CI-NEB calculations, all supercell volumes were fixed. ----------------------------------------------------------------
Therefore, 4 parameters will be calculated to obtain self-diffusion coefficients through DFT calculations: !'!H f , !'!S f , !'!Hm, and v * . Of those, !'!H f and !'!Hm contribute to the thermal activation energy Q, while !'!S f and v * contribute to the prefactor D0. We will describe their calculation in Sec. II 2. Before that, the vacancy structures of Mo and /3-Ti need to be relaxed, which is a straightforward task for Mo, but requires some extra attention for Ti, which we will discuss in Sec. II 4. The basic parameters for the DFT calculations are then discussed in Appendix II 8, the large-displacement method in Sec. II 3, and the calculation of lattice expansion to introduce vibrational anharmonicity in Sec. II 5.
The large displacement method (LDM) treats the anharmonicity through large atomic displacements into the high-temperature anharmonic range of the atoms' energy wells. Within LDM, phonons are calculated from the cur-vature of the harmonic envelope of the quartic energy-vs.-displacement curve (red curve in fig. 1). At high tem-peratures, or above the transition temperature, the atoms actually sample at the large vibrational amplitudes, which makes them carry across the "hump" at zero displace-ment (black circles in fig. 1), i.e. the local energetic bar-rier, and thus overcomes the mechanical instability of the bcc phase, therefore eliminating all imaginary fre-quencies from the phonon dispersions. Technically, this approach thus works like any standard finite-difference phonon calculation with the exception that the displace-ments are one or two orders of magnitude larger. The only question there is how large the displacement should be. In [30], a large displacement value of 0.88 Å was pro-
Formation quantities such as !'!G f , !'!S f or !'!H f can be calculated by the respective difference between perfect supercell and the relaxed cell with one vacancy. As an example,
where X is the lattice atom, V indicates a vacancy and N indicates the number of lattice sites in the supercell.
In order to calculate the vibrational contribution to the free energies of formation at finite temperatures, Ff (T ), and the entropy of formation, S f (T ), we use the quasiharmonic approximation (QHA) [36] based on supercell this problem, we first perform an approximate, nearestneighbor (NN) only, relaxation based on free energies for Ti, where free energy calculations based on the LDM (Sec. II 3) at three different temperatures between phase transformation temperature (1150 K) and melting temperature (1940 K) were performed for differently displaced NN shells in an otherwise unrelaxed vacancy cell. The minimum-energy displacement was then linearly extrap- derivative of energy on phonon displacement (
Ti at zero temperature. A harmonic potential well (red curve) is obtained by matching energy and curvature of the energy data at a large displacement value of 0.88 Å as determined in Ref. [30]. This represents the high temperature regime (black dashed line), where the atomic vibration exhibits large amplitude (blue marker) with high thermal energy, resulting in the envelope harmonic potential governing its vibration and real phonon fre-
In the quasiharmonic version of transition state theory we employ in this paper, the four variables that determine vacancy self-diffusion according to Eq. 4 are !'!H f , !'!S f , !'!Hm and v * . They are evaluated at the minimum-lattice constant in the Helmholtz free energy for each temperature, determined in turn from the traditional quasi-harmonic approximation to vibrational entropy. The Helmholtz free energy of a system as a function quencies ( d 2 E > 0). du 2 and experimental phonon dispersion. The sensibility of this large-displacement value for /3-Ti proposed is further confirmed by the comparison of the theoretical thermal expansion coefficients with the experimental values (section III 1). As was also shown in Ref. 30, regular moderately small displacements (we use 0.05 Å) are the most sensible choice for Mo, which is mechanically stable at zero temperature. of volume (V ) and temperature (T ) is expressed as f
where E is the total (internal) energy at 0 K, in addition to the quantities defined for Eq. 6. As described in Sec. II 3, we used the LDM for the Ti phonon calculations and traditional displacements for Mo.
In order to determine the thermal expansion coefficient, we fit the lattice constant data vs. temperature with a modified expression derived from Due to /3-Ti's mechanical instability, one cannot simply relax its supercell with one vacancy at zero temperature like is commonly done for other crystals to determine the ground state structure. If one tried, one would end up with unphysically large atomic displacements and a negative formation energy as shown in result Sec. III 2. Because of the same reason, the CI-NEB method also cannot be used to find the saddle point structure. To circumvent
where a0 is the lattice parameter at 0 K, 0D is the Debye temperature, and < is the Thacher function to approximate the Bose-Einstein weighted integral from the Debye specific heat, Eq. ( 6) from [40]. Other than in the original model, we give Ja a linear temperature de-
(2) Supercells with one vacancy (initial and final position, separated by one bond length) were relaxed to their ground state in the traditional way for Mo and Ta, while the structure for Ti was determined as described in Sec. II 4. Then, the CI-NEB (Sec. II 2) was used to determine the saddle point configuration for Mo and Ta, while Ti was again gained from the interpolation procedure.
(3) The lattice constants calculated in (1) were assigned to the ground state and transition state (saddle point) structures for Mo and Ti determined in (2).
(4) r-point phonon calculations using a finite difference method with displacements as described in Sec. II 3 were performed for all temperature dependent structures, along with zero-temperature total energies.
(5) Self-diffusion coefficients at different temperatures were calculated based on the results from (4) using Eqs. 4-7.
(6) Analysis of the source of the non-linearity of the diffusion coefficient was carried out by quantifying the non-linearity of its different components as defined in Eq. 4.
pendence, Ja (T ) = J0(1 + xT ) with J0 and x as fitting parameters in addition to 0D, since Ja, originally given by Ja = Kyk B/V0 (K is the compressibility; y the Grüneisen constant; and V0 the minimum volume of the cell) should be temperature dependent because both Grüneisen parameter and elastic constants [41] are temperature dependent. Without this, the typically observed increase in thermal expansion with temperature in the high-temperature limit is not described well enough. The linear thermal ex-
pansion coefficient (a) is finally calculated by a = from the temperature derivative of Eq. 10.
In order to quantify the non-linearity in our calculated and previous (experimental) data, we fit the logarithm of the diffusion coefficient in an Arrhenius plot with an equation analogous to bowing parameter dependent quadratic equations used for lattice parameters or elastic constants [42],
First-principles calculations were employed to perform all structural relaxations and CI-NEB runs as well as energy and phonon-frequency calculations. To describe ion-electron interactions, we used projector augmented plane-wave potentials [44] as implemented in the Vienna ab initio simulation package (VASP) version 5.4.1 [44,45]. Exchange and correlation contributions to the total energy were described by the Perdew-Burke-Ernzerhof (PBE) functional [46]. After convergence tests, we chose
(11) There, the first two terms are the linear interpolation between the logarithms of highest (l) and lowest (l) diffusion coefficients (respectively, Dh and Dl can be fitting coefficients), /3 = 1/k BT , and 1 is the "bowing paramea 54-atom supercell (3 × 3 × 3 conventional bcc unit cells) ter", i.e. the measure of the quadratic term, whose form is chosen to be zero at the extremal temperatures. An absolute non-linearity parameter A can then be used to compare the degree of non-linearity between the different elements independent of if the temperature ranges exactly overlap or not. Following Eqs. 3 and 9 in Ref. 43, we calculate A by with 4 × 4 × 4 Monkhorst-Pack k-point mesh for Brillouin zone integration. Cutoff energies of 300 eV for Mo and 700 eV for Ti were found to be necessary for convergence in the phonon calculations. For Ta, the cutoff energy was 300 eV as well. The settings were kept constant for all runs.
The results from the quasiharmonic approximation for the Helmholtz free energy vs. volume for temperatures between 0 and 1800 K in steps of 200 K for Mo and /3-Ti are shown in Fig. 2(a). The energy minima are determined by a fit based on the assumption of a harmonic dependence of energy on the volume and are indicated by crosses in the plot at the minima of the parabolas. The minimum-energy lattice constants for both materials as a function of temperature are shown in Fig. 2(b).
In order to determine the diffusion coefficients for Ti and Mo within the quasiharmonic transition state theory employed in this paper, we performed the following steps:
(1) The thermal expansion, i.e. the temperature dependence of the lattice constants, was calculated within the quasiharmonic approximation as described in Sec. II 5 based on perfect bcc supercells for Ti and Mo.
The results for thermal expansion in Mo, based on traditional phonon calculations, are shown in Fig. 2(a). Fitting Eq. 10 to the DFT data results in parameter values of
Following the procedure described in Sec. II 4, we have performed free-energy calculations of vacancy formation energies in constant-volume cells of /3-Ti (lattice constant the zero-temperature value from Fig. 2(b), 3.237 Å) for three different temperatures (1155 K (experimental transition temperature), 1941 K (experimental melting point) and 1548 K (halfway in-between) as a function of simultaneous nearest-neighbor displacement along the bond direction to the vacancy. For each temperature, the minimum-energy displacement and energy minimum were determined from a parabolic fit. The results are shown in Fig. 3.
Fig. 3c shows the free energy of formation as a function of temperature for the three calculation temperatures. In principle, the calculations performed here result in the Helmholtz free energy of formation at constant volume. In order to estimate the energy error with respect to the Gibbs free energy from the non-zero pressure that comes from neglecting the formation volume of the vacancy [20], we estimate the formation volume from the pressure in the VASP output for the NN-displaced cell and the bulk modulus and pressure derivative reported previously [30] through the pressure expression of the Murnaghan equation [49]. We find typical pressures of < 10 kbar, indicating a volume relaxation in supercell by less than 1%. Multiplying the formation volume that with the pressures, we find that the pV term contributes less than 0.04 eV to the free energy, and thus is small. Fitting a0 = 3.1514, 0D = 327 K, Ja = 5.91 × 10 -6 K -1 , and x = 6.42 × 10 -6 K -1 . The lattice constant at 298 K is calculated with Eq. 10 to be 3.155 Å, while the experimental lattice constant at room temperature is 3.147 Å. The thermal expansion coefficient of Mo is calculated from the temperature derivative of Eq. 10 and is shown in Fig. 2c in comparison to experimental data from Ref. 47. While the (weak) increase with temperature in the calculations is slightly lower than what is found in experiment, the general agreement is rater good with an average deviation in the linear regime between 600 and 1800 K of 5%.
The results for /3-Ti, which are based on phonon calculations from the LDM, are shown in Fig. 2(a) and(b). Since, as described in Sec. II 3, the LDM uses the harmonic envelope of the quartic energy-vs.-displacement curve irrespective of if the kinetic energy of the vibrating atoms is high enough to overcome the energy "hump" at zero displacement, which in the real system only happens above the transition temperature, it can also be used to calculate the energy vs. volume dependence of metastable /3-Ti below the transition temperature and thus is shown here for the entire range. As seen in Fig. 2(a), the calculated values are well described by a quadratic fit for all temperatures, as was the case for Mo. Thus, we can use the entire range for fitting Eq. 9 to get more reliable coefficients and find parameters of a0 = 3.2383, 0D = 474 K, Ja = 7.08 × 10 -6 K -1 , and x = 1.47 × 10 -4 K -1 .
Experimental values compiled in Ref. 48 are for pure Ti only available above the transition temperature from a-Ti to /3-Ti, which we have calculated previously to be at 1200 K with the LDM [30], within 4% of the experimental value of 1155 K. The calculated lattice constant of /3-Ti phase at 1200 K is 3.27 Å, compared to an experimental value at the same temperature of 3.33 Å. The calculated linear thermal expansion coefficient shows excellent agreement with the experimental values between 1200 and 1300 K with a considerable temperature depen-F (T ) = E f -TS to the three values, we find E f = 1.90 eV and S f = 4.3k B. Combining these results with the lattice site density, we find an equilibrium vacancy concentration of C eq = 1.8 × 10 21 cm -3 exp(-1.90[eV]/(k BT )).
For the unrelaxed cell, we had previously found with the same method E f = 2.05 eV and S f = 8.15k B [30]. Overall, the NN relaxation decreases the formation energy by 0.14 eV.
Next, we examine the NN-relaxation as a function of temperature and extrapolate from it a zero-temperature displacement which can be used for a fully relaxed surrogate structure with all neighbor shells relaxed. Having three temperatures, extrapolation can either be done linearly (dashed line in Fig. 3 a) or with a second-degree polynomial. We perform here both and take their average for our zero-temperature structure. For a linear extrapolation, we find a zero-temperature displacement of -0.035 Å, for a quadratic one -0.058 Å, with an average of -0.047 Å. All of these are considerably smaller than the minimum-energy zero-temperature relaxation of -0.101 Å, which is another expression of the mechanical instability of the structure and its strong stabilization by entropy. Following the procedure described in Sec. II 4 of interpolating fully relaxed cells of mechanically stable Mo and Ta vacancies with scaled lattice constants, which have NN relaxations of -0.008 Å and -0.082 Å in cells dence of da/dT = 2×10 -10 K -2 (Fig. 2c). For higher temperatures, it seems that anharmonic effect other than large displacements become important as already discussed in Ref. [30], causing a faster increase in thermal expansion in experiment than captured by LDM alone. Therefore, we restrict our calculations here to this temperature range.
Overall, we have successfully determined the temperature dependence of the lattice constants in good agreement with experiments, which we will use in the following as input for vacancy and diffusion-coefficient calculations. with otherwise clamped other neighbors, we then create a "relaxed" /3-Ti cell which is more or less the straight average of the atomic positions in relaxed Mo and Ta cells, scaled to the /3-Ti lattice constant. Both ground state structure and saddle point configuration were determined this way.
The formation energy from this process agrees with the formation energy from the interpolated structure discussed in Sec. II 4 within 0.01 eV, which is less than the uncertainties in our numbers from the different approximations and numerical calculations such as the 0.04 eV from enthalpy discussed in the present section. This suggests that relaxation beyond nearest neighbors adds only small changes to the calculated formation energies.
shown with previous experimental [50][51][52] and simulation [18] data in Fig. 4. Our qhTST calculations with varying lattice constants that include thermal expansion are represented by the six individual points (1400 -1900 K, with 100 K step). We limit our highest calculated temperature to two third of its melting point, which is the commonly accepted applicability range for QHA [53].
In order to determine the degree of non-linearity in our data, we use Eqs. 11 and 12 and determine a non-linearity parameter of A = 0.051. While our data thus show good absolute agreement with the experimental results, they have a somewhat higher degree of non-linearity than the experimental values of Maier at al. [52], for which we calculate A = 0.021, and also the simulation data of Mattson et al. [18] with A = 0.022.
In order to explore the origin of the bowing anomaly, we now investigate the temperature dependence of formation and migration enthalpies, as well as of formation entropy and attempt frequency. The results are shown in Fig. 6(a), (b) and (c). The calculated vacancy forma-
Our calculated self-diffusion coefficients for Mo calculated within the quasiharmonic transition-state theory are the curvature of the calculated self-diffusion coefficients in Fig. 4. In order to quantify the effect of the different contributions, we perform two fits of Eq. 11, first one for calculated diffusion coefficients where the prefactor is fixed to the average value of all calculated prefactors (which results in a value of D avg = 3.25 × 10 -6 m 2 /s), 0 and secondly an analogous set of diffusion coefficients where the activation energy is fixed to the average calculated value of E avg = 4.268 eV. For the first fit, we find 0 a non-linearity coefficient of which is 13% of the overall non-linearity, while for fixed Ea, we find A = 0.041, which is 87%. Thus, the non-linearity in vacancy diffusion in Mo is clearly dominated by the non-linearity of the diffusion prefactor as a consequence of anharmonicity from thermal expansion. To investigate the effect further, a similar analysis was made with fixed averaged attempt frequency, resulting in a change of non-linearity by less than 1%, making the non-linearity in the formation entropy the by far most important factor. This conclusion is different from the DFT-MD calculation by Mattsson et al. [18], in which the anomaly was mostly attributed to anharmonicity in the vacancy formation enthalpy !'!Hv . However, in their simulations, thermal expansion was not considered and all MD simulations were performed at the zero-temperature volume of vacancy and bulk, which neglected important volume-dependent anharmonic effects from thermal expansion. The fundamental proportionality between Grüneisen parameter, the prime indicator for anharmonicity, and thermal expansion [39] as well as our results here are strong indicators that thermal ex-tion enthalpy !'!Hv increases monotonously with temper-ature from 3.08 eV to 3.17 eV between 1400 and 1900 K (Fig. 6(a)) and is slightly larger than values proposed from analyzing experiments of 3.0 eV [44]. The calculated va-cancy migration enthalpy !'!Hm decreases slightly from 1.167 eV to 1.163 eV (Fig. 6(a)), slightly smaller than the reported migration enthalpy of 1.30 eV from experimental analysis [54]. The calculated attempt frequency v * also from 2.3 × 10 12 Hz to 2.2 × 10 12 Hz as temperature increases from 1400 to 1900 K (Fig. 6(c)). The temperature dependence of the calculated vacancy formation entropy !'!S f on the other hand shows a much more pronounced temperature dependence and is found to be the major rea-son for non-Arrhenius like diffusion in Mo. It decreases from 3.27k B at 1400 K to 2.79k B at 1700 K, before re-covering back to 3. pansion cannot be neglected when trying to understand anharmonic effects. Our findings are further supported by the resistivity measurements by Schwirtlich and Schultz [55], which find that the temperature dependence of the formation enthalpy is small and constant to within 0.1 eV between 2000 K and 2600 K, whereas Mattson et al. predict a much larger increase of 0.5 eV for that interval. Our calculation results of self-diffusion coefficients in /3-Ti are plotted in Fig. 5 in comparison to experimental data [15,56]. Calculations based on TST and large displacement algorithm are applied at five temperatures: 1200, 1230, 1250, 1280 and 1300 K. We choose again the upper limit of 1300 K as two thirds the melting point of Ti (1941 K), while the lower limit of 1200 K the ato-/3 transition temperature from LDM-DFT calculations [20]. Although these choices make the temperature range for the DFT calculations rather small, the anomaly of the self-diffusion coefficient can still be seen in Fig. 5, and the data are sufficient to extract the non-linearity parameter, which is independent of the temperature interval.
From the LDM-DFT calculations, we find for /3-Ti a non-linearity coefficient of 0.091, about 80% larger than what we had found for Mo. A similarly stronger non-linearity by about a factor of two is also found in experiments, where we fit a non-linearity parameter of A = 0.046, which is 120% larger than the Mo value. While thus these trends from experiment are reproduced well by the calculations, and while the values of our diffusion coefficients agree well with experimental values, the calculated non-linearity is considerably stronger than the experimentally observed one, which is true to a similar degree for /3-Ti and for Mo.
In order to explore if the origins of the bowing anomalies in /3-Ti and Mo differ, we again investigate the temperature dependence of formation and migration enthalpies, as well as of formation entropy and attempt frequency. The results are shown in Fig. 6(d),(e) and (f). The vacancy formation enthalpy !'!H f increases from 1.91 eV to 1.92 eV at the temperature from 1200 K to 1300 K (Fig. 6(d)), compared with an experimental estimation of > 1.50 eV [57]. This value from interpolated Mo and Ta zero-temperature cells (Secs. II 4and III 2) agrees well with the results from the nearest-neighbor only freeenergy minimization, where the nearly perfect linearity with temperature indicated a very weak temperature dependence of the formation enthalpy (Fig. 6). Similarly, the vacancy migration enthalpy !'!Hm increases only slightly from 0.146 eV to 0.147 eV as temperature increases (Fig. 6(d). Indeed, when we fit the calculated diffusion again only 15% of the overall bowing. Thus, it is once again the prefactor that we identify as the major source (85%) of the non-linearity. Unlike Mo, however, both vacancy formation entropy and attempt frequency contribute to a more comparable degree to the non-linear Arrhenius curve in /3-Ti. Specifically, the vacancy formation entropy !'!S decreases from 5.08k at 1200 K to a f B minimum of 4.92k B at 1230 K before increasing back to 4.95k B at 1300 K (Fig. 6(d)) and the attempt frequency v * decreases from 3.1×10 12 Hz at 1200 K to 2.9×10 12 Hz at 1250 K, then increases back to 3.1 × 10 12 Hz at 1300 K (Fig. 6(f)). With that, the formation entropy carries 64% of the non-linearity, leaving the remaining 36% to the attempt frequency. Without the anharmonic attempt frequency, the calculated non-linearity in /3-Ti would only be 25% larger than that in Mo, and thus the effect of the formation entropy is rather comparable in both materials.
In this paper, we have demonstrated that the framework of quasiharmonic transition state theory allows reproducing anomalous self-diffusion in bcc metals for the examples of Mo and /3-Ti by including thermal expansion into the calculations. In order to be able to calculate phonon frequencies needed for attempt frequency and formation entropy in the mechanically unstable /3-Ti, we use the previously proposed large-displacement method [30] and introduce two approaches to relax vacancies in mechanically unstable crystals, which is free-energy relaxation coefficients once again with an averaged prefactor (which has a value of D avg = 3.403 × 10 -5 m 2 /s for /3-Ti), we 0 only find a non-linearity parameter of A = 0.013, once Self-diffusion coefficient D (m 2 /s) is somewhat cumbersome. We find that phonon frequencies, at the heart of the two quantities that define the diffusion prefactor -vibrational entropy and attempt frequency -change more sensitively with volume changes than formation enthalpy and migration energy and thus dominate the non-linearity in the Arrhenius plots, which is in contrast to previous suggestions. In order to quantify the degree of non-linearity, we have proposed to use the non-linearity parameter from Ref. 43 and find that for both metals studied 85% of the non-linearity comes from the diffusion prefactor. While these values are similar, we find a considerable difference in the distribution of the non-linearity on formation entropy: While nearly all of the bowing in Mo comes from the formation entropy, the non-linearity in /3-Ti is divided 2/3:1/3 between formation entropy and attempt frequency, respectively, again in contrast to previous suggestions. Our calculated diffusion coefficients agree well with measured values, although we calculate a somewhat larger non-linearity than what we extract from experimental datasets. The proposed methodology is general enough that it also can be applied to other mechanically unstable crystals.
E (eV/atom)