## Abstract

Single-phase solid-solution refractory high-entropy alloys (HEAs) show remarkable mechanical properties, such as their high yield strength and substantial softening resistance at elevated temperatures. Hence, the in-depth study of the deformation behavior for body-centered cubic (BCC) refractory HEAs is a critical issue to explore the uncovered/unique deformation mechanisms. We have investigated the elastic and plastic deformation behaviors of a single BCC NbTaTiV refractory HEA at elevated temperatures using integrated experimental efforts and theoretical calculations. The in situ neutron diffraction results reveal a temperature-dependent elastic anisotropic deformation behavior. The single-crystal elastic moduli and macroscopic Young’s, shear, and bulk moduli were determined from the in situ neutron diffraction, showing great agreement with first-principles calculations, machine learning, and resonant ultrasound spectroscopy results. Furthermore, the edge dislocation–dominant plastic deformation behaviors, which are different from conventional BCC alloys, were quantitatively described by the Williamson-Hall plot profile modeling and high-angle annular dark-field scanning transmission electron microscopy.

## INTRODUCTION

High-entropy alloys (HEAs) are defined as equimolar or nonequimolar single- and multiphase solid solutions, depending on the number and type of principal/secondary metallic elements with the configurational entropy greater than 1.61*R*, (where *R* is the gas constant) (*1*, *2*). The minimization of the Gibbs free energy by increasing the configurational entropy results in the formation of a single- or multiphase solid solution, such as body-centered cubic (BCC) (*3*), face-centered cubic (FCC) (*4*), and/or hexagonal close-packed solid-solution phases (*5*) instead of intermetallic compounds (*6*). Moreover, the random distribution of multiple elements with different atomic radii leads to a severely distorted lattice with variations in the interatomic distances along the crystallographic direction of the unit cells (*3*). These features, in turn, contribute to desirable mechanical properties, such as high hardness, strength, ductility, and softening resistance at room temperature (RT) and elevated temperatures (*3*).

To identify the origin of these desirable mechanical properties, theoretical efforts have been conducted to verify the deformation mechanisms, by considering the elastic and plastic deformations in the BCC refractory HEAs (*7*). In contrast to conventional BCC alloys that have a narrow distribution of interatomic distances, HEAs have a wider variety of bond lengths and atomic interactions in the unit cell, which may play a role in the unique elastic and plastic deformation behaviors at RT and high temperatures.

Recently, Tian *et al*. (*7*) suggested alloy design strategies to develop elastically isotropic HEAs (TiZrNbMoV* _{x}* systems) using state-of-the-art ab initio simulations. According to their calculations, the TiZrNbMo

*V*

_{y}*[where (*

_{x}*x*,

*y*) are (0, 0.9) and (0.5, 0.8)] refractory HEAs exhibit the elastically isotropic deformation behavior, which is strongly correlated to the value of the valence electron concentration (VEC) (

*7*). Furthermore, these previous investigations provide criteria to evaluate the brittle-ductile behaviors based on single-crystal elastic constants, predicting that a few BCC HEA systems can be ductile (

*7*–

*9*). However, the validity of these theoretical calculation methods to design the elastically isotropic and ductile BCC HEAs are still under debate and are lacking experimental validation (

*7*,

*10*), which are critical issues within the HEA community.

To demonstrate these predicted unusual deformation behaviors, we investigated the elastic and plastic deformation behavior for the NbTaTiV BCC refractory HEA at RT and elevated temperatures, using in situ neutron experiments and theoretical approaches. Different from the conventional BCC alloys (*11*), the lattice strain evolution and measured diffraction elastic moduli indicate the elastically isotropic behavior at RT. The magnitudes of single-crystal and macroscopic elastic moduli, obtained by in situ neutron experiments and resonant ultrasound spectroscopy (RUS), respectively, are compared to theoretical calculations, including first-principles calculations and machine learning (ML) approaches. Furthermore, it is found that the present NbTaTiV refractory HEA indicates the lack of strong temperature-dependent elastic anisotropic deformation behavior at elevated temperatures, which is atypical elastic deformation behavior, compared to other conventional metallic materials (*12*–*15*).

The excellent plasticity of the present refractory HEA can be predicted by the relationship between the polycrystalline elastic moduli, such as the Cauchy pressure (*C*_{12} − *C*_{44}) (*16*), Pugh ratio [bulk modulus (*K*)/shear modulus (*G*)] (*8*), and Poisson’s ratio (ν) (*9*). Moreover, we investigated the type of mobile dislocations using the neutron diffraction (ND) peak-broadening modeling and analysis using the Williamson-Hall plot method (*17*), experimentally verified by the high-angle annular dark-field (HAADF) scanning transmission electron microscopy (STEM) technique. It is found that the dominant type of mobile dislocations during plastic deformation could be edge dislocations, which is different from conventional BCC metals. These new discoveries of unique elastic and plastic deformation behaviors in the NbTaTiV refractory HEA will facilitate the development of possible alloy systems for structural material applications.

## RESULTS

### Formation of the single BCC solid solution and the mechanical properties of the homogenization-treated NbTaTiV HEA at elevated temperatures

The microstructure of the homogenized NbTaTiV refractory HEA was characterized, using the scanning electron microscopy (SEM) backscattered electron (BSE) technique, as shown in Fig. 1A. A typical simple solid-solution HEA microstructure with grain sizes of 200 to 400 μm was identified without the presence of secondary phases or elemental segregations. The chemical composition of the sample was verified (using the point energy-dispersive x-ray spectroscopy analysis) to be Nb_{23.8}Ta_{25.5}Ti_{24.9}V_{25.8}, indicating a homogeneous elemental distribution in a nearly equiatomic percentage. The measured ND patterns of the bulk samples, as illustrated in Fig. 1B, confirm the single BCC solid solution at RT and elevated temperatures. The lattice parameters of the present alloy at RT, 500°, 700°, and 900°C, which were determined by a single-peak fitting approach (*18*), are 3.2318, 3.2445, 3.2511, and 3.2588 Å, respectively. The single BCC phase is stable up to 900°C, indicating the excellent phase stability at elevated temperatures. Figure 1C shows the nominal compressive true stress-strain curves, as obtained from the in situ neutron measurement of the homogenization-treated NbTaTiV alloy at RT, 500°, 700°, and 900°C with an estimated strain rate of 1 × 10^{−4} s^{−1}, which were directly obtained from the in situ neutron measurement. The yield strength (σ_{y}) and work-hardening exponent (*N*) of the present alloy gradually decreased from 1064 MPa (RT) to 595 MPa (900°C) and from 0.2732 (RT) to 0.1199 (900°C), respectively. The work-hardening exponent (*N*) was estimated for the different experimental temperatures using the Hollomon’s equation (*19*)* _{t}* is the true applied stress,

*k*is the strengthening coefficient, and ε

*is the true plastic strain. During plastic deformation, the applied load relaxes slightly, while the neutron recording time was set to approximately 20 min, with a constant displacement control mode at 700°C. Moreover, this relaxation behavior becomes more pronounced at 900°C (about 440 to 580 MPa), as compared to 700°C (about 130 to 260 MPa), presented in Fig. 1C.*

_{t}### Lattice strain evolutions of the homogenized NbTaTiV alloy at RT and elevated temperatures

In situ ND with compression tests on the homogenization-treated NbTaTiV alloy were carried out to investigate the elastic deformation behavior for grains with different orientations. Figure 2 shows the corresponding lattice strain (ε* _{hkl}*) evolution along the axial and transverse directions of the {110}, {200}, {211}, and {310} planes, during elastic deformation at various temperatures (RT, 500°, 700°, and 900°C). The lattice strains were obtained from the following equation

*d*is the

_{hkl}*hkl*lattice spacing of the corresponding applied stress, and

*hkl*lattice spacing before loading. The slope of the lattice strain evolution versus applied stress, during elastic deformation, determines the diffraction elastic constants for the specific {

*hkl*} planes. It was found that the variation in diffraction elastic constants among all the grains is not clearly observed, indicating the elastically isotropic deformation behavior, as presented in Fig. 2A.

Unlike the lattice strain evolution at RT, different slopes in the lattice strains versus applied stresses at elevated temperatures are observed due to the increased elastic anisotropy, as presented in Fig. 2 (B to D). At higher temperatures, where the overall stiffness is reduced, the anisotropy increases. The {200} grain exhibits the largest elastic lattice strain with the lowest stiffness, and the {310} grains undergo the second largest elastic strain with a nearly linear response. Whereas the {110} and {211} grains present lower lattice strains and higher directional strength-to-stiffness ratios than other orientations, the diffraction elastic moduli, *E _{hkl}*, of different crystallographic planes at elevated temperatures (listed in table S1) show almost the same values of elastic moduli in different orientations. This trend indicates that the NbTaTiV HEA exhibits the unique elastic isotropic deformation behavior at RT. Moreover, the elastic anisotropy gradually increases at elevated temperatures, but the variation of the increase is not notable, which can be distinguished from other conventional alloys (

*12*–

*14*).

### Calculations and measurements of single-crystal elastic moduli for the NbTaTiV alloy via the Kroner model and RUS technique

To better understand the underlying elastic deformation behavior of the NbTaTiV alloy, as a function of temperature, the macroscopic elastic, shear, and bulk moduli (*E _{M}*,

*G*, and

_{M}*K*) are obtained from the diffraction elastic moduli (

_{M}*E*) and Poisson’s ratio (ν) by applying theoretical models. The three most widely applied theoretical models to determine the macroscopic elastic moduli are the Voigt, Reuss, and Kroner models, respectively (

_{hkl}*20*,

*21*).

Figure 2 (E to H) shows the reciprocal diffraction elastic moduli (1/*E _{hkl}* and ν

_{hkl}/

*E*) as a function of the elastic anisotropy factor,

_{hkl}*A*, [

_{hkl}*20*,

*21*). The results, as exhibited in Fig. 2 (E to H), indicate that the Kroner model is in the best agreement with the experimental data for all temperature ranges, as compared to the Voigt and Reuss models. Note that the absolute percentage SEs of Kroner model fitting are 2.46, 6.34, 3.15, and 8.85% for RT, 500°, 700°, and 900°C, respectively, showing the acceptable accuracy of the estimated single-crystal elastic constants. These two theoretical models only consider homogeneous strains (Voigt) or stresses (Reuss) in all the grains during continuous loading. On the other hand, the Kroner model accounts for the variation of stresses and strains among all the grains. Thus, the Kroner’s self-consistent model, which is the most suitable applied calculation method to determine the macroscopic elastic moduli from diffraction elastic moduli for HEA systems (

*20*,

*21*), was applied to the lattice strain evolution results to obtain the macroscopic elastic moduli. The detailed calculation processes for the Kroner model are described in section S1.

The calculated single-crystal elastic constants, *C _{ij}*, macroscopic elastic (

*E*), shear (

_{M}*G*), bulk (

_{M}*K*) moduli, and Poisson’s ratio (ν) at RT and elevated temperatures are listed in Table 1. As the temperature increases up to 700°C, the variation in

_{M}*C*

_{11}is more pronounced, as compared with

*C*

_{12}and

*C*

_{44}. These results suggest that a longitudinal strain notably induces a change in volume. This change in volume is strongly associated with the temperature and, thus, corresponds to a relatively large change in

*C*

_{11}. Consequently, the tetragonal shear modulus,

*C′*= [(

*C*

_{11}−

*C*

_{12})/2], is substantially reduced with the gradually decreasing elastic, shear, and bulk moduli. Moreover, compared with the elastic constants at 700°C, the values of

*C*

_{11}and

*C*

_{12}are drastically reduced at 900°C, while

*C*

_{44}remains nearly constant at RT, 500°, 700°, and 900°C, resulting in a small change in the tetragonal shear modulus,

*C′*. This result implies the elastic stability of this BCC HEA against the tetragonal shear deformation.

Further analysis of the elastic moduli was performed using the RUS. The elastic constants, *C*_{ij}, were determined from ~45 resonances where the frequencies ranged from 100 to 650 kHz using an iteratively inverse fitting process of the measured resonances. A typical RUS resonance spectrum obtained from the NbTaTiV refractory HEA is shown in fig. S1, in which the inset presents an example of Lorentz fitting of a single resonance peak. The obtained values of the elastic moduli for the NbTaTiV alloy are summarized in Table 1, showing good agreement with the in situ neutron data.

### Predictions of elastic moduli for the NbTaTiV refractory HEA via first-principles calculations

On the basis of the experimental observations of the elastic deformation behavior for the NbTaTiV alloy, first-principles calculations are used to theoretically confirm the values of single-crystal elastic moduli at RT using the density functional theory (DFT). The Young’s, bulk, and shear moduli at −273°C and RT for the current HEA are quantitatively calculated from the elastic constants determined by fitting ab initio–calculated stresses that correspond to the strained structures. The calculated and measured bulk elastic moduli at the absolute zero temperature and RT are summarized in Table 1, which show good agreement between experimental and theoretical results. The calculated elastic moduli at RT were less than the experimentally measured values by approximately −15.9, −16.6, and −6.0% for the Young’s, shear, and bulk moduli, respectively. The above result corroborates the general trend of decreased stiffness in metals with increasing temperature. Furthermore, the calculated elastic constant values for *C*_{11}, *C*_{12}, and *C′* all decrease from absolute zero to RT, while *C*_{44} slightly increases from 31.6 to 33.5 GPa.

The elastic properties of the NbTaTiV HEA were analyzed and plotted with respect to the crystallographic orientation, as illustrated in Fig. 3. To compute the direction-dependent elastic properties, we applied tensorial operations to the fourth-order elastic tensor determined from the DFT calculations. The fourth-order elasticity tensor is transformed in a new basis set to determine elastic properties in any direction using the equation below*22*).

A higher degree of isotropic elastic properties results in the smaller variation in the magnitude of the three-dimensional (3D) surfaces and 2D contours. The 2D contour plots correspond to a cut through the 3D surface plots on the *XY*, *XZ*, and *YZ* planes. The Young’s modulus is represented by a solid surface with a color map in which darker colors denote smaller values, and brighter colors indicate larger values. The maximum-to-minimum Young’s modulus ratio is 1.16, indicating little anisotropy characteristic. The maximum value of 107.7 GPa corresponds to the <100> direction, while the minimum value for the <111> direction is 93.0 GPa (see Fig. 3, A and B). The 3D surface plot for the linear compressibility is represented by a translucent green sphere, which exhibits negligible anisotropy because it has a value of 2.42 TPa^{−1} in all directions (Fig. 3, C and D). The range of the minimum and maximum values of the shear moduli and Poisson’s ratios for each direction of the applied stress is attributed to the variation in the direction of the measurement. In the 3D surface plots of the shear modulus and Poisson’s ratio, the translucent green surface represents the maximal values and envelops the solid surface inside, which represents the minimal values. The minimum surface has a color map applied such that darker colors represent smaller values than brighter colors, as illustrated in Fig. 3 (E and G). In the 2D plots, the maximum values are represented by a green curve, and the minimum values are denoted by a pink curve. The ratio of the maximum to the minimum shear modulus values is 1.17. The maximum is in the <111> direction with a value of 39.33 GPa, whereas the minimum corresponds to the <110> directions with a value of 33.51 GPa. Furthermore, the ratio of the maximum Poisson’s ratio to the minimum ratio is 1.32. The maximum and minimum values both occur in the <110> direction with values of 0.44 and 0.33, respectively (Fig. 3, F and H).

Recently, Ranganathan and Ostoja-Starzewski (*23*) introduced a measurement method of anisotropy using the universal elastic anisotropy index, which can generally be applied to crystal systems that consider the bulk component of the elastic tensor. As opposed to other anisotropy measurement methods, such as the Zener anisotropy index, which is only applicable to cubic crystals, and the Ledbetter-Migliori anisotropy index, which only considers the shear component of the elastic tensor, the universal anisotropy index is given by the following equation*G* is the shear modulus, *K* is the bulk modulus, and the *V* and *R* superscripts refer to the Voigt and Reuss estimates, respectively. As a comparison, the *A ^{U}* values of the constituent elements obtained from the Materials Project database are Nb, 2.71; Ta, 0.07; Ti, 0.27; and V, 3.42. As a result, the

*A*value for the NbTaTiV HEA is 0, which defines the single-crystal elastic properties to be isotropic.

^{U}### Predictions of elastic moduli for the NbTaTiV refractory HEA via ML models

The set of elemental and structural properties that were found to be good predictions in the shear modulus ML model were the Holder mean (fourth power) of the group numbers of the elements, cohesive energy, the Holder mean (second power) of the atomic radii of the elements, the Holder mean (fourth power) of the electronegativity of the elements, and density. The set of elemental and structural properties that were found to be good predictions in the bulk modulus ML model were the Holder mean (fourth power) of the group numbers of the elements, cohesive energy, the Holder mean (first power) of the atomic radii of the elements, Holder mean (negative fourth power) of the electronegativity of the elements, and density. The ML models’ performance was tested on 20% of the dataset, which was never seen by the ML model during the fitting process. The bulk and shear moduli were log normalized for model training and testing to ensure that large outlier values do not have an exaggerated effect on the models. The model performance metric used was the mean squared error. The mean squared errors of the shear modulus and bulk modulus models were 0.0409 and 0.0284, respectively. The ML model predictions, compared with experimentally and DFT-calculated values, are summarized in Table 1. The predicted shear modulus is 36.6 GPa, and the predicted bulk modulus is 146.3 GPa, giving relative errors of 2.5 and 6%, respectively.

There is also the stochastic variability in the model performance, depending on how the training dataset is shuffled. To determine the variability of model performance that is simply due to which training examples are included in the training dataset, we plot the learning curves in Fig. 4. The variability in model predictions due to the random shuffling of data is studied by randomly shuffling the data in 20 different ways. The red curves represent the cross-validation mean squared error of the models during training averaged over the 20 different trials, and the green curves denote the mean squared error of the models calculated on the testing set averaged over the 20 different trials. The shaded areas represent 1 standard deviation (SD) above and below, as determined from the 20 different trials. The *x* axis of the learning curves shows how the size of the training dataset itself affects the model performance. The learning curves indicate that the model prediction performance on the testing set converges after 3000 training samples. The small difference in the training and testing mean squared errors indicates that the model is sufficiently good at predicting the shear and bulk moduli for compounds that were not in the training dataset. This trend also illustrates the generalizability of the model in the sense that it can be used on compounds, which are not included in the original dataset.

## DISCUSSION

### Elastic isotropy and anisotropy deformation behaviors of the NbTaTiV refractory HEA

Typically, the response of lattice strains under applied stresses is strongly associated with the oriented grains due to the effects of elastic anisotropies in polycrystalline metal materials (*11*, *13*), including FCC-structured HEAs (*24*). The elastically anisotropic deformations originate from the orientation-dependent elastic modulus, which is presented by the different slopes in the plots of the lattice strain evolution as a function of applied stresses. In the case of BCC-structured materials, the {110} and {211} grains, which have higher elastic stiffness than other grains, exhibit the larger slopes. The {200} and {310} grains, on the other hand, have smaller slopes with lower elastic stiffness (*11*, *25*). After this elastic stretch, the sequence of grain yielding depends on the value of the Schmid factor, showing varying stress levels for the onset of yielding in differently grains, which is commonly characterized as the grain-to-grain yielding behaviors (*25*). The plastic anisotropy is derived from the occurrence of plastic deformation on certain (*hkl*) grains, on which slip is preferred to occur. An increase in the slope indicates plastic yielding of particular grains, whereas the reduction of the slope represents an intergranular load transfer from the yielding grains to the neighboring grains, which are still in the elastic deformation region. With these lattice strain evolution trends, {110} and {211} grains with the largest stiffness in the elastic region tend to yield at the early stage due to their slip system. Subsequently, a load transfers from the aforementioned grains to the {200} and {310} grains rather than the occurrence of load transfer in various orientations (*11*, *13*, *26*). This inharmonic load partitioning could lead to stress concentrations on particular grains ({200} grains in BCC structured alloys) after a large amount of plastic deformation. Unlike conventional BCC alloys, as discussed above, the NbTaTiV BCC refractory HEA exhibits no orientation-dependent lattice strain evolutions during elastic deformation at RT, leading to nearly similar values of elastic moduli in differently grains, as presented in Fig. 2A. With this unusual elastic deformation, all orientated grains yield at the same stress level (around 1000 MPa), as presented in fig. S2. Furthermore, the features of the load transfer were not clearly presented during plastic deformation at RT, indicating the reduced plastic anisotropy. This trend of the reduced plastic anisotropic deformation behavior could result in the excellent ductility (plasticity) at RT because of the homogeneous load sharing between differently grains rather than the load concentration on particular grains during plastic deformation, as shown in fig. S2.

The elastic isotropy is often quantitatively described by the Zener anisotropy (*A _{z}*) ratio and Voigt and Reuss bounds (

*A*)

_{VR}*G*= [5(

_{R}*C*

_{11}−

*C*

_{12})

*C*

_{44}]/[4

*C*

_{44}+ 3(

*C*

_{11}−

*C*

_{12})] and

*G*= (

_{V}*C*

_{11}−

*C*

_{12}+ 3

*C*

_{44})/5 (

*7*). The calculated values of

*A*and

_{Z}*A*at elevated temperatures can be found in table S2. In the case of ideal isotropic materials, the value of the tetragonal shear modulus (

_{VR}*C′*) will be identical to the cubic shear modulus (

*C*

_{44}), and it results in

*A*= 1 and

_{z}*A*= 0 (

_{VR}*7*,

*10*). Typically, the Zener ratios for most of the conventional FCC and BCC structural materials are well over 2.0 or 3.0, indicating the substantial elastic anisotropy (

*12*,

*24*). Recently, Tian

*et al*. (

*7*) suggested an alloy design strategy to achieve the elastic isotropy in HEAs using ab initio calculations. It was reported that the Mo-containing BCC refractory HEAs have the same polycrystalline shear and Young’s moduli along all orientations, indicating elastic isotropy. Furthermore, these isotropic HEAs have a similar VEC of about 4.7 (in general, HEAs tend to form BCC solid-solution phases when VEC becomes less than 7.55) (

*7*,

*27*), which is similar to other elastic isotropy alloys, such as the Ti-V alloy (gum metals) (

*28*). However, their reported works have strongly focused on only theoretical calculations, due to the lack of experimental verifications. According to our calculated results of

*A*and

_{Z}*A*in table S2, which are primarily obtained from the in situ neutron experiments, the NbTaTiV refractory HEA manifests almost perfectly isotropic elastic deformation at RT (

_{VR}*A*≈ 1.2,

_{Z}*A*≈ 0.0056, and VEC ≈ 4.68). Moreover, the elastic isotropy was gradually increased at elevated temperatures. Many conventional alloys show the strong temperature-dependent elastic anisotropic deformation behavior, indicating a substantial increase of the Zener anisotropy parameters at higher temperatures (

_{VR}*11*–

*13*). For example, the nanostructured BCC ferritic alloy (NFA) and FCC-structured Ni-based superalloy exhibited the large variation of anisotropy parameters from RT to high temperatures [from 1.7 (RT) to 5.3 (800°C) for NFA and from 2.8 (RT) to 3.6 (700°C) for a Ni-based superalloy] (

*12*,

*13*). Furthermore, the pure iron and Fe

_{90}Cr

_{10}random alloy are also expected to have a notable increase of anisotropy parameter at elevated temperatures (

*14*). It was further found that the CoCrFeMnNi HEA exhibited the strong temperature-dependent elastic and plastic anisotropic deformation from RT to 527°C (

*15*). In case of the present NbTaTiV HEA, the Zener anisotropy parameters were gradually increased at elevated temperatures, but the variation of anisotropy parameters in NbTaTiV (1.2 at RT and 1.6 at 900°C) is much smaller than other conventional alloys (

*11*–

*13*). The detailed comparisons of the Zener anisotropy parameters and yield strengths at elevated temperatures between the NFA (

*13*) and the present NbTaTiV refractory HEA are indicated in fig. S3, manifesting that both alloys show the temperature-dependent elastic anisotropy and the smaller variation of the Zener anisotropy parameters in the NbTaTiV HEA than the nanostructured ferritic alloy at elevated temperatures.

### Mechanical stability and ductility (plasticity) of the NbTaTiV refractory HEA

The mechanical stability of the NbTaTiV refractory HEA can be evaluated, considering the values of single-crystal elastic constants. Born and Huang (*29*) suggested the evaluation criteria of the elastic stability through the dynamical stability conditions, *C*_{11} > ∣*C*_{12}∣, *C*_{11} + 2*C*_{12} > 0, and *C*_{44} > 0 (*30*). According to the elastic constants in Table 1, the present alloy shows good mechanical stability in all temperature ranges.

The intrinsic brittle-to-ductile deformation behavior of HEAs, which dominantly consists of solid-solution phases, can be predicted by the interrelation of polycrystalline elastic moduli, such as the Cauchy pressure (*C*_{12} − *C*_{44}) (*16*), Pugh ratio [bulk modulus (*K*)/shear modulus (*G*)] (*8*), and Poisson’s ratio (ν) (*9*). The Cauchy pressure is often used to categorize the type of atomic bonding. The relatively large positive Cauchy pressure value represents an increased degree of the metallic bonding character with the ductile deformation behavior, while the negative Cauchy pressure stands for covalent bonding with brittle behavior (*16*). Moreover, the *K*/*G* ratio, which has been suggested by Pugh, can be applied to predict the ductility and brittleness of materials (*8*). The bulk modulus (*K*) assesses the resistance to fracture, and the shear modulus (*G*) is defined as the resistance to the plastic deformation of polycrystalline materials. When the alloys have large values of *K* with relatively small values of *G*, alloys tend to exhibit the considerable fracture resistance after the onset of plastic deformation, leading to the ductile characteristics of materials (*8*). From existing experimental results for a large number of studied materials, it was found that the brittle-to-ductile transition point (*K*/*G*)_{0} is 1.75 (*30*). Hence, higher values of Pugh ratios (>1.75) indicate that the materials are ductile, whereas lower values are related to the brittle property. On the basis of the empirical criteria of the Pugh ratio above, the brittle-to-ductile limit can also be assessed by the Poisson’s ratio, considering the relationship between ν and (*K*/*G*)_{0}, as ν = (3*K* − 2*G*)/[2(3*K* + *G*)], in the case of isotropic polycrystalline materials. Therefore, the brittle-to-ductile transition limit, (*K*/*G*)_{0} > 1.75, corresponds to ν > 0.26. Note that the brittle-ductile criterion for bulk metallic glasses is ν > 0.31 (*9*). The values for the Cauchy pressure, Pugh ratio, and Poisson’s ratio (see Table 1 and table S2), which are obtained from in situ neutron experiments, suggest that NbTaTiV is intrinsically ductile.

### Characterization of mobile dislocations for the NbTaTiV alloy at elevated temperatures

To further investigate the plastic deformation behaviors of the NbTaTiV alloy, the line-profile modeling and analyses, considering the evolution of the ND peak width, are carefully performed as a function of the macrostrain at elevated temperatures. In general, plastic deformation often leads to substantial broadening of diffraction peaks, which is mostly caused by the formation of inhomogeneous strain fields (*31*). It is well known that the inhomogeneity of the subgrain structure is usually more related to the accumulation of dislocations than variations of intergranular strains. Thus, it gives rise to strain fluctuations with the strain field around dislocations responsible for diffraction peak broadening (*32*). With the above fundamental understanding of the plastic deformation behavior in metal materials, evolutions of diffraction peak widths at different strain levels for the NbTaTiV refractory HEA are further discussed by the Williamson-Hall peak-broadening analysis approach (*17*). The evolution of the peak widths of the full width at half maximum (FWHM/FWHM_{0}) of {110}, {200}, {211}, and {310} diffraction peaks as a function of macrostrain at RT, 700°, and 900°C, and the detailed processes of the Williamson-Hall peak-broadening analysis during plastic deformation are described in section S2.

Figure 5 shows the modified Williamson-Hall plots for the NbTaTiV alloy corresponding to the macrostrain range from 0 to 15%. The response of linear fitting of the plot, which is obtained from the assumption of specific types of dislocations, refers to the dominant mobile dislocation character during plastic deformation (*33*). Hence, when treating *q* and *C*_{h00} values as pure screw dislocations (*q* = 2.2244 and *C*_{h00} = 0.2415), the nonmonotonous deviation of *∆K* as a function of *K*^{2}*C* gradually increases with an increase in the macrostrain, as exhibited in Fig. 5A. Nevertheless, the data follow a linear trend at all strain levels, when the values of *∆K* and *K*^{2}*C* obtained with the pure edge dislocation consideration (*q* = 0.0673 and *C*_{h00} = 0.1937), as exhibited in Fig. 5B. The deviations from the linear-fitted line are quantitatively described by *R*-squared coefficients in table S4. These results clearly indicate that the type of mobile dislocations of the NbTaTiV refractory HEA is edge dislocations in the initial condition (nondeformed) at elevated temperatures. Moreover, with increasing macrostrains (up to 15% deformed), the <111> {110} edge dislocation remains as a dominant mobile dislocation character at RT and elevated temperatures (Fig. 5, B to D).

The BCC metallic materials typically contain edge dislocations, which move faster than screw dislocations because of their core properties of dislocations during plastic deformation (*34*). The critical stress, which forces the edge dislocation movement, is quite small (<10 MPa for the pure BCC Fe), leading to the easy annihilation of the edge dislocations, which results in the screw/edge anisotropy in the apparent mobility deviation (*34*). Therefore, the remaining dislocations tend to have more screw character in BCC metals. Recently, however, Rao *et al*. (*35*) have reported atomistic simulations of plastic deformation behavior for the BCC HEA in terms of dislocation motion. This in-depth investigation predicts the tendency of dislocation motions in the BCC solid-solution HEA to correspond to the random elemental distribution, which is considerably different from conventional BCC metals. During the plastic deformation, the edge dislocation core can be deviated from their neutral dislocation planes, producing the wavy dislocation lines, which are attributed to the severe lattice distortion. This feature is induced by compositional fluctuations and, thus, leads to not only strengthening but also notable difficulty in the edge dislocation movement. These predictions of plastic deformation behavior for the BCC HEA suggest that the number of edge dislocations is comparable to that of the screw dislocations, which is in contrast to the conventional low-entropy BCC metals. Similar to the above reported prediction, the dominant mobile dislocations of the investigated NbTaTiV BCC refractory HEA, which has severely distorted lattices, are close to <111> {110} edge dislocations with increasing macrostrain, as presented in Fig 5B.

To further identify and provide visible evidence of the type of mobile dislocations, during plastic deformation in the NbTaTiV HEA, the nanoscale microstructure of the deformed sample was carefully investigated with the TEM study. Figure 5E exhibits the filtered HAADF-STEM image of the 15% RT-deformed NbTaTiV HEA, which are acquired along the [110] zone axis. The formation of edge dislocations is clearly observed with extra planes on the image, described by the symbol, “┴,” which is consistent with line-profile modeling, using Williamson-Hall plot results (Fig. 5B). It is found that these edge dislocations were formed on different planes, indicated by various colored symbols of edge dislocations. The high-magnification image of Fig. 5E is illustrated in Fig. 5F to provide a clear description of the edge dislocations.

Overall, this comprehensive study provides evidence for the unique elastic and plastic deformation behavior of the single BCC solid-solution phase refractory HEA at elevated temperatures using in situ neutron experiments coupled with theoretical calculations. The results of the ND and characterization of the microstructure show the excellent BCC single-phase stability up to 900°C. Furthermore, the present alloy has the high yield strength and excellent plasticity at RT and softening resistance at elevated temperatures. The lattice strain evolution results show the distinctive elastic deformation behavior as an elastic isotropy at RT and less temperature dependency in the elastic anisotropic deformation behavior than conventional metallic materials at elevated temperatures.

The single-crystal elastic moduli are obtained from the in situ ND data coupled with the Kroner model approach, showing excellent agreement with the experimental values determined from RUS measurements and first-principles calculations. Furthermore, the good agreement of the ML model predictions of shear and bulk moduli is somewhat unexpected, given that the Materials Project database only contains ordered compounds and does not contain disordered alloys, such as HEAs. The lattice distortion effect may also cause the ML model to produce less accurate results for the HEAs, but given the lack of data for comparison, the good agreement suggests that these models may be useful for conservative estimates.

The NbTaTiV alloy is predicted to be intrinsically ductile, based on the ductility criteria, using the Cauchy pressure, Pugh ratio, and Poisson’s ratio. The verification of the mobile dislocation type was quantitatively investigated using the Williamson-Hall plot profile, experimentally verified by the HAADF-STEM. The dominant mobile dislocation is determined to be the character of edge dislocations during the 15% plastic deformation at elevated temperatures. It is thought that the severely distorted lattice induces the deviation of the formation of dislocations from their neutral dislocation planes, which, consequently, leads to a substantial reduction in the mobility of edge dislocations during plastic deformation. Thus, the edge dislocations are believed to be the dominant dislocation type for the BCC refractory HEA, which is atypical, as compared to the conventional BCC metals and alloys. These new discoveries of unique elastic and plastic deformation behaviors of the BCC NbTaTiV refractory HEA may lead to the overall excellent mechanical performance and also provide a road map toward the discovery and production of newly designed polycrystalline materials for structural materials applications.

## MATERIALS AND METHODS

### Materials preparation and microstructural characterization

The alloys were synthesized by arc melting under an argon atmosphere on a water-cooled Cu hearth from Nb, Ta, Ti, and V elements of 99.99 weight % purity. The nominal composition of the present alloy is Nb_{25}Ta_{25}Ti_{25}V_{25} in atomic %. The ingot was melted over 10 times to ensure a fully homogeneous distribution of elements. From the master alloys, the specimens were fabricated, followed by direct casting into cylindrical rods with a 4-mm diameter and 50-mm length using a drop-casting technique. These alloys were sealed with the triple-pumped argon in quartz tubes and homogenized at 1200°C for 3 days, followed by water cooling. The microstructure was examined by SEM using a Zeiss Auriga 40 equipped with BSE. The chemical analysis of the alloy was studied by the energy-dispersive x-ray spectroscopy after the homogenization treatment.

### In situ ND experiments

The in situ ND measurements were carried out on the VULCAN Engineering Diffractometer at the Spallation Neutron Source facility located at the Oak Ridge National Laboratory (ORNL). The ND instrument uses the time-of-flight measurements, which allows for the development of the ND pattern, covering a wide range of d-spacings without the rotation of detectors or samples. The two detectors, which are designated as banks 1 and 2 at ±90°, were used to collect diffracted neutrons from the polycrystalline grains with lattice planes parallel to the axial and transverse directions, respectively. The compression experiments on the homogenization-treated NbTaTiV alloy with a diameter of 4 mm and a length of 8 mm were conducted to investigate the elastic and plastic deformation behaviors at RT and elevated temperatures using a mechanical testing system load frame. The samples are illuminated by the incident neutron beam of a 3 × 3–mm^{2} slit size and 2-mm receiving collimators. During the elastic deformation period, a stepwise force control sequence was used by collecting the 20-min measurement time of ND data at each stress level. When the stress level reached 1030 MPa for RT, 650 MPa for 500°C, 585 MPa for 700°C, and 550 MPa for 900°C (close to the macroscopic yield strength), a stepwise displacement control with an incremental step of 0.2 mm was used. The collected data were analyzed by single-peak fitting, using the VULCAN Data Reduction and Interactive Visualization software (VDRIVE) program (*18*).

### First-principles calculations

First-principles calculations were performed with the Vienna Ab-initio Simulation Package (VASP) using the projector augmented-wave (PAW) method (*36*). The exchange-correlation energy was described with the generalized gradient approximation in the Perdew-Becke-Ernzehof parameterization. A plane-wave cutoff of 700 eV and Monkhorst-Pack *k*-point grid of 4 × 4 × 4 were used for all calculations. The chemical disorder of the solid solution was modeled using a 64-atom special quasirandom structure (SQS). The SQS was used through a method based on a Monte Carlo–simulated annealing loop where atomic configurations are proposed by randomly swapping atomic positions and accepting changes when the configuration’s calculated correlation functions are a closer match to that of a disordered state.

The elastic tensor of the SQS supercell was calculated using a set of distorted structures based on the methodology described in (*37*). The supercell’s lattice parameters were fixed at the experimentally determined RT values, and the ionic positions were allowed to relax. A set of deformed structures were generated by adopting 3 × 3 Green-Lagrange strain tensors of varying magnitudes at ±0.5 and ±1%. For each distorted structure, the 3 × 3 stress tensor is computed by the DFT. A linear least-squares fit between the stress and strain tensors was performed to calculate the elastic tensor of the SQS. However, because the SQS supercell does not conserve the point-group symmetry, a projection technique is used to derive the closest elastic tensors with a cubic symmetry for the NbTaTiV alloy, which is equivalent to crystallographically averaging tensor elements. The projection method is computationally efficient, and the elastic constant values converge quickly with the size of SQS. The elastic moduli of the current HEA are calculated using the elements of the projected elastic tensor, as described below. The bulk modulus (*K*) and shear modulus (*G*) were obtained from the elastic constant. Then, Young’s modulus was mathematically calculated, substituting *K* and *G* in the simple equation *E*, *K*, and *G**s _{ij}* is the stiffness tensor, and the elastic tensor,

*C*, was calculated, following their correlation,

_{ij}### ML models for bulk and shear moduli

ML models using the gradient-boosted trees algorithm were developed to predict shear and bulk moduli (one model for each). Models were built using the scikit-learn library in Python (*38*). The ML models were trained on DFT-calculated shear and bulk moduli, elemental properties, and compound structural properties from the Materials Project database (*37*). There were 6826 inorganic compounds in the dataset used for training and testing the ML models. The most general approach to combining elemental properties for each compound was taken by computing the compositionally weighted Holder means with a range of powers from −4 to 4. The Holder means are defined in the equation below* _{p}* is the Holder mean with the power,

*p*; and

*w*is the compositional fraction of the element,

_{i}*i*, in a compound with

*n*elements. A Holder mean with the powers of 1, −1, or 2 corresponds to the arithmetic mean, harmonic mean, or Euclidean mean, respectively (

*39*). A feature selection methodology based on a genetic algorithmic optimization scheme was used to determine which Holder means of which elemental properties and structural properties result in the best predictive models for shear and bulk moduli. The models were trained using 80% of the dataset with a fivefold cross-validation scheme to prevent overfitting. The model performance is determined on the remaining 20% of the dataset.

### Resonant ultrasound spectroscopy

Elastic constants were also measured using RUS. A cylindrical sample of the NbTaTiV HEA was corner-mounted between two transducers. In this configuration, one transducer excites the sample, while a receiving transducer detects and relays the sample’s response of resonant frequencies. Assuming that the elastic isotropy and the elastic constants, (*C*_{11}, *C*_{12}, and *C*_{44}), were determined, using sample dimensions, and observed resonances. In our study, the frequency was swept through the 100- to 650-kHz range with large responses detected when the frequency corresponds to an eigenfrequency of the sample. From the elastic constants, bulk modulus (*K*) and shear modulus (*G*) were determined using Hill’s approximations (*40*).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/37/eaaz4748/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We would like to thank J. Brechtl, H. Choo, and Y. Chen for assistance with the experiments, respectively. Neither the U.S. Government nor any agency thereof, nor LRST, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

**Funding:**The research was supported by the U.S. Army Office Project (W911NF-13-1-0438 and W911NF-19-2-0049) with the program managers, M. P. Bakas, D. M. Stepp, and S. Mathaudhu, and the material was developed. P.K.L. also thanks the support from the NSF (DMR-1611180 and 1809640) with the program directors, J. Yang, G. Shiflet, and D. Farkas, and the neutron work was conducted. C.L. and B.L.M. would like to acknowledge the partial support from the Center of Materials Processing with C. Rawn as the director, a Tennessee Higher Education Commission (THEC) Center of Excellence located at The University of Tennessee, Knoxville. A portion of research at the ORNL’s Spallation Neutron Source was sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy. Reference herein to any specific commercial product, process, or service by the trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof. Y.-C.C. thanks the funding from the Ministry of Science and Technology (MOST) of Taiwan under grant no. MOST-109-2636-M-009-002, the core facility support at the National Chiao Tung University (NCTU) from the MOST. Y.-C.C. thanks the partial support from the “Center for the Semiconductor Technology Research” from The Featured Areas Research Center Program within the framework of the Higher Education Sprout Project by the Ministry of Education (MOE) in Taiwan and supported, in part, by the Ministry of Science and Technology, Taiwan, under grant MOST 109-2634-F-009-029. M.C.G. acknowledges the support from the U.S. Department of Energy’s Fossil Energy Cross-Cutting Technologies Program at the National Energy Technology Laboratory (NETL) under the RSS contract no. 89243318CFE000003. Furthermore, the present work was supported by the Basic Research Laboratory Program through the Ministry of Education of the Republic of Korea (2019R1A4A1026125) and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (no. 2020R1C1C1005553). W.C. and G.K. acknowledge the support from the NSF under grant nos. OAC-1940114 and DMR-1945380. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under contract no. DE-AC02-05CH11231. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant no. ACI-1548562. This work was partially funded by the Department of Energy, National Energy Technology Laboratory, an agency of the U.S. Government, through a support contract with Leidos Research Support Team (LRST).

**Author contributions:**All authors contributed extensively to the work presented in this manuscript. C.L., P.K.L., K.A., and G.S. performed the SEM, in situ neutron, and mechanical tests and wrote the main manuscripts and the Supplementary Materials. M.C.G. designed the HEA, conducted the calculation of phase diagram, and wrote the computational prediction parts in the main manuscript. G.K. and W.C. conducted the first-principles calculations and ML and wrote the calculation parts in the main manuscript. Y.C. and Y.-C.C. performed the HAADF-STEM and wrote the atomic-scale structure investigation parts in the main manuscript. B.L.M. and V.K. conducted the RUS measurement and wrote the corresponded parts in the main manuscript. All authors discussed the results and implications and commented on the manuscript at all stages.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).