Theoretical Predictions for the Equation of State of Metal Nickel at Extreme Conditions

1. Introduction

In the design of new materials, theoretical predicting the equation of state (EOS) without experimental data supporting is crucial to understand the thermodynamic properties. In fact, even for existing materials, the predictions are also necessary since experimental measurement of EOS still confronts some problems, especially at high pressures and high temperatures.

Taking nickel (Ni) as an example, it is one of the most abundant elements in the earth's core, where the temperature and pressure reaches 6000 K and 350 GPa, and Ni-based alloys have been widely applied in industries. However, rare experiments concern the EOS of Ni-based alloys, and even for the metal of pure Ni, hydrostatic measurements at room temperature are limited to 98 GPa [1] so far, and at high temperatures, there are only two sets of experiments with pulsed laser to heat the specimen, resulting in the uncertainty of the temperature measurements being as large as ~ 100 K above 1000 K [1,2], and the results of the two experiments are distinctly different.

According to the ensemble theory, the EOS and all the other thermodynamic properties of materials can be theoretically calculated from the partition function (PF) without any experimental measurement. However, this fundamental approach faces formidable computational challenges: for an N-atom system, a 3N-dimensional configuration integral is required to be solved, rendering it computationally intractable for realistic solids containing ~10²³ atoms. In spite of a century of research, this fundamental limitation persists in contemporary materials science.

The quasi-harmonic approximation (QHA) [3,4,5] has emerged as an effective compromise, achieving wide applications by simplifying lattice dynamics through harmonic vibrational modes. However, its inherent assumption of small atomic displacements becomes inadequate at elevated temperatures where anharmonic effects dominate. In fact, a recent study [6] has quantitatively demonstrated that the accuracy of QHA deteriorates progressively with both increasing temperature and atomic volume. Alternative approaches have also been developed to circumvent this limitation. The particle-in-a-cell model (PCM) [7,8,9] reduces the 3N- dimensional integral into a 3-fold one. Nevertheless, achieving the 3- dimensional integral at ab initio level remains a technical challenge. Moreover, PCM's validity is restricted to the conditions above the Debye temperature [10]. A more radical simplification appears in the classical mean field (CMF) [11,12,13,14] method, which reduces the 3-fold integral in PCM into a one-fold spherical approximation. This approach fundamentally assumes isotropic atomic interactions, a significant oversimplification that introduces substantial errors for most crystalline metals. The anisotropic nature of interatomic forces in crystalline systems, particularly along different crystallographic directions, renders CMF's spherical symmetry assumption physically inappropriate for accurate thermodynamic predictions.

In contrast to the methods mentioned above, a recently developed direct integral approach (DIA) [15] introduces a novel perspective by transforming the 3N-fold integration in PF into a 3N-dimensional effective volume, enabling the PF to be solved with ultra-high efficiency at ab initio level. Up to now, DIA has been successfully employed to calculate the EOS of several metals [15,16,17,18], argon [15], 2-D materials [19], and its success in reproducing the phase transitions of vanadium [20], zirconium [21] and aluminum [22] further highlights the accuracy of DIA. DIA without any empirical parameters works universally to any solids.

In the present work, DIA is applied to calculate the EOS of metal Ni up to 3000 K and simultaneously 500 GPa to further check the reliability of DIA and give theoretical predictions for the EOS at extreme conditions. In order to obtain an analytical EOS of the metal Ni, we applied the widely used Mie-Grüneisen-Debye (MGD) formalism [23,24,25] to fit our calculated results, which shows that one of the three parameters in the MGD model, the Debye temperature, could be removed without losing any accuracy of the fitting. Furthermore, a much more simplified and universal model is developed to produce analytical EOS for general solids.

2. Calculation Method

From the ensemble theory it is known that for a system containing N atoms occupying the spatial position r^N = {r₁, r₂, …, r_N}, the PF reads

(1)

in which U(r^N) is the interaction energy of the system, and m, h, k_B is the atomic mass, Planck constant, and Boltzmann constant, respectively. If the configurational integral

(2)

is solved, the pressure (P) as a function of V and T, i.e. the isothermal EOS could be obtained by

(3)

and all the other thermodynamic quantities can also be achieved without any experimental supporting.

To solve the 3N-dimentional integral in Eq. (2), an approach was developed recently by Ning et al. [15], and its essential framework is as follows.

For a crystal in an equilibrium state with the N atoms located at their ideal lattice sites R^N = {R₁, R₂, …, R_N} and the total energy of U₀(R^N), a transformation is firstly introaduced,

(4)

so that r'^N represents the displacements of atoms away from their equilibrium positions. Then, the Q in Eq. (2) is expressed as

(5)

Here, we define a function

(6)

and then the 3N-fold integral in Eq. (5) turns into

(7)

in which

(8)

is defined as the effective length, and r'_ix (r'_iy or r'_iz) is the x (y or z) coordinate of the ith atom away from its ideal lattice site with simultaneously the other two degrees of freedom and all the other atoms fixed at the ideal lattice sites.

The Eq. (6) could be expanded by a Taylor series as

(9)

Clearly, if U^'(r^'N) changes slow enough with respect to r^'N, ∂²U^'(r^'N)/∂r^'_ij∂r^'_ik and further F(r^'N) would approach to zero, so that the Q of Eq. (7) is close to

(10)

To ensure Eq. (10) is a good approximation, the axes of a Cartesian coordinate system should be set as the direction along which U^'(r^'N) changes slowest.

For face-centered cubic metal Ni, if the lattice orientations [100], [010] and [001] are selected as the axes of a Cartesian system, the variations of U^'(r^'N) along all the three axes are the slowest. Thus, L_x equals to L_y and L_z, so that

(11)

The PF calculations for Ni at specified atomic volumes is implemented through the following protocol. As illustrated in Figure 1, we constructed a perfect 3 × 3 × 3 supercell containing 27 atoms with an arbitrary volume and displaced an arbitrary atom along the [100] direction in increments of 0.05 Å until the U^'(x^') in Eq. (4) exceeded 3.5 eV, ensuring convergence of the effective length at 3000 K (<10^-5). By varying volumes, PF as a function of V was obtained from Eqs. (1, 8, 11) and then P-V relationship could be derived from Eq. (3).

First-principles calculations were performed using the Vienna Ab initio Simulation Package (VASP) [26,27] to get the U^'(x^'), and the computational parameters were set as follows: (1) the projector augmented wave (PAW) pseudopotential [28,29] is used for describing the electron-ion interactions with 3d⁹4s¹ selected as the valence orbitals; (2) the Perdew-Burke-Ernzerhof (PBE) [30] exchange-correlation functional within the generalized gradient approximation (GGA) framework is employed to describe the electron-electron interactions; (3) electronic self-consistency criterion is set as 10^-6 eV; (4) a Γ-centered 15 × 15 × 15 k-mesh generated via the Monkhorst-Pack method [31] is adopted to sample the Brillouin zone; (5) the tetrahedron method with Blöchl corrections [28] is used to determine the electron orbital partial occupancy with the plane-wave cutoff energy set as 400 eV. The convergence of all these parameters was verified across all the considered volumes with effective length fluctuating below 10^-5 Å.

3. Results and Discussion

3.1. The Room-Temperature Isotherm

Based on DIA, the room-temperature P-V relationship of Ni up to 500 GPa is obtained and shown as the red solid line in Figure 2, which includes three sets of experimental data [1,2,32] and the results calculated by QHA [33]. It is roughly seen that the calculations of both QHA and DIA agree well with the experiments. To specifically inspect the feasibility of DIA, under the same volume, differences between the pressure calculated by DIA (P_DIA), or QHA (P_QHA) and the experimental measurement (P_EXP) are conducted and presented in Figure 3.

As shown in Figure 3(a), most of P_DIA agree well with the P_EXP of Ref. [2] up to 61 GPa with the difference less than 1 GPa, but at P_EXP of 13.1, 17.6, 21.4, and 22.1 GPa (in the green box in Figure 3(a)), P_DIA presents a large relative deviation (|P_DIA – P_EXP| / P_EXP) of 25.2, 13.2, 11.8 and 25.1%, respectively, which should be attributed to the accidental error in the experiments. Comparatively, QHA manifests a little worse than DIA with nearly all the P_QHA smaller than the P_EXP, resulting in an average difference (|P_EXP - P_QHA|) of 2.85 GPa. In 2015, Pigott et al. [1] performed a larger range of diamond anvil cells compressions of Ni up to 98 GPa. Compared with the experiment of Ref. [2], higher precision is taken in Ref. [1] in terms of zero uncertainty for the controlled room temperature and better hydrostatic conditions. As shown in Figure 3(b), the difference between P_DIA and the P_EXP of Ref. [1] fluctuates around zero in the whole compression condition with average |P_EXP – P_DIA| of 1.35 GPa, while nearly all the P_QHA are larger than the P_EXP with an average deviation of 2.59 GPa.

Recently, Hirao et al. [32] reported a set of volume measurements of Ni up to 368 GPa, but only the ones below 21 GPa are hydrostatic. As shown in Figure 3(c), P_DIA agrees well with these hydrostatic measurements with an average difference of only 0.56 GPa, but most of the non-hydrostatic P_EXP are smaller than P_DIA, which qualitatively indicates the reliability of DIA since the pressures under non-hydrostatic conditions are recognized to be smaller than that under hydrostatic conditions. As the pressure calculated by QHA in Ref. [33] is just up to 180 GPa, comparisons between P_QHA and the experiment of Ref. [32] is only within 0 - 180 GPa. Obviously, nearly all the P_EXP are larger than P_QHA.

The above discussions show that when compared to the experiments, DIA manifests better than QHA, but from the perspective of practical applications, the differences of the room-temperature isotherms supplied by the hydrostatic experiments Ref. [2], Ref. [1] and calculations of QHA and DIA could be neglected.

3.2. The High-Temperature Isotherms

On the P-V relationship of Ni at high temperatures (> 300 K), Campbell et al. [2] using a multi-anvil press and laser-heated diamond anvil cells (LHDAC) completed a group of measurements with the covered P, T reaching 2457 K and 66 GPa. Pigott et al. [1] using LHDAC achieved the volume measurements of Ni up to 109 GPa and 2941 K. Unfortunately, a large uncertainty of the manipulated temperature is noted in both the experiments of Ref. [2] and Ref. [1], so that the measured P-V-T data are isolated points instead of isotherms. To inspect the feasibility of DIA, we calculated the pressure at each pair of the experimentally measured (V, T) of the Ni sample in Ref. [2] and Ref. [1], and then performed specific comparisons with the experiments. To show the differences clearly, we arranged the experimental pressures in Ref. [2] and Ref. [1] in an ascending order mapped by the Arabic numerals shown as the blue symbols in Figure 4(a) and (b), where the red symbols are the corresponding calculated results by DIA. Error bars for our calculated pressures are obtained as follows. Taking the point indicated by the arrow in Figure 4(a) as an example, the experimentally measured V and T corresponding to this pressure is 9.806 ± 0.010 cm³/mol and 1976 ± 92 K, and based on DIA, at 1976 K, when the volume of Ni is 9.806 cm³/mol, our calculated pressure is 40.384 GPa. Making the calculated pressure 41.435 GPa at volume of (9.806 – 0.010) cm³/mol and temperature of (1976 + 92) K as the upper bound and 39.334 GPa at volume of (9.806 + 0.010) cm³/mol and temperature of (1976 – 92) K as the lower bound, the upper error for this calculated pressure (40.384 GPa) is ΔP₊ = 1.051 GPa and the lower one is ΔP_- = 1.050 GPa.

From Figure 4(a) and (b) as well as the respective insets, we can see that our calculated results agree well with the experiments of Ref. [2] and Ref. [1], especially better with the latter one. Considering the uncertainty, if the calculated pressure has an overlap with the experimentally measured result, we think they are consistent. By this way, the degree of the coincidence of our calculations with the experiment of Ref. [2] is 59.55%, and with Ref. [1] reaches 78.28%. The re-arranged experimental data of Ref. [2] and Ref. [1] along with the corresponding calculated pressures by DIA are listed in supplementary material, Table SI and SII, respectively.

For realistic applications, having a good command of the isotherms of Ni at high temperatures is demanded. In Ref. [2] and Ref. [1], the measured P-V-T data were fitted by the Mie-Grüneisen-Debye formalism (MGD) [23,24,25], and from the fitted parameters, the isotherms of Ni could be obtained. However, the results in Ref. [2] and Ref. [1] are effective just to 66 GPa, 2457 K, and 109 GPa, 2941 K, respectively. In the present work, we calculated the isotherms of Ni up to 3000 K and simultaneously 500 GPa, which are listed in Table 1. To discuss the differences among the high-temperature isotherms of Ni reported by different experiments and theoretical calculations, the P-V curves of Ni under 1000, 2000, and 3000 K from Ref. [2], Ref. [1], DIA and QHA are displayed in Figure 5(a) - (c). Since the largest calculated pressure in Ref. [33] is limited to 180 GPa, and according to Ref. [34], solid Ni exists above 7.5 GPa at 2000 K and 50.6 GPa at 3000 K, the comparisons and discussions in this section are just within these conditions. Under the same volume, the differences between P_QHA, P_EXP of Ref. [2], or Ref. [1] and P_DIA at 1000, 2000, and 3000 K are correspondingly displayed in Figure 5(d) - (f), from which we can see that P_QHA [33] is close to P_EXP of Ref. [1] in the whole pressure and temperature range with the largest difference of less than 4.3 GPa. Whereas, both P_QHA and P_EXP of Ref. [1] are consistently smaller than P_DIA. The deviation of P_EXP in Ref. [1] from P_DIA increases gradually with increasing pressure with the largest relative error ((P_EXP – P_DIA) / P_DIA) of 3.5%, 4.1%, and 4.6 % at 1000, 2000, and 3000 K, respectively, while P_QHA almost presents a constant relative deviation ((P_QHA – P_DIA) / P_DIA) of 5.6% within the entire discussed conditions. It is noted that the pressures in the experiment of Ref. [2] were also determined based on the MGD model, but the derived isotherms are far away from the ones reported in Ref. [1], and also distinguish from the calculated results of QHA and DIA. As shown by the gray solid line in Figure 5(d) - (f), the largest relative difference between P_EXP in Ref. [2] and P_DIA reaches 17.1%, 17.2%, and 15.6% at 1000, 2000, and 3000 K, respectively.

The above discussions as well as those included in the section A show that our calculations agree better with the experiment of Ref. [1], in which the accuracy of the measurements should be higher than that of Ref. [2]. In addition, it is noted that QHA is rarely applied to such high temperatures up to 3000 K due to its effectiveness only to small atomic vibrations, but it should be appropriate for the Ni system because as confirmed by one of our investigations [6] that QHA would work well for the system with atomic volume smaller than 22 ų/atom, and under the lower pressure limit for the solid Ni [34], its atomic volume calculated by DIA at 300, 1000, 2000, and 3000 K is 10.95, 11.31, 11.33, and 9.74 ų/atom, respectively. Thus, it is reasonable that the isotherms of Ni calculated by QHA distinguish little from that by DIA with a nearly constant relative error of only 5.6% in a large temperature range from 300 up to 3000 K. As the accuracy of DIA increases with decreasing volume of the material [19], the predictions for the isotherms of Ni up to 3000 K and simultaneously 500 GPa should be reliable.

3.3. The Analytical EOS

For practical application of metal Ni, it is necessary to obtain an analytical isothermal EOS based on the calculated results of DIA to determine the pressure for arbitrary volume V and temperature T. For this purpose, the MGD model [23,24,25] should be a good choice since it has been applied successfully to fit experimental data to obtain analytical EOS at high temperatures [1,2]. In this model, the total pressure P is considered to be the summation of thermal pressure Pth and the room-temperature one P_300K. Pth as a function of V and T is determined by three parameters via the following equations,

(12)

and

(13)

where R is the gas constant and n is the mole number of atoms confined within volume V. The Grüneisen parameter γ and Debye temperature θ are regarded to be only volume dependent via

(14)

and

(15)

in which the quantities with subscript 0 represent the corresponding values at 300 K and one atmospheric pressure. γ₀, q and θ₀ are the three parameters to be decided. The P_300K can be described by the third-order Birch-Murnaghan (BM) [35] EOS

(16)

with K₀ and K´₀ the bulk modulus and its pressure derivative at one atmospheric pressure.

Although K₀ and K₀´ could be easily obtained by fitting experimental data via Eq. (16), determining γ₀, q and θ₀ will cost too much computer hours because an integral equation (Eq. (13)) is involved in the MGD model. Actually, in previous applications of the model [1,2], the parameter θ₀ measured by other experiment was pre-introduced. For example, Ref. [1] employed MGD model fitting the experimental P-V-T data of metal Ni with θ₀ of 415 K [36], and γ₀ and q were then determined to be 1.98 and 1.3. In fact, the fitted results are insensitive to θ₀. When we varied the value of θ₀ from 1 to 500 K without changing the value of γ₀ and q, the obtained pressures deviate from the original ones with θ₀ of 415 K by no more than 0.45 GPa (Figure 6), which is significantly smaller than the uncertainty of the experimental measurements.

It is notable that the internal energy expressed by Eq. (13) is based on the harmonic approximation, and tends to be U_T = 3nRT when the temperature T is larger than the Debye temperature. It is a fact that the Debye temperature for most solids is in the range of 200 - 400 K, and therefore the internal energy at temperature above 300 K should be independent of volume and can be simplified as 3nRT. Accordingly, the expression of the thermal pressure Pth can be simplified as

(18)

Thus, the MGD model is reduced to a simplified one without parameter θ₀ and the integral equation Eq. (13). We apply this simplified model to calculate the pressure of metal Ni with parameter γ₀ and q taken as 1.98 and 1.3 determined in Ref. [1], the results deviate from the predictions under MGD model by no more than 0.42 GPa [Figure 6], which is significantly smaller than the uncertainty of the experimental measurements. Furthermore, we obtained an analytical EOS of metal Ni by fitting the calculated results (Table 1) of DIA via the simplified model, showing that the fitted pressure deviates from the original one by no more than 1 GPa, which is smaller than the precision of DIA calculations [Figure 7]. As shown in Table 2, the parameters of the EOS determined by DIA present some difference from the ones determined by fitting experimental data under MGD model reported in Ref. [2] or Ref. [1]. Relatively, our theoretical parameters are closer to the latter ones.

Applying the analytical EOS of metal Ni, we obtained a map of pressure distribution over T-V plan shown in Figure 8(a), which is similar to the empirical predications (Figure 8(b)) based on the MGD fitting experimental data [1]. The isobars in Figure 8(a) and (b) indicate clearly that the pressure depends mainly upon the volume instead of the temperature varied from 300 to 3000 K. For given pressure, the volume expends linearly with the temperature and the expansion rate gets slower and slower with increasing pressure. When the pressure is larger than 400 GPa, the volume keeps nearly unchanged even if the temperature increases from 300 to 3000 K. For arbitrary volume (or pressure) and temperature, the pressure (or the volume) can be easily read from this pressure map.

The difference between our theoretical predictions and the empirical ones take place mainly in the high-pressure region. As shown in Figure 8(c), there exists a T-V curve along which our theoretical prediction is the same as the ones of the empirical prediction, and the theoretical pressure is larger than the empirical one by 5 GPa when the pressure is 150 GPa. As the pressure increases to 500 GPa, the pressure difference gets larger and larger up to 35 GPa, which can be understood in consideration of the fact that the empirical predictions using MGD model is based on the experimental data below 109 GPa.

4. Conclusions

In this work, we applied DIA at ab initio level to study the EOS of metal Ni up to 3000 K and simultaneously 500 GPa, and proposed an analytical expression of the EOS. The theoretical results are in good agreement with available hydrostatic experiments, suggesting that the predications under extreme conditions will be reliable. The simplified model developed to produce analytical EOS of solid should be universal and will find its vast applications for obtaining analytical EOS of various solids by fitting limited experimental or theoretical data.