Temperature-dependent elastic properties of binary and multicomponent high-entropy refractory carbides

Available information concerning the elastic moduli of refractory carbides at temperatures (T) of relevance for practical applications is sparse and/or inconsistent. We carry out ab initio molecular dynamics (AIMD) simulations at T = 300, 600, 900, and 1200 K to determine the temperature-dependences of the elastic constants of rocksalt-structure (B1) TiC, ZrC, HfC, VC, and TaC compounds as well as multicomponent high-entropy carbides (Ti,Zr,Hf,Ta,W)C and (V,Nb,Ta,Mo,W)C. The second order elastic constants are calculated by least-square fitting of the analytical expressions of stress vs. strain relationships to simulation results obtained from three tensile and three shear deformation modes. Moreover, we employ sound velocity measurements to evaluate the bulk, shear, elastic moduli and Poisson's ratios of single-phase B1 (Ti,Zr,Hf,Ta,W)C and (V,Nb,Ta,Mo,W)C at ambient conditions. Our experimental results are in excellent agreement with the values obtained by AIMD simulations. In comparison with the predictions of previous ab initio calculations - where the extrapolation of finite-temperature elastic properties accounted for thermal expansion while neglecting intrinsic vibrational effects - AIMD simulations produce a softening of elastic moduli with T in closer agreement with experiments. Results of our simulations show that TaC is the system which exhibits the highest elastic resistances to both tensile and shear deformation up to 1200 K, and identify the high-entropy (V,Nb,Ta,Mo,W)C system as candidate for applications that require good ductility and toughness at room as well as elevated temperatures.


Introduction
Elastic moduli and elastic stiffness constants of crystals have primary importance in materials design and materials discovery. For example, the search of novel solid crystal phases with potentially new properties includes verifying the Born von Kármán conditions of elastic stiffness constants for mechanical stability [1]. The temperature variation of the elastic anisotropy has a fundamental impact on the material's microstructure evolution and thus on observed age hardening mechanisms in ceramic coatings [2,3]. Hardness and toughness of refractory ceramics are discussed in terms of the residual stress in the material, which is evaluated through the knowledge of elastic stiffness constants.
Indeed, ab initio investigations have shown that the predicted elastic constants may vary up to ≈20% depending on deformation path, strain range, and order of the polynomial used to fit energy vs. strain curves [4]. In regard to the identification of materials with superior mechanical or structural performance, the observed trends in hardness [5,6], strength [7,8], toughness [9-11] and ductility [7,8,12,13] are often correlated to trends in elastic constants of the materials.
The elastic constants of solid crystal phases can be efficiently screened by high-throughput ab initio calculations at 0 kelvin [14]. In contrast, realistic modeling and evaluation of mechanical properties as hardness and toughness, which would require simulation boxes with volumes of the order ≈10 3 -10 5 nm 3 [15], is currently unfeasible for first-principles methods. For this reason, ab initio calculated polycrystalline elastic moduli, such as shear modulus and Poisson ratio, are widely employed to design empirical indicators of hardness and toughness of compounds and alloys [9, [16][17][18][19][20] which, combined with the progress in machine-learning approaches [20][21][22], can lead to an increasingly more rapid identification of solids with enhanced mechanical performance.
In the case of refractory carbide, nitride, and boride ceramics, properties as hardness and toughness are of paramount importance, especially for high-temperature applications. Despite this, the elastic constants of these refractory ceramics have been primarily evaluated at 0 K. Ab initio estimations of elastic constants at finite temperatures are often done using static calculations at cell volumes determined in the quasiharmonic approximation for the thermal expansion (see, e.g., [23-

AIMD simulations
In this work, we employ tensile and shear stress/strain data within the elastic (prior to yielding) material's responses obtained from AIMD simulations performed at T = 300, 600, 900, and 1200 K [40] and evaluate the SOEC of B1 TMC as a function of temperature (T). The entire (elastic and plastic) stress vs. strain curves calculated for TMC subject to tensile and shear deformation will be reported elsewhere [35]. The procedure for AIMD stress/strain calculations is detailed in Refs. [32][33][34]. The B1 TMC supercells used in our simulations have three different crystallographic orientations (Fig. 1), which are convenient for modeling elongation orthogonal to (001), (110), and (111) (i.e., surface terminations of lowest formation energies [41]) and shearing on the primary slip systems of B1-structure ceramics, namely, {001}〈11 ̅ 0〉, {110}〈11 ̅ 0〉, and {111}〈11 ̅ 0〉 [42][43][44][45]. The stress vs. strain relationships are fitted with 2 nd -order polynomials. Then, the second-order elastic response of the crystal is calculated from the slope of the polynomial in the limit of vanishingly small strain. The approach is analogous to fitting the energy-density vs. strain curves with 3 rd -order polynomials, as illustrated, for example, in Ref. [4].
AIMD simulations are carried out with VASP code [46], which implements the projector augmented-wave method [47]. The Perdew-Burke-Ernzerhof (PBE) electronic exchange and correlation functionals of the generalized gradient approximation (GGA) [48], Γ-point sampling of the Brillouin zone, and planewave cutoff energy of 300 eV are used in all simulations. The ionic equations of motion are integrated at timesteps of 1 fs, using 10 -5 eV/supercell convergence criteria for the total energy during both system equilibration and mechanical testing. Prior to modeling tensile and shear deformation, the supercell equilibrium structural parameters are optimized at each investigated temperature via an iterative NPT [49] simulation procedure, which employs the Parrinello-Rahman barostat [50] and the Langevin thermostat [51]. Subsequently, NVT [52] sampling of the configurational space is used to equilibrate the unstrained structures for 3-5 additional ps using the Nosé-Hoover thermostat [51]. In these simulations, we ensure that the time-averaged |xx|, |yy|, and |zz| stress components are ≲ 0.3 GPa. The structural parameters determined for unstrained (equilibrium) supercells at different temperatures allow us to evaluate lattice constants and average linear thermal expansion coefficients of B1 TMC.

Material synthesis and elastic constants measurements
For fabricating the high entropy carbides used in this study, binary carbide powders were weighted out in 24 gram quantities to achieve an equiatomic mixture of the 5 metallic species, additional graphite was added to yield the stoichiometric monocarbide ratios when Mo2C was included. The binary carbides included: TiC, ZrC, HfC, VC, TaC, Mo2C and WC, all obtained from Alfa Aesar, USA. The powders were blended and hand mixed to achieve the following bulk compositions: (Ti,Zr,Hf,Ta,W)C and (V,Nb,Ta,Mo,W)C, followed by high energy ball milling in 30 minute increments, for a total of 120 minutes, using 10-minute cooling periods between each step.
The powder samples were then compacted and sintered using current and pressure assisted densification (CAPAD) in a well-controlled vacuum environment. The samples were sintered using 20mm diameter graphite dies and plungers. The initial vacuum condition was 20 mTorr, and upon heating, whenever the vacuum decreased to 100 mTorr, the heating profile was paused to allow the vacuum to return to the 20 mTorr level. The heating and pressure profile used in the sintering of both materials was as follows: using an initial pressure of 5 MPa, the samples were heated at 100C per minute to 2200C, then held at temperature for 9 minutes, before increasing the pressure to 80 MPa over one minute, then holding at this pressure and 2200C for 15 minutes. The samples were subsequently cooled to room temperature in vacuum under 5 MPa pressure. Each sample was ground on both circular faces to remove any excess graphite from the die setup, then polished using standard metallurgical techniques to a final polish of 0.04 µm colloidal silica. The carbon-to-metal ratios assessed for our HEC specimens are 0.90±0.03 [55].
To obtain the elastic moduli, both longitudinal and shear, the method outlines in ASTM E494-15 (Standard Practice for Measuring Ultrasonic Velocity in Materials) were followed, using the 20 mm diameter and ~3 mm thick sintered high entropy carbide samples. The samples were placed on an acoustically-isolated substrate for wave speed measurements using separate piezoelectric transducers for longitudinal and shear velocities. A total of 4 separate measurements, in each velocity direction, were obtained in order to compute the respective wave speeds from which to compute the different elastic properties. To compute the various elastic constants, both the respective wave speeds and the sample density must be measured, with the latter measured by Archimedes method. The following equations were used to compute the elastic properties, which assume isotropic elastic responses. For Poisson's ratio: v = (VL 2 -2VS 2 )/[2(VL 2 -VS 2 )]; for Young's modulus: E = 2 VS 2  (1+ v); for shear modulus: G = VS 2 ; and for bulk modulus: B = E/[3(1-2v)], where  is density, VL is the measured longitudinal wave speed and VS is the measured shear wave speed.

Results and discussion
The calculated and experimental lattice parameters and average linear thermal expansion coefficients  of B1 TiC, ZrC, HfC, VC, TaC, (Ti,Zr,Hf,Ta,W)C, and (V,Nb,Ta,Mo,W)C at temperatures ranging from 0 to 1200 K are presented in Table 2 (AIMD data listed in bold).
Previous 0-K ab initio calculations based on GGA slightly overestimate (by less than 1%) the experimental lattice parameters of B1 TMC (see Table 2). The discrepancy is due to known limitations of GGA functionals [56], as well as to the fact that the synthesized transition-metal carbides typically contain anion vacancies, which reduce equilibrium volumes [57,58]. At deviance with ab initio results at 0 K, our AIMD simulations yield lattice constants that are systematically ≈0.5% smaller than experimental values measured at 300 K. This is likely due to the relatively lower accuracy of present simulations (details in Section 2.1) in comparison with static ab initio calculations, for which the structural parameters are converged with respect to cutoff energies and kpoint mesh thicknesses at low computational efforts. Overall, our AIMD results closely match with the lattice constants and thermal expansion coefficients determined by experiments at different temperatures ( (1/E300K)×(dE/dT) ≈ -1.1±0.2  10 -4 K -1 for temperatures increasing from 300 to 1200 K. Note that the decay coefficient varies with the sample stoichiometry [67]. A decrease in elastic stiffness vs. T is expected for crystal phases that are dynamically stable from 0 K to melting points. The monotonic reduction in SOEC with increasing temperature [68][69][70] is also reflected by a corresponding decrease in mechanical strength and hardness [42,71]. Deviations from the typical behavior of SOEC vs. Tsee, e.g., temperature-induced increase in the C44 of some metals [72,73] may arise from strainor thermally-induced modifications of electronic-structures.
All TMC modelled in this work display expected trends in C11 and C44 vs. temperature. Both elastic constants decrease monotonically from 300 to 1200 K ( Fig. 2a and 2c). Conversely, the values of the C12 elastic constants appear to be less affected by temperature. Indeed, the C12 calculated for TaC, HfC, and (Ti,Zr,Hf,Ta,W)C remains nearly unvaried (within statistical uncertainty) with T (Fig.   2b). The trends that we obtain for all other carbides indicate slow C12 reductions with increasing temperature (Fig. 2b). The diminution in C11 and C44 implies that the moduli B, G, and E decrease as a function of temperature for all materials considered here ( Fig. 2e-f).
Consistent with the temperature dependence observed for the polycrystalline moduli G and E, the directional elastic E001, E110, E111 and shear G001, G110, G111u moduli of the investigated carbide systems decrease with increasing T (Figs. 2c, 3a- [42,43,74]. In contrast, {001}〈11 ̅ 0〉 glide in carbides has been experimentally observed in very few cases [44,45]. Nonetheless, a stiffer elastic response to change of shape (quantified by the shear moduli) does not necessarily imply that the shear strength (resistance to slip) is larger. For example, B1 VN0.8 exhibits greater elastic shear stiffness, but much lower shear strength than B1 VN [32]. Analogously, although TaC displays the highest shear moduli G111u among all TMC (Fig. 3b), {111}〈11 ̅ 0〉 plastic deformation of B1 TaC crystals can be observed already at room temperature [74].
At relatively low temperatures, mechanical strength and plastic deformation in B1-structure ceramics is primarily controlled by dislocation mobilities [43][44][45]74]. The mobility is, in turn, dictated by the core structure of the line defect. An atomic-level determination of dislocation core structures in ceramics can be a formidable task. The problem is further complicated by the fact that dislocation cores may exhibit different polymorph variants [75,76] and be affected by lattice stoichiometry [77].  . 3f). Conversely, B1 VC, TaC, and (V,Nb,Ta,Mo,W)C are predicted to become more ductile (and less strong) with increasing temperature (see arrow in Fig. 3f).
Previous ab initio studies conducted on refractory B1-structure ceramics indicated that an The theoretical prediction that TaC and VC should be more ductile than TiC, ZrC, and HfC is consistent with the fact that {111}〈11 ̅ 0〉 slip in Group-VB carbides is activated at lower temperatures than in Group-IVB carbides [42]. The operativity of lattice slip on {111} crystallographic planes is of crucial importance for ductility in B1-structure ceramics, because it provides a sufficient number of independent pathways to induce plastic deformation [78]. In agreement with the conclusions of Ref. [42], other experiments have shown that the concentration of intrinsic {111} stacking faults in B1 TaC is much greater than in, e.g., B1 HfC [79]. The observation has been rationalized by ab initio results, which indicated that these extended defects are metastable in Group-VB, whereas they are unstable in Group-IVB carbides (see [79] with supplemental material). The energetically-favored formation of {111} stacking faults facilitates, in turn, plastic deformation via {111}〈112 ̅ 〉 synchro-shear mechanisms [79].
The results of our AIMD simulations demonstrate that TaC is the system of highest elastic resistances to tensile (C11, E, E001, E110, E111) as well as shear (C44, G, G110, G111u) deformation at all investigated temperatures, Figs. 2 and 3. VC is the crystal with second highest tensile and shear elastic stiffnesses (with a few exceptions for T ≈1200 K). Conversely, we find that B1 ZrC generally displays the lowest values of G and E elastic stiffnesses. With regard to high entropy alloys, Fig. 2a shows that the C11 elastic constants of (Ti,Zr,Hf,Ta,W)C and (V,Nb,Ta,Mo,W)C are nearly identical, and with values intermediate to those of Group-IV and Group-V carbides at the temperatures probed in the present work. Our AIMD results also demonstrate that (V,Nb,Ta,Mo,W)C possesses the lowest C44 and highest C12 values for 300 ≤ T ≤ 1200 K (Fig. 2c). As discussed above and illustrated in for B1-structure carbides and nitrides of similar VEC (≈9.5 e -/f.u). [84,85]. Hindered dislocation motion across the faults provides high hardness and mechanical strength [84,85]. (V,Nb,Ta,Mo,W)C does indeed exhibit rather high hardness values (27±3 GPa [37]). Concurrently, a substantial amount of Mo and W in the host B1 (V,Nb,Ta,Mo,W)C lattice should promote slip on {111} lattice planes [86], which fulfills the criteria for plasticity in B1-structure crystals [78]. The elastic moduli of the two HEC systems investigated in this work have been previously assessed by experiments and ab initio calculations [36]. Our AIMD results for B, G, E, and ν of the HECs at 300 K closely match (deviations < 5%) with first-principles elastic properties at 0 K [36] and are also in good agreement with previous experimental data [36] (see Tables 8 and 9 Before conclusions, we offer a comparison of different computational techniques for the evaluation of finite-temperature elastic properties. As anticipated at the beginning of this section, finite-temperature elastic properties of some refractory nitrides and carbides have been calculated via the QHA [23,27,62] or volume-scaling methods [28,31] in the framework of static first-principles techniques. These, however, generally produce a slower decay in elastic stiffnesses vs. T than methods which account for intrinsic anharmonic effects, as demonstrated for B1 Ti1-xAlxN [28] and CrN [31] refractory compounds. Here we extend the analyses of Refs. [28,31] by showing trends in elastic constants vs. T obtained by AIMD simulations, experiments, and static ab initio calculations. ZrC, and (c) TiN. Note that B1 TiN is included in our analysis because of the sparseness of temperature-dependent elastic properties for refractory carbides. A detailed description of the elastic properties of TiN and other nitrides will be provided elsewhere [35]. In all three cases presented in  Table 3, and B and E vs. T for ZrC, Table 4). This suggests that the inclusion of vibrational effects is of particular importance to achieve accurate temperature-dependent elastic shear stiffnesses C44, which, in turn, strongly affect the values of polycrystalline shear moduli G.

Conclusions
We  . Note that the actual supercells used in our AIMD simulations of tensile and shear deformation contain 576 atoms each (see figure 1 in [33]). Spheres of different colors indicate the metal and carbon sites of the B1 crystal structure. This figure is generated using the VMD software [89].  The solid lines, which serve as guide to the eye, are 2 nd order polynomials obtained by imposing zero slope at 0 K and least-squares differences with the elastic properties determined by AIMD. The statistical uncertainty on the calculated elastic properties is of the order of the difference in pressure (GPa) between a symbol and the corresponding fitted curve. In (f), the empirical criteria suggest that temperature has the general effect of improving ductility while reducing strength. However, some carbides exhibit a non-monotonic behavior in such sense (see Tables 3-9).  [30]. Ab initio results that account for both thermal expansion and anharmonic effects (solid blue line) are from Ref. [28]. Results of static ab initio calculations: purple circles (QHA) [62] and green dotted line (volume scaling) [28]. The experimental results (red diamond) are from Ref. [95]. The figure shows that the temperature-induced decay in experimental shear elastic stiffness is better reproduced by AIMD simulations than by static ab initio calculations that account for thermal expansion but neglect lattice vibrations. (C11-C12+C44)/3 (C11-C12-2C44)/(3√2) 0 Table 1. Correspondences between stress variations [∂ij/∂= obtained during tensile or shear deformation for vanishingly small strains and secondorder elastic constants of a cubic crystal. All ij values equal 0 mean that the derivative of that ij vs. strain component is equal zero for =0. We note that the supercell volume and shape remain constant at each strain step of AIMD/NVT simulations. Thus, the materials' elastic resistances to 〈001〉, 〈110〉, and 〈111〉 elongation do not correspond to the directional Young's moduli E001, E110, E111. A direct evaluation of Ehkl during tensile strain would require lateral supercell relaxation.