Ideal simple shear strengths of two HfNbTaTi-based quinary refractory multi-principal element alloys

Atomistic simulations are employed to investigate chemical short-range ordering in two body-centered cubic refractory multi-principal element alloys, HfMoNbTaTi and HfNbTaTiZr, and its influence on their ideal simple shear strengths. Both the alias and affine shear strengths are analyzed on the {110} and {112} planes in the two opposing ⟨ 111 ⟩ directions. In both quinary alloys, local ordering of NbNb, TaTa, HfNb, HfTa


I. INTRODUCTION
Refractory multi-principal element alloys (RMPEAs) have been identified as promising candidate materials for aerospace propulsion systems, land-based gas turbines, nuclear reactors, heat exchanger tubing, and the chemical process industry. 1 These materials consist of at least three metallic elements with equal-or near-equal molar ratio with most or all constituent elements being refractory metals. 2 Two equal-molar quinary RMPEAs, HfMoNbTaTi and HfNbTa-TiZr, possess outstandingly high room temperature (RT) strengths. Take the compressive yield strength as an example. It is 1369-1713 MPa for of HfMoNbTaTi 3,4 and 900-1073 MPa for HfNbTaTiZr. [5][6][7][8][9] Unlike pure refractory metals, many RMPEAs can retain ultrahigh strengths over a broad range of elevated temperatures. For example, the compressive yield strength of HfMoNbTaTi decreases from 1369 MPa to 699 MPa (i.e., −49%) when the temperature increases from RT to 1473 K. 3 The percentage decrease is smaller than those in pure refractory metals. For example, in Nb and Ta, the ultimate tensile strengths at RT, 290 and 441 MPa, respectively, decrease to 65 and 106 MPa at 1473 K. 10 The relative change is −78% in Nb and −76% in Ta.

A. Interatomic potentials
Embedded-atom method (EAM) potentials 42 are used to describe the interatomic interactions. The six elemental potentials are Hf, 43 Mo, 44 Nb, 45 Ta, 44 Ti, 44 and Zr 44 and the cross interactions among these elements are based on the formulations by Johnson 46 and Zhou et al. 47 The resulting interatomic potentials built for the HfMoNbTaTi and HfNbTaTiZr RMPEAs are referred to hereinafter as the alloy potentials. Based on them, we construct two A-atom potentials, which provide mean-field representations of the RMPEAs, 48 resulting in two artificial pure metals, denoted as HfMoNbTaTiA and HfNbTaTiZrA, respectively. Except the HfNbTaTiZr alloy potential, 43 all four potentials are newly developed here. They can be downloaded at https://github.com/shuozhixu/APLmater_2022.
To assess the accuracy of these potentials, we performed density functional theory (DFT) calculations for the lattice parameters (a 0 ) and effective BCC elastic constants (C † 11 , C † 12 , C † 44 ) for the random structures of the two RMPEAs and compared these quantities with those from the alloy potentials. The atomistic and DFT calculations here are performed via Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 49 and Vienna Ab initio Simulation Package (VASP), 50 respectively, following our previous work. 41 We also compare the alloy potentialbased results with available DFT or experimental data in the literature. In Table I, it is found that the alloy potentials generally well predict the structural parameters compared with the present and prior DFT results. In particular, our HfNbTaTiZr alloy potential, which yields C † 12 = 111.83 GPa, performs better than a modified EAM potential, 34

B. Ideal shear strength
Generally, the ideal shear strength is the minimum stress needed to plastically deform an infinite defect-free material and is an upper bound to the shear strength of the material. 51,52 It is not related to the shear strength of a material that already contains defects. There are two types of ideal shear: the alias shear and the affine shear. 53 The alias shear strength is related to the potential energy to displace two crystalline halves across a given slip plane. The affine shear strength allows the entire crystal to deform to accommodate the homogeneously applied shear. In each type, the applied shear can be either a pure shear or simple shear. These TABLE I. Lattice parameters a 0 (in Å) and effective BCC elastic constants C † 11 , C † 12 , C † 44 (in GPa) of different models in two RMPEAs. Superscripts § and ‡ denote data from experiments and DFT calculations, respectively. Models without superscripts are from atomistic simulations of size D structures in Table II, among which one set of data is based on a modified EAM potential of HfNbTaTiZr. 34 All unreferred data from DFT or alloy potentials are newly calculated here, where a 0 and C † ij are obtained using the energy-volume method 35 and stress-strain method, 36 respectively. Data of three constituent BCC pure metals, based on the elemental interatomic potentials, are also presented as references. Experimental and computational elastic constants were based on polycrystals and single-crystals, respectively.  56 revealing that (i) the simple shear strength is higher than the pure shear strength and (ii) the affine shear strength is higher than the alias one. The first statement was found to hold in pure BCC metals. 54,57 However, the affine and alias shear strengths have not been compared in any pure BCC metals or RMPEAs. Here, we consider both alias and affine shear strengths in simple shear. Compared with the pure shear strength, the simple shear strength is more closely related to the slip resistance for a dislocation when the shear plane is selected as the glide plane. 58 Between the two RMPEAs studied here, dislocations were experimentally investigated only in HfNbTaTiZr, 7,9,59 suggesting that screw dislocations on {110} and {112} planes are prevalent. Our recent atomistic modeling work in six RMPEAs also showed that these two planes are operative during plastic deformation. 36,60,61 Therefore, two shear planes, {110} and {112}, are considered in this work.

Simulation cells for the RMPEAs
We first construct special quasi-random structures (SQSs) in which the five elements are randomly distributed 62 using ATAT. 63 Two size C and one size D cells are employed to generate structures with CSRO. We first determine the chemical potential differences between Nb and the other five elements under the semi-grand-canonical (SGC) ensemble at 1500 K via the hybrid molecular dynamics (MD)/Monte Carlo (MC) technique using a small cubic cell containing 2000 atoms. The values obtained are summarized in Table S1 in the supplementary material. Next, based on cells of size C or size D, we perform hybrid MD/MC simulations under the variance-constrained SGC (VCSGC) ensemble at three annealing temperatures of 300, 600, and 900 K, respectively. Each MC cycle is accompanied by 20 MD steps. After 600 000 MC cycles, the potential energy converges to a constant, and so the equilibrium configurations with the equal-molar composition and thermodynamically correct occupation of atomic sites are achieved. The last step is to quench these structures to 1 K and to minimize the system energy with all three normal stresses being zero. More details of the hybrid MD/MC simulations can be found in our previous work. 65 As a result, for each RMPEA, we obtain nine atomic structures with CSRO, in addition to the eight random ones in Table II. Eight size C atomistic structures in the two RMPEAs are presented in Fig. 1. Note that the parallel algorithm under the VCSGC ensemble requires that each edge length of the simulation cell must be no less than four times the cutoff radius of the EAM potential, 66 28.62 Å, thus excluding size A and size B cells.
CSRO is measured using the nearest neighbor shell-based Warren-Cowley (WC) parameters: where p ij denotes the probability of a j-type atom being around an atom of type i within the shell, cj is the concentration of j-type atom, and δij is the Kronecker delta function. 67,68 Note that αij = αji in equal-molar MPEAs. In a completely random structure, αij = 0 for all ij pairs. A tendency for segregation or local ordering corresponds to a positive αij for pairs of the same species (i.e., i = j) or a negative αij for those of different species (i.e., i ≠ j). In either case, a higher αij magnitude suggests a higher degree of CSRO.
In what follows, we use "number+K+MDMC" to denote the structure obtained by hybrid MD/MC simulations at a specific annealing temperature. For example, "600KMDMC" represents the structure created by hybrid MD/MC simulations at 600 K. The values of structural parameters of size D structures with CSRO are provided in Table I. In general, the presence of CSRO influences C † 11 , C † 12 , and C † 44 much more than it does a 0 .

Alias simple shear strength
The alias simple shear strength, T al , is the maximum stress for the sliding between two adjacent atomic planes along a certain direction. In BCC metals, T al has been related to the generalized stacking fault energy (GSFE) curve γ gsf . 69,70 All eight random structures and three structures with CSRO are used to calculate γ gsf . A 20 Å thick vacuum is added, along the x and y directions, respectively, to the {110} and {112} cells (Table II). Along all directions, periodic boundary conditions (PBCs) are applied. The top atomic planes are displaced by d with respect to the bottom ones along the z⟨111⟩ direction. Following each displacement, energy minimization is performed, in which the top two and bottom two atomic planes are fixed while the remaining ones are allowed to move along the x and y direction, respectively, in {110} and {112} cell. A GSFE curve is then obtained as a function of d from 0 to √ 3a 0 /2. More details of the GSFE calculations in RMPEAs can be found in Refs. 36 and 41, which were focused on random solid solutions. The procedure is repeated for several parallel planes with distinct atomic arrangements. We obtain 80 γ 110 gsf curves and 100 γ 112 gsf curves for each cell in each RMPEA. Each GSFE curve exhibits a single peak value, called the unstable SFE γ usf , i.e., In pure BCC metals, the {112} plane exhibits a geometric asymmetry in shear resistance in the twinning (soft) and anti-twinning (hard) ⟨111⟩ directions. In RMPEAs, the fluctuations in composition across the glide plane can induce a chemical asymmetry in the resistance to shear between positive and negative ⟨111⟩ directions, even for the {110} plane. 41 Here, we use T al-s to denote the soft sense and T al-h the hard sense, as defined by

Affine simple shear strength
To calculate the affine simple shear strength T af , four size D cells-one SQS and three structures with distinct CSRO-are adopted for each RMPEA. Size D, instead of size C, is chosen because it is sufficiently large to accommodate possible defect nucleation and evolution. PBCs are applied to all directions. An affine simple shear is applied to each cell, with the stress-strain response recorded. The shear plane is either {110} or {112}, while the shear direction is one of the two opposite ⟨111⟩ directions. T af corresponds to the peak stress in the stress-strain curve, when defects (e.g., dislocations and twins) or phase transformation (PT) start to emerge within the crystal. Similar to T al , T af possesses a soft and hard sense, i.e., T af-s and T af-h , for each plane in each structure. To assess temperature and strain rate effects, calculations for a few selected cases are repeated over a series of strain rates from 5 × 10 7 s −1 to 5 × 10 11 s −1 and temperatures from 0 to 300 K. Figures S1 and S2 of the supplementary material show that T af converges for strain rates ≤5 × 10 8 s −1 and decreases linearly with increasing temperature. Unless stated otherwise, most simulations in the current work are carried out at 5 × 10 8 s −1 and 1 K, at which the CSRO effects are important. A total of four T af are obtained for each sense on each plane in each RMPEA. For reference, T al and T af of five pure metals-two A-atom ones, HfMoNbTaTiA, HfNbTaTiZrA, and three BCC constituents, Mo, Nb, and Ta-are also calculated using similar procedures. Results of the last three metals are summarized in Table S2 in the supplementary material.

A. Size effects in the alias shear strength
Calculating γ usf and T al requires selecting a cross-sectional area A. Figure 2 shows the effect of A on the average (T al ) and standard deviation in T al (S al ) for two RMPEAs. For the broad range of sizes considered, A is found to have a small effect onT al , with a <3% variation for HfMoNbTaTi and <6% for HfNbTaTiZr. The standard deviation S al , on the other hand, is affected by A and decreases with increasing A. These size effects are similar to those found in a study on the intrinsic SFE in FCC equal-molar FeNi binary and unequalmolar CrFeNi MPEA. 71,72 Here, the size effect is relevant and strong when A ≤ 1 nm 2 , about the size of a dislocation core. In this regime, HfMoNbTaTi has a larger standard deviation, but similar coefficient of variation, with respect to HfNbTaTiZr.  The size effect is also examined for the quinaries with CSRO by repeating the γ usf and T al calculations for structures of sizes C and D. Since the chemical ordering is short-range, the influence of changing CSRO is expected to diminish for the larger A. We find that for the random and 900KMDMC cases, the values are the same between the two sizes. For 600KMDMC and 300KMDMC, the effects of CSRO are slightly subdued with the larger size D. Apart from this, CSRO does not change the fact that the quinaries have lowerγ usf andT al than the three BCC constituents, that HfMoNbTaTi has higherγ usf andT al values than HfNbTaTiZr, and that the standard deviation for size D is lower than that for size C. The last point is demonstrated in Fig. S4 in the supplementary material.

B. CSRO
The WC parameters for all possible pairs in both quinaries are provided in Tables S3 and S4 in the supplementary material. At the lowest annealing temperature of 300 K, results for HfNbTaTiZr show that NbNb, TaTa, HfNb, HfTa, HfTi, NbTa, NbZr, TaZr, and TiZr pairs are likely, while in HfMoNbTaTi, the NbNb, TaTa, TiTi, HfMo, HfNb, HfTa, MoNb, MoTa, and NbTa pairs would tend to segregate. In general, the median of the absolute WC parameter values is higher for HfNbTaTiZr than for HfMoNbTaTi. In the 300KMDMC case, for instance, the median is 0.224 for HfMoNbTaTi compared to 0.242 for HfNbTaTiZr. Figure 3 shows the evolution of the five most significantly favored elemental pairs in the two RMPEAs, as the annealing temperature decreases. The most favored pair is TiTi and HfTi, respectively, in HfMoNbTaTi and HfNbTaTiZr.
Between the two RMPEAs, the CSRO in only HfNbTaTiZr has been previously investigated. Experiments have found segregation of the NbTa pair between 773 and 1073 K 9,13,73-76 and the HfTa pair between 873 and 1073 K. 77 At a high temperature of 2773 K, ab initio MD simulations found that CSRO is absent. 78 Based on a modified EAM potential, 34 Huang et al. 79 showed that four unlike atomic pairs tend to segregate at 800 K: HfTi, HfZr, NbTa, and TiZr. First principles calculations at 1002 and 1298 K 80 and CALPHAD modeling between 273 and 1073 K 75 identified segregation of HfZr and NbTa pairs. First principles calculations also predicted TaTa clustering at 886 K and ZrZr clustering at 1023 K. 80 The only two discrepancies between the present work and first principles calculations pertain to ZrZr and HfZr. These are likely owing to the difference in the interatomic potential or the cell size used.
C. CSRO effects on shear strength Figure 4 presents bothT al and T af of two RMEPAs with CSRO attained at three different annealing temperatures. The evolution of the strengths in the hard and soft senses is tracked individually. We include for reference the strengths of the A-atom versions of these RMPEAs, which are found to provide a good estimate ofT al , but not T af , with respect to the random alloy.
We first discuss the alias shear strength.T al values for HfMoNbTaTi are higher than those for HfNbTaTiZr, a reflection of the fact that Mo is stronger than Zr. For both quinaries, in the case of A-atom pure metal and random alloy, which have no CSRO,T al on the {112} plane is higher than that on the {110} one. In HfNbTaTiZr, with an increasing CSRO,T al only slightly increases in the hard sense, while decreasing in the soft sense. In HfMoNbTaTi, the influence of CSRO is more noticeable, yet it is not consistent for the two shear planes. For example, when the annealing temperature decreases from 600 to 300 K,T al on the {110} plane decreases, while that on the {112} plane increases. The anisotropy in CSRO effects, exhibited for different senses or planes, is intriguing. It may be a result of the anisotropic clustering/segregation during the annealing. For both quinaries, increasing the CSRO levels tends to increase the asymmetry between the hard and soft senses. Overall, the CSRO effects onT al are relatively weak. On the other hand, prior atomistic simulations showed that T al is positively correlated with the local slip resistance in pure refractory metals 70 and RMPEAs. 36,60 These suggest that the CSRO effects on slip resistance are weak in RMPEAs, consistent with prior atomistic studies. [20][21][22][23][24] Next, we discuss the affine shear strength. T af values for the A-atom pure metals and RMPEAs exhibit the expected anisotropy with the {112} strength in the hard sense being greater than the {110} strength and both greater than the {112} strength in the soft sense. Moreover, compared with T al , T af is less sensitive to the bulk composition. However, both quinaries display similar trends in T af with increasing CSRO level. For example, when the annealing temperature decreases from 900 to 300 K, T af increases in both RMPEAs.
The mechanism of affine shear deformation is not concentrated to a single plane in the crystal and in fact could be unrelated to dislocation motion. Thus, its deformation mechanisms are completely different from that presumed by alias shear. For example, subject to the affine simple shear, both Mo and Nb yield via the homogeneous nucleation of twins, as shown in Fig. 5(a). The same yielding mechanism was reported in tensile loading of a Mo single crystal using the same interatomic potential employed here. 32 In Ta, as well as the two A-atom materials, yield occurs first by a BCC-to-FCC phase transformation (PT) followed by homogeneous nucleation of twins, as shown in Fig. 5(b). A similar BCC-to-FCC PT was reported in shearing of Fe using atomistic simulations. 82 Both RMPEAs, with or without CSRO, yield first by a BCC-to-FCC PT in local regions followed by twinning in those same regions. In HfMoNbTaTi, Ti atoms tend to form clusters and initial local PT takes place at the boundaries of the Ti-rich regions, as shown in Fig. 5(c). Likewise, in HfNbTaTiZr, the HfTi pair tends to segregate and initial local PT occurs in these HfTi-rich regions. Between the two RMPEAs, deformation twinning was experimentally observed in HfNbTaTiZr, 5 but not yet in the other alloy.
CSRO effects on deformation of RMPEAs have been studied previously only for polycrystals using experiments (exp) or MD simulations, where CSRO was found to increase the tensile strength in HfNbTaTiZr (exp) 9 and MoNbTaW (MD) 20 as well as hardness and compressive strength in HfNbTaZr (exp). 83 Although shear deformation was not explicitly studied, results from uniaxial deformation indicate that CSRO likely increases T af . Recent atomistic simulations in FCC CoCrNi MPEA 84  FCC-to-BCC PT, which suggests that it may decrease T af . CSRO was also shown to increase γ usf in two BCC RMPEAs, MoNbTi and TaNbTi, 85 suggesting that it might promote T al . Here, our simulations show that higher levels of CSRO have different effects on T al and T af in RMPEAs. For both strengths, the CSRO effects are not influenced by either the cross-sectional area or temperature, a conclusion that can be reached if one compared Fig. 4 with Fig. S5 in the supplementary material. Finally, for the same sense on the same plane in the same RMPEA, T al is plotted with respect to T af in Fig. S6 in the supplementary material. It is shown that the two strengths are loosely correlated, with higher T al corresponding to higher T af and vice versa. The weak relationship may be expected since the mechanisms accommodating deformation in each type are different, in either pure metals or RMPEAs.

IV. CONCLUSIONS
In this article, we conduct atomistic simulations to calculate the alias (T al ) and affine (T af ) ideal simple shear strengths in two quinary RMPEAs, HfMoNbTaTi and HfNbTaTiZr, with or without CSRO. The strengths in two opposing ⟨111⟩ directions on both {110} and {112} planes are calculated. For reference, the same strengths of the constituent elements Mo, Nb, and Ta, as well as hypothetical A-atom homogeneous versions of HfMoNbTaTi and HfNbTaTiZr, are also calculated. The main findings are summarized as follows: 1. The interatomic potentials for both RMPEAs are validated by comparing the corresponding lattice parameters and elastic constants against those from DFT and experimental data. 2. The size of the cross-sectional area A does not affect the average T al much. However, a smaller A increases the standard deviation in T al . 3. The two RMPEAs possess much lower T al and T af than the constituent BCC metals: Mo, Nb, and Ta. 4. In both RMPEAs, five atomic pairs, NbNb, TaTa, HfNb, HfTa, and NbTa, tend to segregate. The degree of CSRO increases with decreasing annealing temperature. The most likely pairing is TiTi in HfMoNbTaTi and HfTi in HfNbTaTiZr. 5. Subject to affine simple shear deformation, the two RMPEAs yield by two steps: first a BCC-to-FCC PT followed by twinning. In HfMoNbTaTi, the PT preferentially occurs at the boundaries of Ti clusters, whereas in HfNbTaTiZr the PT tends to initiate in the HfTi-rich regions.

SUPPLEMENTARY MATERIAL
See the supplementary material for information as referred to in the main text.