A critical study of the elastic properties and stability of Heusler compounds: Phase change and tetragonal $X_{2}YZ$ compounds

In the present work, the elastic constants and derived properties of tetragonal and cubic Heusler compounds were calculated using the high accuracy of the full-potential linearized augmented plane wave (FPLAPW). To find the criteria required for an accurate calculation, the consequences of increasing the numbers of $k$-points and plane waves on the convergence of the calculated elastic constants were explored. Once accurate elastic constants were calculated, elastic anisotropies, sound velocities, Debye temperatures, malleability, and other measurable physical properties were determined for the studied systems. The elastic properties suggested metallic bonding with intermediate malleability, between brittle and ductile, for the studied Heusler compounds. To address the effect of off-stoichiometry on the mechanical properties, the virtual crystal approximation (VCA) was used to calculate the elastic constants. The results indicated that an extreme correlation exists between the anisotropy ratio and the stoichiometry of the Heusler compounds, especially in the case of Ni$_{2}$MnGa.


I. INTRODUCTION
Heusler-type intermetallic compounds X 2 Y Z (X, Y = transition metals, and Z = main group elements) have become of particular interest due to their fascinating thermal, electrical, magnetic, and transport properties 1 . The Heusler compounds crystallize in a face-centered cubic (fcc) lattice. They are distinguished in two groups: regular or inverse Heusler compounds. Regular Heusler compounds belong to F m3m (space group no. 225) symmetry, and inverse Heusler compounds belong to F 43m (space group no. 216). Both cubic phases may undergo a cubic-tetragonal phase transition, in which the regular Heusler compounds transform from F m3m to the tetragonal I4/mmm (no. 139), and the inverse Heusler compounds transform from F 43m to tetragonal I4m2 (no. 119). Thus, the parent cubic and obtained tetragonally distorted phases obey a supergroup-subgroup relation.
Due to a simple feature of Heusler compounds, it is critically important to have an instrument for phase prediction. For example, cubic ferromagnetic Heusler compounds follow the Slater-Pauling rule for localized moment systems. Their magnetic moment m depends simply on the valence electron concentration n v with m = n v − 24. Further, prospective candidates for su-perconductivity include certain Heusler compounds with 27 electrons that exhibit a saddle point at the L point close to E F in the band structure according to the van Hove scenario 2,3 . On the other hand, a high density of states at the Fermi level causes instability and a phase transition to lower symmetry forced by a band Jahn-Teller distortion 4,5 . This competition is one example that shows the importance of phase prediction in the Heusler compounds. However, both tetragonal and cubic phases have their own importance for industrial as well as fundamental research.
Tetragonally distorted Heusler compounds have attracted interest in the field of spintronics, in particular, for spin-torque applications, owing to their magnetic anisotropy in the perpendicular axes [6][7][8][9] . Therefore, the theoretical prediction of new materials with suitable design properties is active research in this field 10 .
In fact, processing and designing new materials requires knowledge of physical properties, such as hardness, elastic constants, melting point, and ductility. The calculation of elastic constants is an efficient and fast tool used for elucidating physical properties as well as the mechanical stability and possible phase transitions of crystalline systems. Applied strains, such as shear or elongation, provide not only valuable information about the instability itself but also the directional depen-dence of instabilities in crystals. The directional dependence of instabilities becomes important when not only cubic-tetragonal but also cubic-hexagonal, tetragonalhexagonal, and lower symmetry phase transitions are relevant, such as those observed in Mn 3 Ga 6,7 . Unlike mechanical instability, the determination of elastic constants is essential for applications of magnetic shape memory alloys. The elastic constants also provide valuable information on the structural stability, anisotropic character, and chemical bonding of a system [11][12][13] . Moreover, other measurable properties can be estimated using the elastic constants, such as the velocity of sound, Debye temperature, melting point, and hardness. This information is an essential requirement for both industrial applications and fundamental research. For example, these properties are essential for studying superconductivity and heavy fermion systems in which a drastic change of elastic constants and related properties have been reported upon the phase transition 14,15 . The elastic properties are so important that Gilman 16 concluded: "the most important properties of a crystal are its elastic constants".
In the present study, some well-known tetragonal and cubic Heusler compounds are examined and compared with the available experimental and theoretical data [17][18][19][20][21][22][23][24][25] . Starting from the cubic phase, cubic L2 1 Ni 2 MnGa and Rh 2 FeSn are considered for detailed studies. For the tetragonally distorted systems, Ni 2 MnGa (in the nonmodulated tetragonal (c/a > 1) structure), Mn 2 NiGa, Fe 2 MnGa, and Mn 2 FeGa Heusler compounds are examined. The intermetallics Mn 2 Y Ga (Y = Fe, Ni) and X 2 MnGa (X = Fe, Ni) undergo tetragonal magnetostructural transitions that may result in half metallicity and magnetic shape memory or magnetocaloric effects. In the case of Ni 2 MnGa, the composition dependence (chemical disorder effect) of the phase transition is studied using the virtual crystal approximation (VCA). Calculating the mechanical and elastic properties of off-stoichiometric compounds in the tetragonal phases elucidates phase transformations. Elastic constants and mechanical properties of some Rh-based Heusler compounds, reported by Suits 26 , are calculated. The dependence of the elastic constants and the number of used k-points and plane waves (defined in full-potential linearized augmented plane wave (FPLAPW) by R MT k max , where R MT is the muffin tin radius and k max is the largest k vector) are discussed in detail. The importance of using sufficiently large numbers of k-points and plane waves for a reliable estimation of the elastic properties is demonstrated.
The present work concentrates on the elastic properties of metastable cubic, tetragonal, and phase change materials that exhibit magnetic shape memory or magnetocaloric effects. The results for the famous series of half-metallic Co 2 -based Heusler compounds that has a high impact on magnetoelectronics will be published elsewhere 27 .

A. Computational details
In this section, the basic equations for calculating the elastic constants are presented (for more details see Appendix A). The most easily determined quantity is the bulk modulus B, which provides the behavior of the crystal volume or lattice parameters under hydrostatic pressure. There are several ways to calculate the bulk modulus from the energy-volume E tot (V ) relation (see References [28][29][30]. In the present work, the bulk modulus B is determined by fitting the total energy calculations to the Birch-Murnaghan equation of state 29,30 . According to this model, the dependence of the energy on the change in the crystal volume V under hydrostatic pressure p is given by where B ′ = dB/dp is the pressure derivative of the bulk modulus and η = v 2 − 1, with the ratio v = (V 0 /V ) 1/3 of the actual volume V under pressure to the relaxed volume V 0 at the lowest total energy E 0 (see Reference 31). The related pressure is given by For tetragonal crystals, the dependence of the bulk modulus B on the elastic stiffness is given by B = 1 9 (2c 11 + 2c 12 + 4c 13 + c 33 ).
In the case of cubic crystals, c 12 = c 13 , c 11 = c 33 , and c 44 = c 66 (see also Appendix A). Therefore, the equations of the elastic constants for cubic systems are easily obtained from the tetragonal equations.
The remaining elastic properties are determined by applying different types of strain (e i ) to the tetragonal lattice and by applying proper relations between the total energy and the strain components. The energy E(e i ) of the strained lattice is calculated using Hooke's law (see equation (A3)). According to Wallace 31 , if the strains e i are small, the change of the energy is given by Here, ∆E = E(e i ) − E 0 with equilibrium energy E 0 at volume V 0 without strain. The linear terms vanish at equilibrium or if the strain causes no change in the volume of the crystal. The elastic constants c i,j are obtained from the second-order terms and are calculated from the second derivatives of the energy with respect to the strains: The second derivative relation of elastic constants (c ij ) with total energy highlights the importance of an accurate calculation of the total energy. Therefore, the choice of the density functional theory solver in the calculation of elastic constants and related properties is significant.
For cubic crystals, three independent elastic constants need two strains for calculations, while in tetragonal systems (space group Nos. 89-142), six independent elastic constants need five different strains. In fact, there are numerous ways to apply the six different strains and their combinations to the crystal. The necessary side condition of equation (5) is that the volume must be conserved when applying the strain. Therefore, the use of linear strain components (e i = δ in all possible cases and combinations in the strain matrix of equation A2 of Appendix A) would lead to large uncertainties because they are not always volume conservative, or they make the use of additional derivatives necessary (for example, ∂E ∂V , ∂V ∂ei , and higher orders). Table I and Figure 1 summarize the applied strains that are used to determine the elastic constants in the present work. The applied strains are the same as those reported by Kart et al 32 . The isotropic strain (0) is not used directly for the calculation of the elastic constants, as it gives the same information as discussed for the bulk modulus B (see discussion above). The five strain types (Equations (1)- (5) in Table I) are chosen to be volume conservative. The last strain type does not conserve the volume, but it keeps the same symmetry as the crystal and thus can be calculated from the energy versus c/a relation 32 , where c/a is the ratio of the two independent lattice parameters of the tetragonal crystals.
In the present work, six distortions of each type in the range of −3% ≤ δ ≤ +3% were applied to the relaxed structure with V 0 from the structural optimization using the Birch-Murnaghan equation of state. For tetragonal systems, the energy E(δ) versus applied strain curves were fitted to a fourth-order polynomial E(δ) = E 0 + a 2 δ 2 + a 3 δ 3 + a 4 δ 4 .
Here, an additional method of verifying the convergence as well as the accuracy of the results is introduced. In principle, it is sufficient to use either Equation (4) or (5) of Table I to calculate all six elastic constants. However, the elastic constants and all related properties are calculated with both equations to ensure that they provide the same results. This happens, indeed, only if the results are well converged (see also Section III A). The system is overdetermined by using both types of strains; however, in this way, the accuracy of the calculated quantities can be estimated. In fact, the values reported here have an error below 0.5%. The combination of different strains as well as different types of equations of state  Table I. allow determination of the uncertainty of the calculated results, which is expected for a reliable computer experiment.

B. Electronic structure calculations
The ab-initio electronic structure calculations were performed using the Wien2k code 33 . The all-electron full-potential method, FLAPW, with an unbiased basis covers all elements of the periodic table with any spin configuration. This feature is essential for Heusler compounds because they may contain diverse types of atoms, including lanthanide and actinide atoms, together with exotic magnetic ordering. The accuracy of this method makes it suitable for the studied systems. For example, Co 2 TiAl fails with spherical potentials or full symmetry potentials together with bare exchange-correlation functionals neglecting gradient corrections [34][35][36] . Since elastic constants are calculated from the second derivatives of the total energy, an accurate calculation of total energy is extremely important.
The exchange-correlation functional was taken in the generalized gradient approximation of Perdev, Burke, and Enzerhoff (GGA-PBE) 37,38 . The number of plane waves was restricted by R MT k max = 9, and the number of k-points was set to 8000 k-points in the full Brillouin zone. As discussed in Section III A, these criteria ensure the convergence of the calculated elastic properties for the investigated systems.
The lattice parameters were optimized before calculating the elastic constants. The results of the structural optimizations are summarized in Table II along with some  (5) are volume conservative. Only components with ei = 0 are given. Type (0) corresponds to calculation of the bulk modulus.

Type
Strain previously reported experimental and theoretical values.
Here, the c/a ratio was obtained by a full optimization of the Heusler compounds in tetragonal space groups 119 or 139. In other words, to find the energy minimum, not only the c/a ratio changed (the elongation of c) 32,39 (see also Figure 2) but also the volume of the structures was relaxed. The assignment of lattice parameters should be performed carefully since the optimization was performed in the tetragonal symmetry. When reducing the cubic f cc cell to a tetragonal f ct cell, the cubic lattice parameter a = a c becomes c, and the tetragonal parameter a = a t becomes a c /sqrt (2). To better understand the distortion of the cubic Heusler structure to the tetragonal Heusler structure, the distortion parameter ε is defined by where a t is the tetragonal a parameter. At ε = 0, the structure is cubic; at ε < 0, the cubic cell is compressed; and at ε > 0, the cubic cell is elongated along one of the principal axes.
The magnetic state was verified using different settings of the initial magnetization: ferromagnetic (all initial spins parallel) or ferrimagnetic (initial spins partially antiparallel). From the studied systems, all Mn 2 YZ compounds exhibit ferrimagnetic order in which one Mn has majority and the other Mn has minority orientation. For all other compounds, the ferromagnetic ground state has the lowest total energy.

III. RESULTS
Prior to discussing the tetragonal Heusler compounds, the elastic constants of the unstable and metastable cubic systems are discussed. Here, Ni 2 MnGa with L2 1 structure is selected as a metastable system. This compound is one of the most investigated materials owing to its shape memory behavior and its potential applications in actuator devices. In fact, in Section III A, only the cubic phase is addressed, and the tetragonal phase of Ni 2 MnGa is discussed in Section III B. Moreover, Ni 2 MnGa was used as an example case to study the significance of increasing the number of k-points and plane waves and their relations to the convergence of the elastic constants. In the second part of this section, the tetragonal phase of Heusler compounds are discussed. The elastic constants together with the corresponding measurable properties for selected tetragonal Heusler compounds are investigated. The role of the stoichiometry on the phase transition of Ni 2 MnGa is also explored.
A. Elastic constants and metastability in cubic and tetragonal compounds.
Stoichiometric Ni 2 MnGa undergoes a structural phase transition from the austenite into the martensite phase 42 . Depending mostly on the composition, the martensite structure is characterized by the tetragonal 5M modulated structure with c/a ≈ 0.94, the orthorhombic 7M structure with c/a ≈ 0.9, and the non-modulated tetragonal structure with c/a ≈ 1.2 [39]. Figure 2 shows the appearance of different stable and metastable phases with varying c/a elongation. To focus on the cubic phase, Figure 2(b) shows only the small range of strains with c/a < 1 that cover the cubic phase. The deepest energy minimum is located at a strain of about ε = 0.27, corresponding to c/a ≈ 1.26 (non-modulated phase). A shallow minimum (see Figure 2(b)) appears at a strain of about -0.05 (c/a ≈ 0.94, 5M phase). The elastic constants of the structure with c/a > 1 will be discussed in Section III B. As shown in Figure 2(c), the metastable cubic phase exists only under an infinitesimal strain and exhibits a very low energy modulation. The optimized lattice constant of cubic Ni 2 MnGa with the L2 1 structure is 5.81Å, in excellent agreement with the experimental value of 5.82Å [41]. This phase is only stable within ±1 meV energy changes, which confines the lattice distortion to < ±1%. This distortion relates to the tetragonal distortion, providing the c 11 − c 12 combination of elastic constants. To observe such a non-trivial change in energy and to have a smooth dependence on strain, the results need to be precisely converged with very high precision. This does not imply, however, that the results do not need to be converged for a wide energy window with a deep minimum. As an example, the importance of converged results is demonstrated in Figure 3. The c 44 shear modulus is stable for a lattice distortion of about ±3% and an energy change of more than ±30 meV. In this case, the rough calculation provides a smooth curve, but the calculated elastic constants significantly deviate from the converged results. The convergence of the results have the same importance for the calculation of the bulk modulus (B) (see Figure 3(c,d)).   Figure 3 shows the convergence of c 44 and B with respect to R MT k max and number of k-points. As shown in (a,c) for constant R MT k max = 7, increasing the number of k-points increases the c 44 and B values. These results converge at 8000 k-points. In contrast, increasing R MT k max decreases c 44 and B. Note that a similar result could be obtained at a more relaxed criterion, for example at 2000 k-points and R MT k max = 7, due to error cancellations. Therefore, R MT k max = 9 and 8000 k-points are the minimum criteria to converge the results for the systems studied in this work.
The calculated elastic constants of Ni 2 MnGa with L2 1 structure are given in Table III. The calculated results show a reasonable agreement with the experimental results and coincide with previously reported theoretical  44 . Moreover, the employed experimental method may result in different measured elastic constants 43 . As an example, the c 44 value deviates by about 60 GPa based on the experimental method. In general, the measured elastic constants are inversely related to temperature 45 . Hence, a higher value should be expected for the calculations. Fortunately, c ′ , which is a difference between two constants (c 11 and c 12 ), is argued to be less dependent on temperature 45 . In fact, the calculated c ′ = 4.5 GPa exhibits a better agreement with the experiment (4.1 GPa) than previously reported values 32,39 . In addition, c 12 and c 11 do not have any physical basis; in other words, no phonon mode directly corresponds to these constants. Mixing with other stiffness (c ij ), however, results in a meaningful combination. For example, the tetragonal shear modulus c ′ = (c 11 − c 12 )/2 corresponds to cubic-tetragonal distortion. Moreover, it is well established that c ′ -associated with slow transverse acoustic waves 43 -plays an important role in the occurrence of structural transformations. Another important quantity is the Cauchy pressure c p = c 12 − c 44 . A negative value of c p (c 12 < c 14 ) may indicate covalent bonds, where the angular dependence of the inter-atomic forces becomes important.
Furthermore, detailed analysis of elastic constants sheds light on the stability and phase transition in Heusler compounds. Cubic Ni 2 MnGa with soft c ′ is on the border of the phase transition. With the similar interpenetration (see Appendix B), the large elastic anisotropy A e of Ni 2 MnGa (A e =28) hints on its tendency to deviate from the cubic structure. Anisotropy is another indicator for the instability of cubic structures. The elastic anisotropy of crystals is also an important parameter for engineering since it correlates to the possibility of micro-cracks in materials. Unlike the mechanical properties, anisotropy shows the tendency of a sys-tem toward phase transitions as it inversely relates to the c ′ parameter. In fact, an illustrative way to show the anisotropy is to visualize the rigidity modulus G(r) or Young's modulus E(r). Figure 2(a) shows that the rigidity modulus is largest in the 111 -type direction that is along the tetragonal axes. Such a significant deviation from spherical shape indicates that the moduli of Ni 2 MnGa exhibit a large degree of anisotropy. In principle, when c 11 −c 12 → 0, the rigidity distribution exhibits a stronger directional dependency, as shown for Ni 2 MnGa in Figure 2(a). In the next step, a compound that is not stable in the cubic structure was examined for comparison. In the case of Rh 2 FeSn shown in Figure 4, the cubic structure exhibits a maximum of the total energy, and any tetragonal distortion will lead to a different stable structure. Based on Figure 4, the appearance of energy minima are expected for two different tetragonally distorted systems with c/a > 1 and c/a < 1. Larger distortions show that c/a > 1 is the stable phase, while c/a < 1 is a metastable phase. Previous works only reported the structure with ε > 0 26,46 . However, the calculated elastic constants also supported the instability of the cubic phase from negative values of tetragonal shear modulus c ′ and anisotropy A e . The metastable structure with ε < 0 appears when expanding the in-plane lattice parameter a while keeping c at the cubic lattice parameter. Such a situation may be artificially initialized by epitaxial or pseudomorphic thin film growth on a substrate with an appropriate lattice parameter. Similar metastable situations may exist in many other tetragonal Heusler compounds and will open the field of lattice parameter engineering to enlarge the number of properties on demand.

B. Tetragonal Heusler compounds
The tetragonal Heusler compounds studied in this work along with their elastic constants are summarized in Tables IV and V. The Heusler intermetallics Mn 2 Y Ga (Y =Fe, Ni) and X 2 MnGa (X =Fe, Ni) undergo tetragonal magneto-structural transitions that result in halfmetallicity, magnetic shape memory, or magneto-electric effects. In this section, the off-stoichiometric compositions are briefly discussed, and then, the elastic constants and related properties of the Heusler compounds are analyzed. Calculating the elastic properties of the tetragonal phases illuminates the structural transformations, chemical bonding, and mechanical stability of these intermetallic compounds for applications. Likewise, the elastic properties of Rh-based Heusler compounds synthesized by Suits 26 are calculated. Although Ni 2 MnGa has been widely studied experimentally at different phases, there is no experimental measurement of the elastic modulus of the non-modulated tetragonal phase, and only theoretical works on this phase have already been reported 32,39 .
As shown in Table IV, in the case of Ni 2 MnGa, the results show a qualitative agreement with the previous theoretical report (values in brackets). However, a quantitative comparison of the results reveals some significant deviations. These differences can be traced back to the calculation method and the method of performing structural optimization. To address this problem, the cubic phase is briefly considered. In the cubic phase, the present as well as other calculations are performed for the same lattice parameters using different calculation schemes. Here, the results of FPLAPW calculations are slightly larger than those of projected augmented wave (PAW) calculations, and the deviations range from 1% for c 11 up to 7% for c 44 . These small differences are expected because of the selected methods and convergence criteria. In contrast, in the case of the tetragonal system (see Table IV), the deviations between calculations range up to 30%, such as the case of c 44 and c 66 . Indeed, the different calculation methods should not lead to such a large discrepancy (if all factors are set carefully), and these observed differences mainly arise from the underlying structural optimization. Here, all initial tetragonal structures are fully optimized at their relevant symmetries (I4/mmm or I4m2), and thus, they differ from the simply elongated cubic structures (see also Section II B for more details about the calculations).
The elastic constants of all studied tetragonal Heusler compounds follow the inequality B > c 44 > G > c ′ > 0, so that the tetragonal shear modulus c ′ is the main constraint on the stability and properties. Pugh's and Poisson's ratio (see Appendix A) supply valuable information about the malleability and the type of bonding in crystals. Small values of Pugh's ratio indicate low malleabil- ity of crystals 47 , meaning they are brittle. Pugh's ratio indicates the type of bonding, namely covalent or metallic bonding, because changes in the angle of a covalent bond require more energy than stretching-stressing the bond, meaning G becomes larger compared to B; this leads to smaller values of k. On the other hand, Poisson's ratio for covalent bonding is about ν = 0.1 and increases for ionic crystals up to ν = 0.25. For instance, Pugh's and Poisson's ratios of strong covalent compounds such as diamond are assumed to be k = 0.83 and ν = 0.069, respectively, while in the strongly ionic KCl, Pugh's and Poisson's ratios are k = 1.19 and ν = 0.27, respectively 48,49 . Despite having different types of bonding, both KCl and diamond are brittle. On the other hand, strongly malleable gold exhibits ratios of k = 6.14 and ν = 0.42 47,50 , and gold is the most ductile metallic element. According to Christensen 51 , the critical value for the ductile-brittle transition appears at a Pugh's ratio of which corresponds to a critical Poisson's ratio of ν B/D = (3 √ 2 − 1) −1 ≈ 0.31. This criterion will be analyzed in more detail in a forthcoming work on cubic Heusler compounds 27 . Based on the calculated elastic constants, the studied tetragonal compounds exhibit Pugh's ratios between 1.38 and 2.2 and Poisson's ratios between 0.21 and 0.3, indicating that they have an intermediate behavior between ductile-and brittle-type behavior. The values of k and ν indicate the covalent or metallic character of these systems. More details about the bonding type may be found from a Bader analysis of the charge densities 52-54 .
The shear anisotropy factor A e provides a measure of the degree of anisotropy of the bonds between atoms in different planes. Tetragonal systems are described by two different shear anisotropic factors A 001 and A 100 (or equivalent A 010 ). Young's moduli of some of the studied Heusler compounds are shown in Figure 6. The shear anisotropy for {100} planes is considerably higher compared to {001} planes for Mn 2 NiGa and Rh 2 FeSn. This behavior is similar for Mn 2 FeGa and other Rh 2based compounds, as shown in Table IV. In contrast, in the case of Fe 2 MnGa and Mn 2 FeGa, the anisotropy for {001} is higher than that for {100}. Moreover, Rh 2 FeSn has an interesting distribution of Young's modulus: it is isotropic in the square x-y planes of the tetragonal structure, arising from the value of A 001 = 0.95, which is close to unity.

C. Virtual Crystal Approximation
One interesting feature of Heusler compounds, including Ni 2 MnGa, is their sensitivity to stoichiometric compositions. An infinitesimal deviation may lead to a phase transition or to changes in the electronic structure properties. Here, VCA was applied to explore offstoichiometric systems. This approximation is valid for  small changes of components with nearly the same radius, which holds for the considered systems. Moreover, VCA is only valid for neighboring elements and for small differences in the number of valence electrons (∆e − < 0.1).
In VCA, the atom Z A with charge Z is replaced by atom Z ′ A with the virtual charge Z ′ = Z ± ǫ to reflect that the average charge deviates from the original value at a certain position. In Ni 2 MnGa, if the site where 28 Ni resides is partially occupied by 25 Mn, the charge at that site will be lower; accordingly, the charge will be higher when 30 Ga occupies the same site. In parallel to the change in the nuclear charge, the number of electrons in the primitive cell changes to remain neutral. Thus, for 28±ǫ Ni 2 25 Mn 31 Ga, the number of valence electrons will be n v = 30 ± 2ǫ in the primitive cell, whereas the number of core plus semi-core electrons stays fixed at 72. The following cases are used in the present work: All elongated structures of off-stoichiometric Ni 2 MnGa have been fully optimized within the VCA approximation. As shown in Figure 5, a small change in the number of electrons at the Ni site has the most drastic effect on the energy landscape. The largest difference appears between Ga-rich ( 28.1 Ni) and Ni-poor ( 27.9 Ni) compounds.
Increasing Ni at Ga and Mn sites lowers the energy minimum at c/a > 1 compared to the stoichiometric compound. Conversely, increasing Ga at the Ni site increases the energy of c/a > 1 with respect to the stoichiometric compound. These results are in agreement with previously reported calculations 55 . Differences in the E(δ) dependence are more pronounced when the distorted structure is far from the initial structure. In the next step, the elastic constants of the off-stoichiometric Ni 2 MnGa were calculated using VCA for the tetragonal distorted structures with c/a > 1 at the lowest total energy (see Figure 5). Table IV summarizes the results of the elastic constant calculations for stoichiometric and off-stoichiometric Ni 2 MnGa. An extreme effect of the off-stoichiometry is reflected in the anisotropy ratio (A B ). As shown in Table IV, the 27.9 Ni and 30.9 Ga compounds exhibit negative anisotropies of -0.02 and -0.4, respectively. In fact, the negative anisotropy highlights the extreme instability of these systems. In particular, 30.9 Ga has a large negative value. The results explain why 30.9 Ga (Ga-poor) and 27.9 Ni (Mn-rich) Ni 2 MnGa synthesis is difficult. Moreover, the 28.1 Ni (Ga-rich) and 25.1 Mn (Mn-poor) compounds have large anisotropies, which are twice the value of the stoichiometric anisotropies. The large value also explains the tendency of phase transitions in these materials. Therefore, off-stoichiometric Ni 2 MnGa -nearly all synthesized samples are slightly off-stoichiometric -is expected to have pronounced phase transitions depending on the composition 56 . Thus, deficiency of valence electrons at the Mn or Ga sites in these systems leads to a negative anisotropy ratio and thus structural instability. Among the elastic constants, c 11 and c 22 are more strongly influenced when the composition changes compared to c 66 and c 44 , which remain nearly constant. Therefore, a small excess of each of the elements (small change of the valence electron concentration in the vicinity of that site) will not change the shear in the [100] direction. However, the asymmetries A B and A 100 significantly change compared to A 001 , reflecting a change of the in-plane chemical bonding. Here, Ni 2 MnGa is used as an example of the sensitivity of Heusler compounds on their stoichiometry. Disorder-induced phase transitions have been reported for other Heusler compounds such as iron-based compounds 57 .

D. Derived properties of tetragonal Heusler compounds
Finally, some physical properties and material parameters of the compounds are derived from the calculated elastic properties. The velocity of sound is an important quantity. Its averaged values can be directly determined from the calculated elastic constants. In experiments, on the other hand, the sound velocities can be used to measure elastic constants. Therefore, the sound velocities v are nearly synonymous with the elastic stiffness constants c. Further, sound velocities have been used to study various solid-state properties and processes 45 . Therefore, having the sound velocities predicted by calculations in advance could be quite important for experimental measurements. Usually, the directionally dependent acoustic properties are analyzed in terms of the slowness that is the inverse of the phase velocity. The group velocities are found from the derivatives of the slowness. Figure 7 compares the slowness surfaces of Fe 2 MnGa and Ni 2 MnGa. The slowness surfaces reflect the elastic anisotropy in comparison to Figure 6, which shows the distribution of Young's modulus. Three slowness surfaces appear in both cases, representing different polarizations of the sound wave. The pressure (p) wave is longitudinal polarized. Moreover, p has the highest phase velocity and thus the smallest slowness (Figure 7(a) and (d)). The remaining two surfaces belong to the fast (s 1 ) and slow (s 2 ) shear waves that are transversely polarized. The slowness surfaces of the p waves have a similar shape for both materials, and their maxima are found along the {001}-type principle axes. The shapes of the slowness surfaces of the shear waves differ between the two compounds. The observed differences reflect the differences in the anisotropy of both materials, and it is clearly seen that Ni 2 MnGa has a much lower anisotropy in the x − y plane.
Further material parameters are derived from the average sound velocities as described in Appendix C. At low temperatures, where only acoustic vibrational modes contribute, the Debye temperature θ D can be estimated from the average sound velocity 58 . The values estimated in this way are generally larger than Debye temperatures determined from phonon calculations or in experiments 59 because the optical phonon branches are neglected when the elastic constants are used for calculating the "acoustical" Debye temperature θ acc D . In a similar way, the average sound velocities can be used to estimate the "acoustical" Grüneisen parameter ζ ac 60 .
The calculated average sound velocities together with the estimated Debye temperatures and Grüneisen parameters are listed in Table V. The average sound velocities are rather similar for all compounds, ranging from about 3100 to 3800 m/s. The acoustical Debye temperatures are all above room temperature, ranging from 376 to 490 K. As expected, compounds consisting of heavier elements tend to have lower values. With the exception of Mn 2 NiGa, all ζ ac values are about 2. This demonstrates that the anharmonicity of the lattice vibrations is nearly the same for all compounds.
As further shown in Table V, the theory and experimental values have about a 20% discrepancy for the case of Ni 2 MnGa. As reported in Reference [61], however, the stoichiometry of the compound has a large effect on the measured Debye temperature. As shown in Table IV, the Debye temperature decreases by about 20% in 27.9 Ni. However, in the calculations, the changes in stoichiometry are extremely small. The changes in the calculated values may be more evident with larger variations in the stoichiometry. For example, in experiments, the results change from 261 K in the case of Ni 49.6 Mn 21.9 Ga 28.5 to 345 K in the case of Ni 53.1 Mn 26.6 Ga 20.3 . Therefore, having an ideal 2:1:1 system, like that assumed in most theories, is not easily possible or may even be impossible from an experimental standpoint. However, the estimation of measurable properties should provide information about the studied system and its potential for applications.

IV. SUMMARY
In the present work, the elastic constants of tetragonally distorted Heusler compounds were determined. The full-potential LAPW method together with the gradient-corrected PBE exchange-correlation functional were employed for all calculations. The relation between the calculated elastic constants and convergence criteria were discussed. Increasing only one of the parameters, such as the k-points or R MT K max , while keeping the other parameter low led to large errors in the calculated elastic constants. Therefore, to calculate both elastic constants accurately, R MT K max and k-points must be sufficiently large to guarantee convergance. Structural optimization was shown to have an important effect on the elastic constants for tetragonal Heusler compounds. The method was used to investigate the crystalline stability of materials based on the calculation of their elastic properties.
Based on the calculated results, the considered tetragonal Heusler compounds are intermediate materials, between brittle and ductile. Elastically, they exhibit mainly metallic rather than covalent bonding. The structural instability, mechanical properties, structural anisotropy, and other mechanical properties were also explored. Us- The calculated material properties are useful for applications focusing on bulk materials. However, the appearance and prediction of metastable tetragonal structures allows lattice parameter engineering with artificial c/a ratios initialized by epitaxial or pseudomorphic thin film growth. Thus, Heusler thin films could be designed to have specific properties. The equations that describe the elastic properties of solids have been described in detail by Nye 62 ; this discussion is summarized here and compares tetragonal, hexagonal, and cubic systems with focus on the tetragonal case. The strain matrix ǫ transforms lattice A with basis vectors X, Y, Z into the deformed lattice The elastic relations (Hooke's law) between the strain (ǫ) and stress (σ) matrices are mediated by the elastic compliance (S) or the elastic stiffness (C) matrices: From the elastic equations, the relations between the compliance matrix and the stiffness matrix are and vice versa C = S −1 . These relations imply that SC = CS = 1.
In the most general case, the elastic matrix is symmetric and of order 6 × 6. In triclinic lattices, the elastic matrix contains 21 independent elastic constants. This number is largely reduced in high symmetry lattices. For example, in an isotropic system, it contains only the two constants c 11 = c 22 = c 33 and c 12 = c 13 = c 23 , and the remaining diagonal elements of the matrix are determined by c 44 = c 55 = c 66 = (c 11 − c 12 )/2.
In cubic lattices, the three elastic constants c 11 , c 12 , and c 44 are independent. There are five independent elastic constants for hexagonal structures (c 11 , c 12 , c 13 , c 33 , and c 44 ), while tetragonal structures have either seven (classes: 4, 4, or 4/m) or six (c 11 , c 12 , c 13 , (c 16 = −c 15 ), c 33 , c 44 and c 66 ) elastic constants. The elastic matrix for all classes of cubic and hexagonal crystals as well as the classes 4mm, 42m, 422, or 4/mmm of tetragonal crystals have the form where zero elements are assigned by dots and additional tetragonal elements for classes 4, 4, or 4/m are given in brackets. Moreover, the elastic matrix has restrictions c 33 = c 11 , c 66 = c 44 , c 13 = c 12 in cubic systems and c 66 = (c 11 − c 12 )/2 in hexagonal systems.
The relations between the elastic constants c ij and the elements of the compliance matrix s ij are found from Equation A4. In all classes of hexagonal systems or in tetragonal systems belonging to the classes 4mm, 42m, 422, or 4/mmm, the relations between c ij and s ij are given by where c 66 appears only in tetragonal systems. Indeed, the number of equations is much less in cubic systems as shown from the restrictions given above.
The elastic properties of single crystals are completely determined by the elastic matrices C and S. In reality, polycrystalline materials are considered more often than single crystals. Polycrystalline materials consist of randomly oriented crystals, and thus, a description of their elastic properties requires only two independent elastic moduli: the bulk modulus (B) and the shear modulus (G). The relationships between the single-crystal elastic constants and the polycrystalline elastic moduli are given by the Voigt 63 or Reuß 64 averages. Voigt's approach uses the elastic stiffnesses c ij , while Reuß's approach uses the compliances s ij . Voigt's moduli 63 are given as function of the elastic constants by the equations:  66 .
For cubic or isotropic crystals, the bulk moduli in Voigt's (B V ) and Reuß's (B R ) approach are equal, as shown by using the restrictions on c ij given above. In cases other than isotropic or cubic, G R cannot be easily rewritten in terms of the elastic constants.
Finally, the mechanical properties of polycrystalline materials are approximated in the Voigt-Reuß-Hill 65 approach, where the bulk and shear moduli are given by arithmetic averages: The bulk modulus B of a material characterizes its resistance to fracture, whereas the shear modulus G characterizes its resistance to plastic deformations. Therefore, ratios between the elastic moduli B and G are often given for characterization and comparison of different materials. Pugh's modulus k is the simple ratio of the bulk and shear moduli 47 : Poisson's ratio ν also relates the bulk and shear moduli: Further, Poisson's ratio bridges between the rigidity modulus G and Young's modulus E, which is given by Appendix B: Elastic stability and representation of elastic properties The set of elastic moduli and their ratios allows characterization of the elastic behavior of materials. However, the mechanical stability is still open. As a first criterion, the elastic moduli all must be positive. Born and coworkers developed a theory on the stability of crystal lattices 66,67 . For tetragonal crystals at ambient conditions, the seven elastic stability criteria are given by Note that the number of criteria is reduced for the lower number of elastic constants in hexagonal or cubic crystals to 5 or 3, respectively. The last condition is used to define the tetragonal shear modulus c ′ = (c 11 − c 12 )/2. In some works, the direct difference C ′ = c 11 − c 12 is used. If external, hydrostatic pressure p is applied, then the crystal becomes unstable when 2p > C ′ , that is, at p > c ′ .
The linear compressibility β is the crystal response to hydrostatic pressure by a length decrease. For cubic systems, the linear compressibility is isotropic, that is, a sphere of a cubic crystal under hydrostatic pressure remains a sphere. The situation is different in non-cubic systems where β = β(r) becomes directionally dependent. In hexagonal, trigonal, and tetragonal systems, the directional dependence is given by β(r) = (s 11 +s 12 +s 13 )(x 2 +ŷ 2 )−(s 11 +s 12 −s 13 −s 33 )ẑ 2 .
For the tetragonal classes 4, 4/m, and 4, an additional term is present such that The shear anisotropic factors provide a measure of the degree of anisotropy in the bonding between atoms in different planes. The number of different shear anisotropies depends on the crystal system. In both -hexagonal and tetragonal -systems, the shear anisotropic factors A 100 (or equivalent A 010 ) for the {100} shear planes between the 011 and 010 directions and A 001 for the {001} planes between 110 and 010 are: In cubic crystals, both factors are the same A e = A 001 = 2c 44 /(c 11 − c 12 ), as mentioned above. In hexagonal systems, c 66 = (c 11 − c 12 )/2, and thus, A 001 = 1. For isotropic crystals, all A factors must be unity, while any value smaller or greater than unity is a measure of the degree of elastic anisotropy possessed by the crystal.
Comparing the equations (B3,B6) for the elastic anisotropies with the Born-Huang 67 criteria, these equations can clearly be used to show the elastic stability. Most obviously, crystals with one negative anisotropy are not stable. Further, crystals with large anisotropies also tend to instabilities; in particular, crystals are not stable for A → ∞ when one of the denominators becomes zero. This behavior makes the anisotropies important parameters, even though they may not cover all possible causes for Born-Huang instabilities.
In general, it is calculated from logarithmic derivatives of the vibrational frequencies with respect to the crystal volume. However, full phonon calculations as function of crystal volumes are demanding tasks, and fast estimates are thus welcome. Belomestnykh 60 derived an "acoustical" Grüneisen parameter ζ ac that is directly related to the sound velocities. Therefore, ζ ac is given by . (C8)