High energy barriers for edge dislocation motion in body-centered cubic high entropy alloys

Recent theory proposes that edge dislocations in random body-centered cubic (BCC) high entropy alloys have high barriers for motion, conveying high strengths up to high temperatures. Here, the energy barriers for edge motion are computed for two model alloys, NbTaV and MoNbTaW as represented by interatomic potentials, using the Nudged Elastic Band method and compared to theoretical predictions. The average magnitude of the barriers and the average spacing of the barriers along the glide direction agree well with the analytical theory, with no adjustable parameters. The evolution of the barriers versus applied stress is modeled, and the mean strength is in reasonable agreement with the predicted zero-temperature strength. These findings validate the analytic theory. A reduced analytic model based on solute misfit volumes is then applied to Hf-Mo-Nb-Ta-Ti-Zr and Mo-Nb-Ta-Ti-V-W alloys, rationalizing the observed significant strength increases at room temperature and 1000 ∘C upon addition of solutes with large misfit into a base alloy. The analytic theory for edge motion is thus a powerful validated tool for guiding alloy selection.


INTRODUCTION
High-entropy alloys (HEAs) are multicomponent alloys with nondilute concentrations of most or all of the alloying elements. Various HEAs have impressive mechanical properties such as high yield strength at room temperature, high ultimate strength, high ductility, high fracture toughness, or high strength retention at very high temperatures [1][2][3][4][5][6] . Of primary interest here is the high strength retention in the refractory BCC HEAs such as MoNbTaW and MoNbTaVW. High strength retention might suggest a solute drag mechanism, but the high strengths starting from low T and the high vacancy formation (see Table 2) and migration energies preclude standard solute drag at typical experimental strain rates (estimates can be made using theories such as that found in ref. 7 ). The high strength and strength retention at high T thus requires that the barriers to dislocation glide are substantial. Recent theory 8 has predicted that edge dislocations in these BCC HEAs encounter such high barriers due to the random fluctuations in local environments in the complex alloy. This prediction is unusual because it is commonly assumed that screw dislocations control plastic flow in BCC alloys. However, simulations of edge dislocation motion on model alloys at T = 0 K show high strengths comparable to those predicted by the parameter-free theory. Application of the edge dislocation model to real alloys also shows good agreement for the strength versus temperature, including the observed high strength retention at temperatures up to 1600 ∘ C. Recent experiments support the key role of edge dislocations in strengthening of some BCC HEAs 9 . The theory makes other predictions that have been validated by atomistic simulations. Here, we report detailed transition-state computations of the atomistic barriers for edge dislocation glide in two model alloys, NbTaV and MoNbTaW, as represented using interatomic potentials that have been shown to be accurate for this problem, and show broad agreement with the theory. We then apply a reduced analytic theory to several further HEAs recently reported in the literature. We show that the significant strengthening due to the addition or substitution of Mo into HfNbTaTiZr and of V into MoNbTaW and MoNbTaTiW can be quantitatively captured using the theory for edge dislocations. The loss of strength in the former alloys at 1200 ∘ C is correlated with the estimated vacancy formation energy, but no mechanism has yet been identified. Overall, the analytic theory based on large barriers to edge dislocation motion, validated here, provides a path for preliminary alloy selection for high strength and high strength retention at high temperatures.

RESULTS AND DISCUSSION Solute strengthening theory
We first briefly review the theory for yield strength for edge dislocations in a random BCC alloy, as presented in ref. 8 . A dislocation in a random alloy minimizes its total energy by becoming wavy over some characteristic length scales ζ c (wavelength~5.6ζ c ) and w c /2 (amplitude). The waviness forms because the dislocation exists in a rough potential energy landscape and so finds local regions of favorable random concentration fluctuations that lower the dislocation energy, but at the energetic cost of an increased dislocation line length. In the wavy structure, the dislocation resides in lower energy in segments of length ζ c and faces barriers created by unfavorable regions at distance w c along the glide direction. Stress and/or thermal activation are needed to overcome these energy barriers and generate plastic flow by dislocation glide. The detailed derivations in the theory are given in ref. 8 and only the relevant results are stated here.
The theory considers an N-component random alloy composed of atom types n = 1, . . . N each at concentration c n . The solute/ dislocation interaction energy of a type-n solute at position (x i , y j ) relative to an edge dislocation aligned along z is denoted as U n (x i , y j ). The theory then considers a generic wavy dislocation characterized by unknown length scales ζ and w. The fluctuations in dislocation energy caused by the interactions with the random solutes are characterized by the standard deviation of the energy of a unit Burgers vector b as the dislocation glides a distance w along the x-direction, ΔẼ p ðwÞ ¼ X i;j;n c n U n ðx i À w; y j Þ À U n ðx i ; y j Þ 2 " # 1 2 : (1) From this solute/dislocation energy per unit length and the elastic energy cost to create the waviness due to the dislocation line tension Γ, the characteristic wavelength scale ζ c (w) as a function of the amplitude w is obtained by minimizing the total energy with respect to ζ. The required stochastic analysis is discussed in ref. 8 with the final result Here, the coefficient 1.73 is derived, not fitted, in the analysis 8 . The characteristic amplitude parameter w c is then computed by minimizing the total energy versus w, which reduces to the solution of The characteristic energy barrier facing a segment of length ζ c that lies in a local favorable energy minimum is then derived to be where again the coefficient 1.11 is derived. At an applied resolved shear stress τ, this characteristic barrier is overcome by a combination of thermal activation and the work −τbζ c x done on the length ζ c segment as it glides a distance x relative to the minimum energy position. For a sinusoidal energy landscape, the stress-dependent energy barrier is well approximated within 1% error over the whole stress regime by where τ y0 is the zero-temperature flow stress, given as The above results are all analytical, requiring only the solute/ dislocation interaction energies, the dislocation line tension, and the Burgers vector. The line tension is taken to be Γ ¼ 1 12 μ f110gh111i b 2 where μ {110}〈111〉 is the alloy shear modulus for shear on the {110} glide plane in the 〈111〉 direction of glide.
In this paper, the theory is quantitatively assessed against carefully-executed direct numerical simulations of the barriers encountered by dislocation motion through true random alloys. However, one cannot simply insert a dislocation of arbitrary line length into a simulation cell and follow its motion under stress. The characteristic waviness associated with ζ c must be captured, which requires the use of a very long line length L >> ζ c . This is computationally costly. A second strategy is to use a line length of precisely L = ζ c . According to the theory, such a segment should remain essentially straight because the line tension prevents the dislocation from becoming wavy on smaller scales. A dislocation segment of length ζ c is then predicted to encounter an energy landscape along the glide direction that consists of local minima and maxima at an average spacing of w c and with an average energy difference (barrier) of ΔE b . We adopt the second strategy below and use the Nudged Elastic Band method to explicitly compute the energy landscape and configurations of a dislocation of length ζ c in the random alloy (see "Methods"). This study thus reveals the statistical distribution and average value for both w c and ΔE b , from which the distribution of barriers and the strength τ y0 can then be determined.

Material properties
We study two model BCC random alloys, NbTaV and MoNbTaW. In these alloys, we compute the solute/dislocation interaction energy U n (x i , y j ) of a type-n solute using the Zhou et al. EAM-type interatomic potentials 10,11 . The Zhou et al. potentials have been reasonably validated for this problem 8 . More importantly, their accuracy for real alloys is irrelevant for the present paper-we will compare simulation and theory on these well-defined model alloys. For the given alloy composition, we first construct the homogenized reference average-atom (A-atom) interatomic potential from the underlying atom-specific potentials 12 . The A-atom potential accurately and automatically encodes, with no fitting, all the average properties (elastic constants, lattice constant, stacking fault energies, etc.) of the true random alloy. Any individual atom type n can be substituted into any position in any configuration of the A-atom material to accurately obtain the energy of the n atom at that position averaged over all possible configurations of atoms in all surrounding positions. The interaction between the n and A atoms thus encodes all the average solute-solute interactions in the alloy. Using this feature of the A-atom potential, an edge dislocation of Burgers vector a 〈111〉/2 is created in the homogeneous A-atom material using standard methods in a large simulation cell. A solute of type n then replaces an A-atom at position (x i , y j ) near the dislocation, and the energy of the system is measured. The interaction energy U n (x i , y j ) is this energy minus the energy of the type n solute in the dislocation-free crystalline A-atom BCC material. Use of the A-atom neglects explicit solute-solute interactions in the alloy but recent theory 13 that includes such solute-solute interactions shows this effect to be small especially in strong alloys. Moreover, explicit results for MoNbTaW using both Zhou EAM potentials and first-principles methods show that solute-solute interactions have a negligible effect on strength in this alloy 13 .
The solute/dislocation interaction energies for two solutes in each of the two alloys studied are shown in Fig. 1, as examples. The individual energies are not large, typically 0.15 eV or less, and are close to the elasticity estimate U n (x i , y j ) = − p(x i , y j )ΔV n where p (x i , y j ) is the dislocation pressure field at (x i , y j ) and ΔV n is the misfit volume of a type-n solute in the alloy. Although the individual interaction energies are small, the collective random fluctuations of the solutes over the scale (ζ c , w c ) creates large barriers for the dislocation motion as shown below.

Results from NEB simulations
From the interaction energies for all solutes in each alloy, the theory predicts the characteristic length scales (ζ c , w c ), the energy barrier ΔE b , and the zero-temperature yield stress τ y0 as shown in Table 1. Direct simulations of long edge dislocations inserted into the random alloy and relaxed at T = 0 K were previously performed by Maresca and Curtin 8 and the length scales (ζ c , w c ) were then derived from examination of the correlation function along the long wavy dislocation. These values are shown in Table  1 and agree well with the theoretical values. Here we compare the direct NEB results to the theory predictions.
We first present results for the NbTaV random alloy (see "Methods"). Figure 2 shows the entire energy landscape. Many local minima and maxima are found, as expected, indicated by blue and red dots, respectively. From the two random realizations of total length 2500 Å, N = 71 maxima are found. The average barrier ΔE b is found to be 1.95 eV, matching the theoretical value very well (see Table 1). The relatively large standard deviation for the barriers is expected. The average barrier spacing that corresponds to w c is w NEB including these decreases ΔE b and w c by 7% and 6.5%, respectively, thus having little effect on τ y0 . As seen in Table 1, the value of w c is somewhat larger than both the predicted value and the value deduced from the correlations along an individual relaxed dislocation. We have no explanation for the differences in w c , which affects the NEB-derived strength (see below), but have performed further convergence checks in the NEB to verify that no significant barriers were missed in the calculations. We do note that the standard deviation of w c is rather large and encompasses the theoretical value.
The theory predicts that the dislocation line at length L = ζ c should remain essentially straight. Figure 3 shows the atomistic configurations of 12 successive minimum-energy dislocations along the glide plane using Common Neighbor Analysis in OVITO for visualization 14,15 . The cores are spread on the glide plane, as expected (see Fig. 1), and variations from a straight dislocation are quite minimal. As noted by Varvenne et al. 16 , the theory based on line tension does not apply for fluctuations at the scale of the atomic spacing. So roughness on the atomic scale can exist but is outside the scope of the theory. It was argued in Varvenne et al.
that such roughness leads to small energy barriers that are easily overcome at finite temperature, so that strength at finite T is mainly controlled by the larger-scale barrier over glide distance w c .
The theory also assumes that the dislocation remains essentially straight as it moves through the energy landscape from minimum to maximum and then back to the next local minimum. Figure 4a, c shows typical examples of the dislocation configurations corresponding to different average dislocation position along the minimum energy path between two adjacent minima, for two different barriers. The dislocation again remains fairly straight as it moves over the distance 2w c , in agreement with the theory.
Similar results are found for the MoNbTaW alloy. Each individual NEB simulation provides a distance between adjacent minima and the path between them, with an associated energy barrier. Over 100 NEB simulations, the mean of the local minima/maxima spacing w c and mean of the energy barrier ΔE b are both shown in Table 1. Very good agreement is found between our NEB simulations, the parameter-free theory, and direct atomistic simulations of long wavy dislocations. A typical example of the NEB path and atomistic configurations along the path are shown in Fig. 5 for MoNbTaW. As found for NbTaV, the dislocation segment remains quite straight as it moves through the barrier from initial to final configuration (see Fig. 5). Similar behavior is found for the other individual configurations.
The above results for both NbTaV and MoNbTaW demonstrate both qualitatively and quantitatively the features emerging from predictions of the theory. For a dislocation segment length fixed at the theoretical value of ζ c , the NEB studies find that (i) the average energy barriers are very large (≈2 eV) in very good agreement with the predicted barriers, (ii) the dislocation energy landscape corresponds to minima and maxima separated by a typical spacing w c that is comparable to the theory, and (iii) the dislocation segments of length ζ c remain quite straight throughout their entire migration through the rough energy landscape. The essentially 1-dimensional energy landscape postulated and   predicted by the theory is thus fully confirmed by the present simulations.
We now turn to the final fundamental aspect of edge dislocation motion: the evolution of the energy barrier under stress and the associated alloy flow stress as a function of temperature and strain rate. It is computationally demanding to perform direct NEB simulations over a range of stresses for each of the hundreds of barriers studied here and so we take an accurate analytic approach as follows. Starting with the zero-stress energy landscape E(x) as a function of mean dislocation position x (see Figs. 2, 4, and 5), the work done by an applied resolved shear stress τ on the glide plane is −τbζ c x. Even if the dislocation does not remain perfectly straight, the area swept during the motion is still exactly bζ c x when x is the mean dislocation position and the work done remains the same. The total energy landscape at stress τ is then E tot (x) = E(x) − τbζ c x. In this new landscape, the individual energy barriers between local minima and adjacent maxima are computed numerically from the NEB landscape E(x). With increasing applied stress, each barrier is reduced steadily, and we can compute the spectrum of barriers ΔE b (τ) at any τ. For each individual barrier, there is a zero-temperature "flow" stress at which the stable local minimum disappears and the barrier becomes zero; this occurs when dE tot (x)/dx = 0 locally. To explicitly validate the analytic approach, we performed direct NEB simulations versus stress for a few typical barriers in MoNbTaW and find that the barrier height versus stress matches the analytical value nearly perfectly in all cases. Figure 6a, b shows the cumulative distribution of energy barriers over a range of applied shear stresses for both the NbTaV and MoNbTaV alloys. At zero stress, the barriers span a wide range from near 0 to 6 eV. With increasing stress, all barriers are steadily reduced, leading to a leftward shift. Furthermore, an increasing number of barriers are reduced to zero; these are included in the cumulative distribution and so are reflected in the non-zero initial value of the cumulative probability at zero energy.
The theory is developed in terms of a characteristic (mean) barrier and hence mean strength, so that we should compare the theory predictions to the stress at which 50% of the barriers have been reduced to zero. The original formulation of the theory rationalized this choice is as follows. A long dislocation moves through a rough (stochastic) potential energy landscape with characteristic length scales ζ c along the line direction and w c along the glide direction, and characteristic energy barriers between adjacent minima of ΔE b . Individual segments of length ζ c encounter the spectrum of barriers as the entire dislocation moves forward. However, individual segments cannot advance much beyond nor lag much behind the average dislocation position. In advancing beyond, segments are pulled back by the remaining line. In lagging behind, segments are pulled forward by the advanced line. The motion of the long dislocation is thus controlled by neither statistically weaker nor stronger local barriers. This behavior can be observed in direct simulations, but is difficult to quantify. At 50% probability, there are an equal number of advanced and lagging segments, i.e., the dislocation retains the roughness of the relaxed dislocation in the zero-stress landscape, and there are no net forces due to segment interactions. Since our goal is to compare to the theory, we use this 50% level here.
The characteristic 50% stress for the zero-temperature strength as obtained from Fig. 6a, b is shown in Table 1 for each alloy. The mean strength for NbTaV is 293 MPa. This value only~60% of the theoretical prediction, a difference that stems solely from the larger typical value of w c found in the NEB simulations while the average energy barrier ΔE b is very similar to the theory. The w c deduced from the wavy dislocation scale in NbTaV was, in contrast, slightly smaller than the theory prediction. As discussed earlier, the origins of these differences in w c are unclear. The NbTaV alloy does show a comparatively large fraction (~15%) of rather low barriers (<0.4 eV) that are eliminated at very low stress; this may contribute to the reduced average strength. The mean strength for MoNbTaW is 410 MPa, which compares quite well to Finally, in solute strengthening theories, the energy landscape is assumed to have a locally sinusoidal shape between the minimum and the maximum, which leads to the stress-dependent energy barrier given by Eq. (5). Here, we assess this assumption by extracting the values ΔE b and w c from the direct NEB results and use the analytic result of Eq. (5) based on a sinusoidal model to predict the energy barriers at finite applied stress. The results are shown in Fig. 6, and the distribution agrees very well with the full analysis based on the actual (non-sinusoidal) energy landscapes found in the NEB at all stress levels. The sinusoidal approximation of the energy barrier, and its associated T = 0 K strength, is thus quite accurate.
Application to real alloys Both the analytical theory and direct NEB show that there are very large zero-stress energy barriers of 1.95 eV for NbTaV and 2.45-2.65 eV for MoNbTaW. These values are significant, meaning that high stresses are needed to overcome these barriers, even at high temperatures, and thereby cause plastic flow in the alloys. The strain-rate and temperature-dependent flow stress of an alloy follows from the thermally-activated Arrhenius model that connects the plastic strain-rate _ ϵ to the energy barrier via _ ϵ ¼ _ ϵ 0 expðÀΔEðτÞÞ=kT where _ ϵ 0 ¼ 10 4 s −1 is a reference strain rate (see refs. 8,17 ). Combining this relationship with Eq. (5) and inverting gives the finite-temperature, finite strain-rate flow stress τ y ðT; _ ϵÞ. This equation holds for low temperatures and/or high stresses (τ y /τ y0 > 0.5) 8 . At higher temperatures and/or low stresses, the long dislocations (L << ζ c ) can explore higher wavelength fluctuations with larger barriers that have a similar scaling 18 . To cover the full range of temperatures, Maresca et al. 8 proposed the approximate form τ y ðT; _ ϵÞ ¼ τ y0 exp À 1 0:55 The uniaxial tensile yield strength in a polycrystalline alloy is then computed as σ y = Mτ y , where M = 3.067 is the Taylor factor for an untextured BCC polycrystal deforming by edge dislocation motion. It is clear that the large values of ΔE b lead to significant strength retention with increasing temperature. To enable application of the edge theory for alloy design, Varvenne et al. 16 developed a simplified elasticity theory for FCC alloys that was later extended to BCC alloys 9 . The theory uses the elasticity approximation U n (x i , y j ) = − p(x i , y j )ΔV n , where ΔV n is the misfit volume of a type-n atom in the alloy and p(x, y) is the edge dislocation pressure field. Assuming a common dislocation structure (spreading of the BCC edge dislocation) then led to analytic estimates of the T = 0 K strength and energy barrier as  Here, μ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi C 44 ðC 11 À C 12 Þ=2 q , ν ¼ 3BÀ2μ 2ð3BþμÞ are the isotropic elastic constants with Bulk modulus B ¼ ðC 11 þ 2C 12 Þ, all computed using a simple rule-of-mixtures of the elemental values. b is the alloy Burgers vector computed using Vegard's law for the alloy atomic volume, and α = 1/12 is a numerical coefficient related to the dislocation line tension. The strength versus temperature and strain rate can then be computed using Eq. (7). The numerical coefficients above were determined by fitting to results of the full theory across a range of alloys in the Mo-Nb-Ta-V-W family, with the fit being most accurate for the high-barrier/high-strength alloys. The simplified theory depends only on elastic moduli and atomic misfit volumes, and was implemented for preliminary design across 10 7 different compositions in the Al-Cr-Hf-Mo-Nb-Ta-Ti-W-Zr family 9 to identify promising high-T alloys.
Here, we use the reduced analytic theory to understand the strengths in several alloys recently reported. The first group of alloys consists of the base alloy of HfNbTaTiZr with additions of Mo. Specifically, Mo was substituted for Nb and Ta to create HfTaTiZrMo and HfNbTiZrMo and was also added to the base alloy to create HfNbTaTiZrMo; note that our labeling of the alloys now does not follow conventional alphabetical order of elements but instead appends the new element (Mo) to the end of the elemental list. Table 2 shows the atomic volumes and elastic constants used to make predictions via Eq. (8). In this family of alloys, Mo has the smallest atomic volume, leading to the largest misfit volume, and also the largest elastic moduli. Mo substitution or addition to HfNbTaTiZr is thus predicted to give significant strengthening. Figure 7a shows the experimental strengths as reported in refs. 3,19 versus the predicted strengths, at temperatures T = 0, 1000, and 1200 ∘ C, and strain rate _ ε ¼ 10 À3 s −1 . Agreement between theory and experiment is very good at both T = 0 and 1000 ∘ C. Most notably, the theory captures the very significant strengthening achieved due to the addition or substitution of Mo. We discuss the deviation at 1200 ∘ C below.
The second group of alloys consists of the base alloys MoNbTaW and MoNbTaTiW to which V is added to create MoNbTaWV and MoNbTaTiWV, respectively, again with the added element V appended to the end of the alloy label. Table 2 shows the atomic volumes and elastic constants for these elements, see ref. 9 . V has the smallest misfit volume among all these elements while its elastic constants are roughly comparable to those of Nb, Ta, and Ti; the V-containing alloys should thus be stronger than the base alloys but with less of an effect than Mo in HfNbTaTiZr. Figure 7b shows the experimental strengths versus the predicted strengths, again at T = 0, 1000, and 1200 ∘ C and strain rate _ ε ¼ 10 À3 s −1 . While the quantitative agreement is not as good, the theory captures the notable strengthening due to the addition of V in both alloys. We note that the enhanced strength of MoNbTaTiW relative to MoNbTaW, i.e., the role of the Ti addition to MoNbTaW, is not captured by the present theory; one possible origin is O or N contamination, well-established to significantly increase alloy strength.
The Hf-Mo-Nb-Ta-Ti-Zr alloys in Figure 7a show a strong loss of strength between 1000 and 1200 ∘ C while the Mo-Nb-Ta-Ti-V-W alloys do not show such a loss of strength and instead continue to agree with the theory trend. High-temperature mechanisms that can defeat edge strengthening are unclear at this time. The random alloy presents edge barriers throughout the material on the scale of w c and ζ c and so avoiding these barriers seems difficult. The solute/dislocation interaction energies scale with alloy elastic moduli, but their decreases with temperature are modest for the refractory elements Mo, Nb, Ta, V, and W 20,21 . We thus postulate that there is some vacancy-related mechanism, such as climb, by which the large edge barriers can be circumvented at sufficiently high temperatures. Such a The list shows the elemental BCC atomic volume, the elastic moduli, and the experimental vacancy formation energy E vf . mechanism would then scale with the vacancy formation energies in these alloys families. In the four Hf-Mo-Nb-Ta-Ti-Zr alloys, the vacancy formation energies estimated by a simple rule-ofmixtures (see Table 2) are in the range of 2.2-2.3 eV. In contrast, the estimated vacancy formation energies for the Mo-Nb-Ta-Ti-V-W alloys are higher, in the range of 2.65-3.1 eV. Thus, a mechanism operating at 1200 ∘ C in the Hf-Mo-Nb-Ta-Ti-Zr alloys that is related to the vacancy concentration ( c V ¼ e ÀE vf =kT % 2 10 À8 ) would require a temperature of 1600 ∘ C in the Mo-Nb-Ta-Ti-V-W alloys to achieve the same vacancy concentration. The MoNbTaW and MoNbTaWV alloys do not show a significant drop in strength up until this range of temperature 22 . This analysis can also be related to the alloy melting points. The Hf-Mo-Nb-Ta-Ti-Zr alloy melting points have been estimated to be in the range of 2270-2360 K 3 while the melting points of the Mo-Nb-Ta-Ti-V-W alloys are much larger, in the range of 2620-3110 K.
Hence the latter alloys are expected to retain strength up to higher temperatures 3 . In both families of alloys, the ratio of estimated vacancy formation energy to melting temperature is similar, E vf /k b T m ≈ 11.5+/0.5, so that we cannot distinguish between vacancy-related mechanisms and any other high-T mechanisms that might operate. It remains possible that the strength loss at high T arises because the strength is controlled by screw dislocations. Screw dominance has been observed in some BCC HEAs 23 , and is the accepted situation for dilute BCC alloys. Theoretically, Maresca et al. 8,24 have shown that some BCC non-dilute alloys and HEAs are controlled by screw dislocations and others by edge dislocations. The high-temperature strength of screw dislocations is controlled by cross-kinks that can be defeated by thermal vacancies 24,25 . Hence, loss of screw strength would also be correlated with E vf and/or T m . While this mechanism is clear, it is difficult to determine reliable alloy parameters for application of the screw theory. It is thus not possible to determine whether the quantitative strength increase upon addition of Mo to HfNbTaTiZr could be attributed to strengthening of screw dislocations. In contrast, the strength at 0 and 1000 ∘ C, and the role of Mo, are well-captured by the parameter-free edge theory. The edge theory should thus be a lower bound for strength and, being analytic, serves as a very useful theory for alloy design/selection. The mechanism(s) of strength loss at high T remain to be resolved; this is a critical open question for BCC HEAs.
In summary, the recent theory of strengthening of edge dislocations in BCC high entropy alloys has been quantitatively supported by detailed atomistic NEB simulations of the complex energy landscape for edge dislocations in two random alloys, NbTaV and MoNbTaW. The analysis reveals the high energy barriers encountered by edge dislocations, and the magnitude and spacing of these barriers are in good agreement with the analytic predictions of the theory. The simulations here further show that the dislocations remain quite straight over the characteristic scale ζ c , as predicted by the theory. The stresses required for dislocation motion are also comparable to those observed in limited direct simulations. The present work thus supports the theory for critical features of strengthening in two different refractory BCC HEAs. A reduced version of the theory has then been applied to two families of alloys. The theory rationalizes the observed significant increases in strength caused by the addition of solutes with large misfit volumes that is a hallmark of the edge theory. The high-temperature mechanisms of strength loss remain unknown but appear correlated with vacancy formation energy and/or melting temperature. The reduced analytic theory (Eq. (8)) provides a validated mechanistic path for selection of BCC HEA compositions for high strength that can be reaffirmed by application of the full theory for the most promising compositions. The strength theory can also be combined with other performance metrics to identify new compositions that are strong at high temperatures, ductile, single-phase, and/or oxidation-resistant, fulfilling the multiperformance requirements of many applications.

Molecular simulations of random alloys
For both model alloys, a rectangular simulation cell oriented with glide direction X||[111], glide plane normal Y||[101], and line direction Z||[121] is first created using the lattice parameter of the corresponding "average atom" alloy. The dimensions of the simulation cell are (L x = 1327, L y = 123) Å and (L x~1 90, L y~9 0) Å for NbTaV and MoNbTaW alloys, respectively. The huge difference between the sample sizes L x for the two alloys stems from different methodologies used to find the edge dislocation minimum energy positions as discussed later. As explained above, the dislocation line length is fixed at the characteristic length of L z = ζ c~2 9 Å for NbTaV and L z = ζ c~5 0 Å for MoNbTaW alloys (see Table 1). The surfaces normal to the Y-direction are traction-free while those along the X-and Z-directions are periodic.
To insert a single edge dislocation, the simulation cell created above is first populated with the corresponding average atoms. Two adjacent YZatomic planes in the lower half of the simulation cell are then removed. The average atom at each lattice site is then randomly replaced with a solute atom of the alloy according to the solute concentration in the alloy. The created simulation cell (with two missing atomic planes in the lower half) enables us to introduce the edge dislocation at any desired position x*. First, the atomic positions in the lower half of the simulation cell are shifted in the x-direction (keeping the same randomness) such that the gap due to the missing atoms is located at x*. The empty space in the lower half at location x* is then eliminated by displacing the atoms on both sides by u x = b (2x ± L x )/(2x* ± L x ) where ± reflects positions to the left and right of x*. The true random alloy with the dislocation is finally created by relaxing the atomic positions. In order to minimize computation time, the Conjugate Gradient algorithm and then the FIRE algorithm are used in sequence, and the domain is relaxed until convergence is achieved (force tolerance < 10 −5 eV per atom). All simulations are performed using molecular statics as implemented in LAMMPS 26 .
For the long NbTaV random alloy sample, the minimum energy positions are first found by placing an initial dislocation at every possible location (spacing of one Burgers vector) along the slip plane and then minimizing the energy. During minimization, the dislocation moves from its initial position to the local minimum energy position within the local energy basin 27 . Many close initial positions thus have the same final position. NEB simulations are then performed to find the path and barrier between each pair of adjacent minima, using a sub-domain of L x = 200 Å centered around the pair of adjacent minima. Two large random realizations were studied, with the minimum energy paths (MEP) obtained using 32 and 50 replicas, respectively, between each pair of initial and final configurations. No significant statistical differences are found between the two realizations, and their results are aggregated. The inter-replica spring constant is set to 1 eV Å −1 (force units) with "ideal-parallel" option in LAMMPS which calculates nudging forces as described in ref. 28 and results in more equally-spaced replicas. Using the same spring constant used for the MoNbTaW NEB computations (see below) results in very similar barriers. Each stage of the NEB calculation (convergence to MEP and barrier climbing) is executed for 10,000 steps and convergence is checked afterwards by observing the MEP progression. These NEB computations were performed for every pair of adjacent minima along the entire length L x~1 300 Å.
For the MoNbTaW alloy, after inserting an edge dislocation in the system and relaxing to the local minimum configuration, a dislocation is again introduced into the same initial simulation cell (same random configuration) at a distance 2w c from the previously-found minimum position. This is the distance around which the adjacent minimum is expected. Full relaxation then enables the dislocation to reach the actual minimum energy position. This procedure establishes the initial and final dislocation positions for an NEB computation of the path and energy barrier between them. This procedure was repeated for over 100 MoNbTaW random realizations of the alloy to generate a statistical distribution of barriers. In no case was another minimum found in between the initial and final positions of any realization. NEB simulations were performed using 100 replicas and inter-replica spring constant is set to 10 −2 eV Å −2 . Convergence was assumed when the maximum force acting on all of the atoms across all replicas was <1 × 10 −3 eV Å −1 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The open-source software LAMMPS was used for molecular dynamics simulations.