Mechanical and Thermal Conductivity Properties of Enhanced Phases in Mg-Zn-Zr System from First Principles

In this paper, the mechanical properties and minimum thermal conductivity of ZnZr, Zn2Zr, Zn2Zr3, and MgZn2 are calculated from first principles. The results show that the considered Zn-Zr intermetallic compounds are effective strengthening phases compared to MgZn2 based on the calculated elastic constants and polycrystalline bulk modulus B, shear modulus G, and Young’s modulus E. Meanwhile, the strong Zn-Zr ionic bondings in ZnZr, Zn2Zr, and Zn2Zr3 alloys lead to the characteristics of a higher modulus but lower ductility than the MgZn2 alloy. The minimum thermal conductivity of ZnZr, Zn2Zr, Zn2Zr3, and MgZn2 is 0.48, 0.67, 0.68, and 0.49 W m−1 K−1, respectively, indicating that the thermal conductivity of the Mg-Zn-Zr alloy could be improved as the precipitation of Zn atoms from the α-Mg matrix to form the considered Zn-Zr binary alloys. Based on the analysis of the directional dependence of the minimum thermal conductivity, the minimum thermal conductivity in the direction of [110] can be identified as a crucial short limit for the considered Zn-Zr intermetallic compounds in Mg-Zn-Zr alloys.


Introduction
Grain refinement is a metallurgical phenomenon that has been exploited in magnesium alloys to achieve desired microstructure and mechanical properties. Zirconium, a powerful grain refiner, has been widely used in magnesium alloys [1][2][3]. Magnesium-zinc-zirconium (ZK) alloys mainly refer to those containing zirconium or grain refined by zirconium, such as ZE41, ZK60, WE43, ML12, and ML10, as well as OS-1-3, and such like, and these commercial Mg-Zn-Zr alloys comprise the basis of the current magnesium alloy business. It was previously agreed that grain refinement of magnesium alloys by Zr was noticeable at low levels of soluble Zr [4], but the subsequently detailed examinations showed that both insoluble zirconium particles and zirconium dissolved in the melt played a role in grain refinement [5]. These conditions require that the magnesium alloys contain maximum soluble and undissolved Zr content. Meanwhile, a large amount of Zr content could lead to the formation of Zn-Zr intermetallic compounds in Mg-Zn-Zr alloys.
Based on the study of the composition and distribution of zirconium, Morozova and Mukhina [6] proposed that the highly dispersed particles of Zn 2 Zr 3 , ZnZr, and Zn 2 Zr intermetallic compounds were the determining factors for the nano-structural mechanism of strengthening in the Mg-Zn-Zr system. Li et al. [7] proposed that the Mg-5Zn-2Gd-0.4Zr alloy (wt.%) showed a significant hardening response during aging, thus forming different morphologies, including ZnZr and Zn 2 Zr phases. Although the ZnZr 2 and Zn 2 Zr 3 phases do not belong to the ground state of the Zn-Zr system, they were dynamically

Computational Details
All calculations in this work were performed by using the Vienna ab initio simulation package code (VASP) [18] within the generalized gradient approximation (GGA) [19] of Perdew-Burke-Ernzerhof (PBE) [20] exchange correlation density functional. The electron configuration treated 3s 2 as a valence state for Mg, 3d 10 4s 2 as a valence state for Zn, and 4s 2 4p 6 5s 2 4d 2 for Zr, respectively. Extensive convergence tests suggested that the cutoff energy of 400 eV was enough for all phases in the calculations. The special points sampling integration was used over the Brillouin zone with 7 × 7 × 3 and 8 × 8 × 8 k-points using the Gamma-centered Monkhorst-Pack method [21] for MgZn 2 and the Zn-Zr system (including ZnZr, Zn 2 Zr and Zn 2 Zr 3 ), respectively, in geometry optimization. The convergence criterion of the Hellmane Feynman force was 0.01 eV/Å for complete relaxation of the atomic positions within the maximum stress on the atom of 0.02 GPa. The electronic iterations convergence was 1.0 × 10 −5 eV for the total energy calculated together with first-order Methfessel-Paxton smearing with a width of 0.2 eV. Considering the unfilled electron of the 4d shell of the transition metal Zr, the spin polarization was considered in the calculation with the initial magnetic moment 3 µ B according to Hund rules.
For obtaining the equilibrium bulk modulus B 0 of the spin state at 0 K, the ground state energy E 0 as a function of the cell volume within the Birch-Murnaghan equation of states (EOS) [22] was applied. Meanwhile, in order to investigate chemical stability, the contribution of the lattice vibrations F vib to the total Helmholtz free energy (F = E 0 + F vib ) was evaluated (it is worth noting that the contribution of the thermal electrons is negligible compared to lattice vibrations, and is therefore ignored in the current work). For the sake of computational efficiency, the vibrational free energy was derived by using the Debye-Grüneisen [23] model as follows: where k B is the Boltzmann constant and n is the number of atoms per formula unit. The Debye temperature Θ was obtained as proposed by using [24]: where M is the molecular mass per primitive cell, and B 0 and σ are the static bulk modulus and Poisson ratio at the equilibrium geometry, respectively. The f (σ) function is: The elastic coefficients were determined by applying a set of given deformation with a finite value fitting the total energy of the crystal, as implemented by Mayer et al. [25]. In order to remain within the elastic limit of the selected phases, small strains up to ±2% at 0.5% interval were used. Figure 1 shows the energy-volume fittings for ZnZr, Zn 2 Zr, Zn 2 Zr 3 at both no-spin and spin states, as well as MgZn 2 at a nonmagnetic state. It is clear that for energy that was dependent on the volume of the two states, the spin state was more energetically stable than the no-spin state for Zn-Zr intermetallic compounds, especially for the Zn 2 Zr phase. Thus, in the following discussion, we focus on properties at the magnetic state. The optimized lattice constants at the ground state and the corresponding fitting B 0 for each Zn-Zr intermetallic compound are given in Table 1. As seen, the calculated lattice parameters of Zn-Zr intermetallic compounds and the MgZn 2 phase are in good agreement with the available calculated results. When comparing the equilibrium bulk moduli of these four intermetallic compounds, which could describe the stiffness of the crystal to the applied strain, it can be observed that the B 0 of the selected Zn-Zr compounds is interestingly larger than that of MgZn 2, which has generally been considered as the main strengthening phase in Mg-Zn alloys. We estimated the chemical stability based on the calculated Helmholtz free energy F (eV/atom) of Zn-Zr compounds. Generally speaking, these phases are thermodynamically stable due to negative Helmholtz free energy, which is considered to be a key factor for the alloys' synthesis and stabilization; the more negative it is, the more stable the structure. Furthermore, it has been shown from the Helmholtz free energy F that the thermodynamic stability sequence is Zn 2 Zr 3 > ZnZr > Zn 2 Zr > MgZn 2 , and Zn 2 Zr 3 is the most thermodynamically stable compound.  [8]. b From Reference [26].

Elastic Constants, Polycrystalline Moduli
The calculated elastic constants of the various phases are shown in Table 2. In reality, Zn-Zr intermetallic compounds exhibit much better resistance to deformation than MgZn2 due to larger elastic constants; not only regarding C11 and C33 under uniaxial stress along the x or z axes, respectively, but also other compression moduli (C12 and C13) and shear moduli (C44, and C66). At this point it may be clear that ZnZr, Zn2Zr, and Zn2Zr3 are effective strengthening phases in Mg-Zn-Zr alloys, aside from MgZn2. For hexagonal MgZn2 and tetragonal Zn2Zr3 crystals C11 = C22 ≠ C33, the difference between C11 (C22) and C33 indicates that the two crystals have relatively strong anisotropic elastic constants, resulting in the directional dependence of the moduli. Interestingly, for hexagonal MgZn2 and tetragonal Zn2Zr3 crystals, the values of C33 is larger than that of C11, implying that the chemical bonds in the direction of [001] are stronger than those along the direction of [100]. The stronger chemical bonds result in hard compressing under uniaxial stress along the z axes. Moreover, the relatively large difference between C11 and C33 of MgZn2 implies that there is greater anisotropy on the directional dependence of the moduli than Zn2Zr3.

Elastic Constants, Polycrystalline Moduli
The calculated elastic constants of the various phases are shown in Table 2. In reality, Zn-Zr intermetallic compounds exhibit much better resistance to deformation than MgZn 2 due to larger elastic constants; not only regarding C 11 and C 33 under uniaxial stress along the x or z axes, respectively, but also other compression moduli (C 12 and C 13 ) and shear moduli (C 44 , and C 66 ). At this point it may be clear that ZnZr, Zn 2 Zr, and Zn 2 Zr 3 are effective strengthening phases in Mg-Zn-Zr alloys, aside from MgZn 2 . For hexagonal MgZn 2 and tetragonal Zn 2 Zr 3 crystals C 11 = C 22 = C 33 , the difference between C 11 (C 22 ) and C 33 indicates that the two crystals have relatively strong anisotropic elastic constants, resulting in the directional dependence of the moduli. Interestingly, for hexagonal MgZn 2 and tetragonal Zn 2 Zr 3 crystals, the values of C 33 is larger than that of C 11 , implying that the chemical bonds in the direction of [001] are stronger than those along the direction of [100]. The stronger chemical bonds result in hard compressing under uniaxial stress along the z axes. Moreover, the relatively large difference between C 11 and C 33 of MgZn 2 implies that there is greater anisotropy on the directional dependence of the moduli than Zn 2 Zr 3 . Table 2. The calculated independent elastic constants (GPa) of Zn-Zr intermetallic compounds and MgZn 2 using the strain-energy method with other calculated (Cal.) and experimental (Exp.) data.

Species
Reference In order to synthetically estimate the mechanical properties, the polycrystalline bulk modulus B, shear modulus G, and Young's modulus E were calculated via Voigt-Reuss-Hill approximations [30][31][32]. Figure 2 summarizes the calculated mechanical performance parameters. Notably, the B values of ZnZr, Zn 2 Zr, and Zn 2 Zr 3 are close, and are all larger than that of MgZn 2 . The larger B values responded to the stronger capacity of the resist deformation, reflecting the good resistance of the selected Zn-Zr intermetallic compounds to deformation. Meanwhile, the B values of all considered intermetallic compounds are in good agreement with the fitting bulk modulus, B 0 . The shear modulus G of the system in descending order is: Zn 2 Zr > ZnZr > Zn 2 Zr 3 > MgZn 2 . It is clear that the Young's modulus E has the same order as the shear modulus G, suggesting that Young's modulus E of the considered polycrystalline materials is more sensitive to the shear modulus than the bulk modulus. Relatively larger mechanical parameters, including bulk modulus B, shear modulus G, and Young's modulus E, prove that Zn 2 Zr has outstanding mechanical properties and pronounced strengthening effects among all strengthening phases. In contrast, the lowest B/G value 1.76 reveals its brittle characteristics relative to other phases, although the material behaves as ductile when the B/G ratio >1.75 [33]. Analysis of Figure 2 thus allows us to conclude that ZnZr and Zn 2 Zr 3 serve to combine the natures of high strength and great ductility. In order to synthetically estimate the mechanical properties, the polycrystalline bulk modulus B, shear modulus G, and Young's modulus E were calculated via Voigt-Reuss-Hill approximations [30][31][32]. Figure 2 summarizes the calculated mechanical performance parameters. Notably, the B values of ZnZr, Zn2Zr, and Zn2Zr3 are close, and are all larger than that of MgZn2. The larger B values responded to the stronger capacity of the resist deformation, reflecting the good resistance of the selected Zn-Zr intermetallic compounds to deformation. Meanwhile, the B values of all considered intermetallic compounds are in good agreement with the fitting bulk modulus, B0. The shear modulus G of the system in descending order is: Zn2Zr > ZnZr > Zn2Zr3 > MgZn2. It is clear that the Young's modulus E has the same order as the shear modulus G, suggesting that Young's modulus E of the considered polycrystalline materials is more sensitive to the shear modulus than the bulk modulus. Relatively larger mechanical parameters, including bulk modulus B, shear modulus G, and Young's modulus E, prove that Zn2Zr has outstanding mechanical properties and pronounced strengthening effects among all strengthening phases. In contrast, the lowest B/G value 1.76 reveals its brittle characteristics relative to other phases, although the material behaves as ductile when the B/G ratio >1.75 [33]. Analysis of Figure 2 thus allows us to conclude that ZnZr and Zn2Zr3 serve to combine the natures of high strength and great ductility.

Electronic Structures
To gain further insight into the reasons for the Zn-Zr system strengthening, the electron localization function (ELF) [34] was applied to assist in identifying the distribution of the charges and the bonding condition. As a general rule, the ELF value is on the range 0 ≤ ELF ≤ 1, where ELF = 0, 1 corresponds to the completely delocalized state and the perfect localization, respectively. The electron localization functions on the (110) plane for all selected intermetallic compounds are

Electronic Structures
To gain further insight into the reasons for the Zn-Zr system strengthening, the electron localization function (ELF) [34] was applied to assist in identifying the distribution of the charges and the bonding condition. As a general rule, the ELF value is on the range 0 ≤ ELF ≤ 1, where ELF = 0, 1 corresponds to the completely delocalized state and the perfect localization, respectively. The electron localization functions on the (110) plane for all selected intermetallic compounds are presented in Figure 3. Clearly, there are obvious ionic characteristics in ZnZr, Zn 2 Zr, and Zn 2 Zr 3 intermetallic compounds due to the delocalization around Zn and localization around Zr. Meanwhile, Zn-Zn and Zr-Zr present typical metal bond characteristics because of the even distribution of charges between the component atoms. In contrast, Mg-Zn bonds show covalent characteristics based on the apparent accumulation of charge distribution between Mg and Zn atoms in MgZn 2 alloys. This result is consistent with the investigation of Reference [35], where more hybridized peaks between Mg p and Zn p appear near the Fermi level, indicating the presence of strong covalent bonding. Based on the above discussion, the bonding characteristics in Zn-Zr and MgZn 2 intermetallic compounds play a role in determining a ductile or brittle nature, meaning that ionic bonds in Zn-Zr intermetallic compounds cause them to have lower ductility than MgZn 2 . In addition, for Zn-Zr intermetallic compounds, the strength of the ionic bond was also compared based on the result of charge transfer using the Bader charge analysis. For a reasonable and intuitive comparison, the average charge transfer amount of per Zr atom (e/atom) in Zn-Zr intermetallic compounds could be used as a basis, and the descending order is: Zn 2 Zr (1.01) > ZnZr (0.95) > Zn 2 Zr 3 (0.72). From the perspective of ionic bonds, the strength of the Zn 2 Zr phase is stronger than ZnZr and Zn 2 Zr 3 , which is entirely consistent with the above elastic moduli results.

Minimum Thermal Conductivity and Anisotropy
It is well-known that thermal conductivity is inversely proportional to temperature. At elevated temperatures, the thermal conductivity will decrease to a limit value considered as the minimum thermal conductivity, which can be developed to identify candidate materials for hightemperature applications [36,37]. For the purpose of precisely calculating the minimum thermal conductivity of selected Zn-Zr and MgZn2 intermetallic compounds with anisotropic chemical bonds, the modified Clarke relation by Liu et al. [16] was used as defined: (4) where kB is the Boltzmann's constant, vm is the average sound velocity, NA is Avogadro's number, ρ is the density, M is the molecular weight, and n is the number of atoms in the molecule. The average sound velocity vm is given by [38,39]:

Minimum Thermal Conductivity and Anisotropy
It is well-known that thermal conductivity is inversely proportional to temperature. At elevated temperatures, the thermal conductivity will decrease to a limit value considered as the minimum thermal conductivity, which can be developed to identify candidate materials for high-temperature applications [36,37]. For the purpose of precisely calculating the minimum thermal conductivity of selected Zn-Zr and MgZn 2 intermetallic compounds with anisotropic chemical bonds, the modified Clarke relation by Liu et al. [16] was used as defined: where k B is the Boltzmann's constant, v m is the average sound velocity, N A is Avogadro's number, ρ is the density, M is the molecular weight, and n is the number of atoms in the molecule. The average sound velocity v m is given by [38,39]: where B and G are the bulk modulus and shear modulus, respectively. Using Liu's model, the minimum thermal conductivity of ZnZr, Zn 2 Zr, Zn 2 Zr 3 , and MgZn 2 is 0.48, 0.67, 0.68, and 0.49 (W·m −1 ·K −1 ), respectively. Since the values of the minimum thermal conductivities of the Zn-Zr intermetallic compounds such as Zn 2 Zr and Zn 2 Zr 3 are far larger than MgZn 2 , it can be proved that the thermal conductivity of the Mg-Zn-Zr alloy will be markedly improved as the precipitation of Zn atoms from the α-Mg matrix form Zn-Zr intermetallic compounds other than MgZn 2 . Meanwhile, Zn 2 Zr 3 , with maximum thermal conductivity, can be considered as the most important contribution to the total thermal conductivity due to its minimum Helmholtz free energy, which is considered to be a key factor for the alloys to be formed. However, for ZnZr, the difference in the minimum thermal conductivity between them could still be ignored. Furthermore, an important question to ask is how the minimum thermal conductivity along a different direction can affect the overall minimum thermal conductivity. To clarify this point, the directional dependence of the minimum thermal conductivity can be computed from the quasi-transverse or quasi-longitudinal sound velocities and the number density of atoms per mole (n) of the compound, according to Cahill's model [17]: for tetragonal Zn 2 Zr 3 crystal structure symmetry, the acoustic velocities can simply be written as:   In addition, the theoretical minimum thermal conductivities along different principle directions can be obtained using the acoustic velocities, as shown in Figure 5. For all considered intermetallic compounds, the directional dependence of the minimum thermal conductivity obtained in the present calculation is around 30% higher than those obtained by the Liu relation. Nevertheless, by comparing the values of minimum thermal conductivity in different directions, we can still understand the anisotropy of the minimum thermal conductivity to some extent. As can be observed, the minimum thermal conductivities of the ZnZr, Zn2Zr, and MgZn2 intermetallic compounds are the same in the three principle axis directions, indicating that the limit of thermal conductivity along the x, y, and z axes are directionally insensitive. Compared with ZnZr, Zn2Zr, and MgZn2, the minimum thermal conductivity of Zn2Zr3 along the z axis is higher than the x and y axes. What is striking in this figure is the minimum thermal conductivity along the [110] direction. It can be seen that the value of minimum thermal conductivity along the [110] direction is the smallest for all considered directions, indicating that Zn-Zr ionic bonding in the [110] direction will greatly affect the thermal conductivity at high temperatures, which will result in slower heat dissipation in the [110] direction. The minimum thermal conductivity in the [110] direction will be a crucial factor when considering the minimum thermal conductivity for the selected Zn-Zr intermetallic compounds, especially for Zn2Zr3. This fact definitely indicates that the bonding anisotropy, reflected in the elastic constant anisotropy, leads to anisotropy in the high-temperature limit of thermal conductivity-that is, the minimum thermal conductivity. In addition, the theoretical minimum thermal conductivities along different principle directions can be obtained using the acoustic velocities, as shown in Figure 5. For all considered intermetallic compounds, the directional dependence of the minimum thermal conductivity obtained in the present calculation is around 30% higher than those obtained by the Liu relation. Nevertheless, by comparing the values of minimum thermal conductivity in different directions, we can still understand the anisotropy of the minimum thermal conductivity to some extent. As can be observed, the minimum thermal conductivities of the ZnZr, Zn 2 Zr, and MgZn 2 intermetallic compounds are the same in the three principle axis directions, indicating that the limit of thermal conductivity along the x, y, and z axes are directionally insensitive. Compared with ZnZr, Zn 2 Zr, and MgZn 2 , the minimum thermal conductivity of Zn 2 Zr 3 along the z axis is higher than the x and y axes. What is striking in this figure is the minimum thermal conductivity along the [110] direction. It can be seen that the value of minimum thermal conductivity along the [110] direction is the smallest for all considered directions, indicating that Zn-Zr ionic bonding in the [110] direction will greatly affect the thermal conductivity at high temperatures, which will result in slower heat dissipation in the [110] direction. The minimum thermal conductivity in the [110] direction will be a crucial factor when considering the minimum thermal conductivity for the selected Zn-Zr intermetallic compounds, especially for Zn 2 Zr 3 . This fact definitely indicates that the bonding anisotropy, reflected in the elastic constant anisotropy, leads to anisotropy in the high-temperature limit of thermal conductivity-that is, the minimum thermal conductivity.

Conclusions
In summary, the structural stability, mechanical properties, bonding characteristics, and minimum thermal conductivity for the intermetallic compounds ZnZr, Zn2Zr, Zn2Zr3, and MgZn2 have been investigated by first-principles calculations. Moreover, the crystal direction/minimum thermal conductivity relationship of the materials were also established.
Based on the difference between C11 and C33, MgZn2 was found to possess greater anisotropy on the directional dependence of the modulus than Zn2Zr3. ZnZr and Zn2Zr3 was found to combine the natures of high strength and great ductility. Strong ionic bonds in ZnZr, Zn2Zr, and Zn2Zr3 was found to lead to the characteristics of a higher modulus but lower ductility than MgZn2. Furthermore, the strength of the Zn2Zr phase was found to be stronger than ZnZr and Zn2Zr3 based on the maximum charge transfer using Bader's charge analysis. Based on the calculated minimum thermal conductivities of the intermetallic compounds ZnZr, Zn2Zr, Zn2Zr3, and MgZn2, we conclude that the thermal conductivity of the Mg-Zn-Zr alloy will be markedly improved as the precipitation of Zn atoms from the α-Mg matrix help to form Zn-Zr binary alloys. However, the minimum thermal conductivity along the [110] direction may serve to be a crucial limit.

Conclusions
In summary, the structural stability, mechanical properties, bonding characteristics, and minimum thermal conductivity for the intermetallic compounds ZnZr, Zn 2 Zr, Zn 2 Zr 3 , and MgZn 2 have been investigated by first-principles calculations. Moreover, the crystal direction/minimum thermal conductivity relationship of the materials were also established.
Based on the difference between C 11 and C 33 , MgZn 2 was found to possess greater anisotropy on the directional dependence of the modulus than Zn 2 Zr 3 . ZnZr and Zn 2 Zr 3 was found to combine the natures of high strength and great ductility. Strong ionic bonds in ZnZr, Zn 2 Zr, and Zn 2 Zr 3 was found to lead to the characteristics of a higher modulus but lower ductility than MgZn 2 . Furthermore, the strength of the Zn 2 Zr phase was found to be stronger than ZnZr and Zn 2 Zr 3 based on the maximum charge transfer using Bader's charge analysis. Based on the calculated minimum thermal conductivities of the intermetallic compounds ZnZr, Zn 2 Zr, Zn 2 Zr 3 , and MgZn 2 , we conclude that the thermal conductivity of the Mg-Zn-Zr alloy will be markedly improved as the precipitation of Zn atoms from the α-Mg matrix help to form Zn-Zr binary alloys. However, the minimum thermal conductivity along the [110] direction may serve to be a crucial limit.