Ab initio modeling of the energy landscape for screw dislocations in body-centered cubic high-entropy alloys

In traditional body-centered cubic (bcc) metals, the core properties of screw dislocations play a critical role in plastic deformation at low temperatures. Recently, much attention has been focused on refractory high-entropy alloys (RHEAs), which also possess bcc crystal structures. However, unlike face-centered cubic high-entropy alloys (HEAs), there have been far fewer investigations into bcc HEAs, specifically on the possible effects of chemical short-range order (SRO) in these multiple principal element alloys on dislocation mobility. Here, using density functional theory, we investigate the distribution of dislocation core properties in MoNbTaW RHEAs alloys, and how they are influenced by SRO. The average values of the core energies in the RHEA are found to be larger than those in the corresponding pure constituent bcc metals, and are relatively insensitive to the degree of SRO. However, the presence of SRO is shown to have a large effect on narrowing the distribution of dislocation core energies and decreasing the spatial heterogeneity of dislocation core energies in the RHEA. It is argued that the consequences of the mechanical behavior of HEAs is a change in the energy landscape of the dislocations, which would likely heterogeneously inhibit their motion.


INTRODUCTION
Previous investigation of the fundamentals of deformation in body-centered cubic (bcc) transition metals have revealed that the core properties of the ½〈111〉 screw dislocations play an essential role in their plasticity 1 , especially at low temperatures where the deformation is thermally activated through the kink-pair nucleation mechanism 2 , and is expected to be strongly temperature dependent. The high lattice friction associated with such screw dislocation motion is a result of nonplanar core structure 1,3 and is related to the height of the Peierls potential 4 .
Due to the importance for plastic deformation, extensive atomistic simulation studies have been devoted to computing core structures and corresponding mobilities of screw dislocations in bcc transition metals 3,[5][6][7][8] . In these studies, one of the significant challenges has been the variation in properties derived from different models for the interatomic potentials. For example, early studies based on classical potential models often predicted a metastable split core structure [9][10][11] , which leads to a camel-hump shape in the Peierls potential. Later density functional theory (DFT) calculations produced symmetric and compact dislocation cores in Mo, Ta, and Fe [12][13][14][15][16] ; similar compact cores have been found in other bcc transition metals, such as W, Nb, and V 17,18 . In DFT studies of the energy landscape of screw dislocations in bcc transition metals [18][19][20] , it was found that nondegenerate cores lead to a single humped curve in the Peierls potential, implying that the split core structure might not be metastable. Alloying effects on the Peierls potential of W have also been explored 21 . Recently developed machine learning-based potentials [22][23][24] and new embedded atom method potentials that consider quantum effects on lattice vibrations 25 and extra constraints 26 all lead to predictions of a single humped curve in the Peierls potential. Due to the dependence of the results for screw dislocations in bcc transition metals on the model for interatomic bonding, DFTbased approaches are of interest to provide benchmarks for subsequent modeling at higher scales.
During the past 15 years, a new class of alloys known as highentropy alloys (HEAs) 27,28 has drawn extensive research interest. These alloys involve multiple principal elements (typically five) in nominally equimolar ratios, and were originally presumed to crystallize as a single-phase solid solution. As a new class of structural materials, some types of HEAs, in particular the CrCoNibased alloys, have been shown to possess exceptional damage tolerance and improved strength at cryogenic temperatures 29,30 . Theoretically, mechanistic, first-principles-based predictive theories for the temperature, composition, and strain rate dependence of the plastic yield strength have been developed and applied to such face-centered cubic (fcc) alloys [31][32][33] . Indeed, most HEA research to date has been focused on these fcc "Cantor-type" alloys 34,35 , whereas a second distinct family of HEAs, comprising mostly refractory elements, has been far less studied. Such refractory high-entropy alloys (RHEAs), which are sometimes termed Senkov alloys 36,37 , invariably crystallize in bcc solid solution phases that have been designed for elevated temperature applications 38 .
For example, RHEAs such as MoNbTaW with single-phase bcc crystal structures have been produced by vacuum arc melting 37 or direct metal deposition 39 with exceptional microhardness 36 , as well as excellent compression yield strength and good ductility at high temperatures 37 . Transmission electron microscopy (TEM) studies on RHEAs have shown a dominant role of screw dislocations with increasing plastic strain 40,41 , similar to traditional bcc metals. Additionally, strong intrinsic lattice resistance has been found in certain RHEAs 41,42 . To model such behavior, molecular dynamics (MD) simulations have been used to study dislocation behavior in bcc RHEAs 43 . For example, screw dislocation core structures in NbTiZr, Nb 1.5 TiZr 0.5 , and Nb 0.5 TiZr 1.5 alloys were recently explored using MD simulations, and significant core structure variation was found along the dislocation line 44 . Recent theory has revealed the potential importance of edge dislocations in controlling the strength of bcc HEAs at high temperatures 45 and the correlation between atomic distortions and the yield strengths of HEAs 46 . However, there are still only very limited studies on the deformation behavior of this new class of bcc alloys, as compared to single-phase bcc transition metals.
Another important aspect of HEAs is the presence of local chemical short-range order (SRO). Although these alloys can be described as "topologically ordered yet chemically disordered", the local chemical environments are unlikely to be characterized by a perfectly random distribution for every atomic species [47][48][49][50][51] . Indeed, their disordered multiple-element compositions lead to a strong possibility of SRO, for example, the preference for certain types of bonds within the first few neighbor shells. This is not particularly rare in conventional alloys 52,53 and glasses 54 ; however, it could be argued that its existence would be even more likely in multiple principal element alloys 49,51,55 due to large number of elements and their equimolar concentrations. Recent DFT and MD simulations on the fcc CrCoNi alloy suggest that SRO can have a profound effect on critical parameters, notably the stacking-fault energy 55 and dislocation mobility 56 ; accordingly, such local order could be an important factor in controlling mechanical properties.
In spite of extensive studies on the bcc transition metals, there are relatively few published studies of dislocation core structures, dislocation mobility, or the effect of chemical SRO for bcc RHEAs. Accordingly, the objective of the current paper is to employ DFTbased methods to compute the dislocation core structures in refractory HEAs and to explore the distribution of dislocation core energetics and its potential effect on Peierls barriers, focusing on the MoNbTaW system.

Dislocation core structures in RHEAs
To compute the core structures and Peierls potential for ½〈111〉 screw dislocations in the refractory MoNbTaW HEA, we employ DFT calculations, making use of the Vienna ab initio simulation package (VASP) [57][58][59] ; details of the DFT calculations are provided in the "Methods" section. For screw dislocations in refractory HEAs, we employ a periodic supercell that contains 462 atoms, as illustrated in Fig. 1a. The simulation cell contains a pair of dislocations with opposite Burgers vectors, in a nearly square quadrupolar arrangement 16 with triclinic symmetry to minimize any effects of periodic boundary conditions and image stress. This dipole approach was first introduced by Bigger et al. 60 and has been widely used in DFT calculations on dislocations 16,17,20,61 . The supercell adopted in current work was previously described by Weinberger et al. 17 and Li et al. 61 , and was used to calculate dislocation core structures in pure bcc transition metals. In addition, since the size of supercell is fixed for all the simulations, the short periodic length might have some influence on the dislocation dipole energy due to its effect on the nature of the SRO. In our current model, we consider an equimolar MoNbTaW bcc RHEA 37 , and doubled the periodic length along the dislocation line direction of the original 231-atom model (see "Methods" section for further details) to minimize as much as possible correlations in the chemical order, as described in the following section.
The initial atomic configuration was generated by creating a special quasi-random structure (SQS) on the 462-atom supercell shown in Fig. 1a. The SQS was generated using the Alloy Theoretic Automated Toolkit (ATAT) program 62 . The SQS methodology was used to minimize chemical correlations, and thus to provide a reference configuration corresponding to random substitutional disorder (i.e., minimizing chemical SRO). This reference configuration was used in Monte-Carlo (MC) simulations to generate supercells with varying degrees of SRO, as described below. For each of the configurations with different level of SRO, we shifted the dislocation dipole over all the possible sites within the simulation cell, to statistically sample dislocation properties. The atomic positions in the system with the dislocation dipole were then relaxed to enable interrogation of the core structures and energies in different lattice sites within the RHEA supercell.
For each configuration representing a different degree of chemical SRO, we calculated 231 different supercells with the dislocation dipole, in which the position of the cores initialized in different local environments. We found that the screw dislocations in bcc MoNbTaW HEAs maintain a compact core structure in most of the resulting relaxed structures, as illustrated by Fig. 1b, which is similar to the case in pure bcc elements 17,18 . In very few a b c [111] [112] [110] situations, the core can be extended on the (110) plane as shown in Fig. 1c 63,65 , as discussed below), chosen to enable the development of dislocation supercell models with representative degrees of chemical SRO. Our focus is specifically on the effect of SRO on the dislocation properties. For generating supercells with different degrees of SRO, similar to previous studies in fcc HEAs 47,55 , we applied a DFT-based lattice MC approach to our 462-atom supercell model; details are described in the "Methods" section. The supercell initiated with an SQS configuration was used as input for the MC simulations. The MC simulation samples swaps of atom types, following the Metropolis algorithm, and the entire simulation considers~2100 such swaps, leading to the evolution of the energy shown in Fig. 2a. Due to the limited number of MC steps and the lack of sampling of atomic displacements, the final configurations may differ from the true equilibrium state of SRO at the simulation temperature, although they appear to be quite close to the state of SRO as calculated by Kostiuchenko et al. 65 for high temperatures (~1200 K). However, the algorithm does lead to appreciable lowering of the energy, as shown in Fig. 2a, and the pair-forming tendencies shown in Fig. 2b are consistent with previous work on SRO in the same system using more comprehensive methods 63,65 , as discussed below. Thus, this method is used to generate representative samples with varying degrees of chemical SRO to explore the resulting effect on dislocation properties.
Similar to the conventional Warren-Cowley description 66 and the previous study for fcc HEAs 55 , we characterize the state of SRO using the so-called nonproportional number of local atomic pairs, Δδ ij , as described in more detail in the "Methods" section. Based on our calculations, the evolution of total potential energy and the overall chemical SRO ð P i;j Δδ ij Þ in the sample during the MC relaxations are plotted in Fig. 2a. With respect to axes, the abscissa is the total potential energy change of the system and the ordinate is the overall chemical SRO of the system. As the MC simulation proceeds, the potential energy of system decreases monotonically, while the chemical SRO increases at the same time. This clear trend indicates that chemical SRO is occurring in the system with the MC simulations. To quantify the effect of SRO on dislocations, three different samples from the simulation (s1, s2, s3) were chosen for further calculation of core structures and energies, indicated by the red arrows in Fig. 2a. State s1 represents the nearly random solid solution configuration with lowest magnitudes of the SRO parameters; s2 represents an intermediate configuration with a medium level of SRO, and s3 represents the configuration with the highest degree of SRO. Figure 2b shows the quantitative values of Δδ ij between all the species in the MoNbTaW alloy; the red dots show that the local SRO in state s3 clearly deviates from the random solid solution. Preferred atomic pairings between Mo-Ta, Mo-Nb, and Ta-W were observed as the Δδ ij values are 0.308, 0.196, and 0.112, while unfavorable pairings between Mo-W and Ta-Nb were also apparent as the Δδ ij values are −0.392 and −0.294. This result confirms the energetic preference for SRO in MoNbTaW alloys; moreover, the tendency to form SRO that we see here is consistent with previous studies using other methods 50,51,63,65 that have shown the Mo-Ta pairs are the most dominant contributors to the SRO, followed by Ta-W and Mo-Nb pairs.

Distribution of dislocation core energies in bcc RHEAs
After the introduction of SRO through MC relaxations, the dislocation dipole described in Fig. 1a was created in samples s1, s2, and s3. To sample over the distribution of local chemical environments for dislocation cores, the dislocation dipole was shifted over all the possible sites within the simulation cell leading to 231 different configurations for each of the three states of SRO. All the configurations with the dislocation dipole were then minimized, following the procedures described in the "Methods" section. Figure 3 shows histograms of the supercell excess energies, that is, the energy difference between the supercell with and without the dislocation dipole, of all the configurations minimized at different SRO states. The histograms for the three SRO states are fit well by normal distributions (the fitted lines are also shown in Fig. 3). The green dash-dot line represents the energy distribution of the nearly random solid solution sample s1. The blue dash line represents the sample s2 with a medium degree of SRO and the red solid line represents the sample s3 with highest degree of SRO.
The mean values of the excess energies for the two samples with SRO differ by 0.38 eV (s2) and 0.87 eV (s3) from that for the most disordered sample (s1).
To compute dislocation core energies from these energies, we consider the components contributing to the supercell excess energy. The excess energy is the sum of the two dislocation core energies, the elastic energy arising from the dislocations, and a contribution from the diffuse antiphase boundary (DAPB) energy between the two cores created by the relative shift of the crystal by a Burger's vector across the planar "cut" region between the dislocations. This excess energy can thus be written as: We note that in previous studies, it has been shown the excess energies of the types of supercells used in this study can also be affected by the residual stress in the simulation box [67][68][69] . In Supplementary Fig. 2, we plot the distribution of this residual stress on all the simulation cells, and the results rule out the correlation between the change in variance in the excess energies with these residual stresses. For what follows, we thus focus on the decomposition of the excess energies into core, elastic, and DAPB contributions.
To first order, the elastic energy can be estimated using the continuum theory as described by Clouet 67,68 and Clouet et al. 69 , using the elastic constants and dislocation Burgers vector and a reasonable assumption for the core radius. For the simulation supercell used here, the DFT-calculated elastic contribution E elastic is estimated to be~6.0 eV. Importantly, for the analysis that follows, we found that the SRO has only an~3% effect on the calculated elastic moduli (see further details in Supplementary Note 2), such that this local order is estimated to contribute only a 3% change (~0.17 eV between s1 and s3) in the elastic energy contribution to the calculated excess energies. Details of the calculations of the elastic constants and elastic energy contribution in the dislocation dipole cell are provided in Supplementary Tables 1 and 2.
Another contribution to the average and variance in calculated excess energies for the supercells is associated with the cut plane between the two dislocation cores. When SRO is present, this cut plane leads to a contribution to the energy of the supercell arising due to the shift of adjacent planes, which disrupts the state of SRO and causes an excess energy E DAPB . Following the convention in the literature, this planar defect is referred to as a DAPB and can be quantified through the so-called DAPB energy per unit area (γ DAPB ). We have calculated the DAPB energy in our current system (see Supplementary Note 3), with the following results: for the state s1, which represents the random solid solution, γ DAPB is~3 mJ/m 2 , that is, essentially zero within the accuracy of our statistical sampling. With increasing SRO, γ DAPB increases to 29 mJ/m 2 in state s2 and 59 mJ/m 2 in state s3 with the highest degree of SRO. E DAPB associated with the cut plane gives rise to an increasing contribution to the excess energy of the supercell: from 0.015 eV in state s1 to 0.59 eV in state s3. Further, due to the important role of the variance in the core energy distribution, which will be discussed below, the variation in E DAPB due to the position of the cut plane as the locations of the dislocation cores are shifted are also calculated through DFT simulations. The standard deviation σ EDAPB is~0.15 eV for s1, 0.17 eV for s2, and is increased to 0.26 eV for s3. Based on these data, we can further decouple the contribution of the variance in excess energies due to the two dislocation cores and the DAPB. The details of these calculations are shown in Supplementary Note 3 and Supplementary Table 4.
Assuming that the excess energy shown in Fig. 3 can be decomposed as E ¼ 2E core þ E elastic þ E DAPB , we can extract the distribution of core energies by subtraction of the contributions from elastic energy (see Supplementary Note 2) and the mean and variance of the DAPB energy (see Supplementary Note 3). Further details are given in Supplementary Note 4 and Supplementary Table 5. Figure 4 shows the average and variance values for the dislocation core energy in MoNbTaW for different SRO states. The average values are compared with the value in pure bcc transition metals from a previous DFT study 18 . The averaged core energy in MoNbTaW HEA is the highest compared with all its constituent pure elements. In addition, the SRO has only a marginal impact on the averaged dislocation core energy since it is a one-dimensional (1D) line defects; this result is in contrast to the effect of SRO on 2D planar defects, such as the stacking-fault energy 55 or DAPB energies. One important feature of Fig. 3 is that the dislocation dipole energy follows a Gaussian distribution, which is an intrinsic feature of an HEA that differs from the pure element metals. Although the averaged core energy is not sensitive to SRO, the variance of the distribution is found to decrease with the increase of SRO. The standard deviation of the excess energy in Fig. 3 for SRO state s1 is 0.72 eV, which decreases to 0.36 eV in s2 and to 0.37 eV in s3, that is, with lower degrees of SRO, the variance becomes more significant. The variances of dislocation core energies, decoupling the effect of the DAPB energy, are illustrated in Supplementary  Table 4 and show similar trends. As described above, we conclude that the dominant contribution to the variance in supercell energy shown in Fig. 3 arises from the variations in dislocation core energies; the results thus also demonstrate the role of SRO in changing the dislocation core energy distribution.
To illustrate the local spatial variation in core energies in the RHEA, and the effect of SRO on these variations, we plot 3D contours of the supercell excess energies and their 2D projection in the supercell in Fig. 5. For simplicity, each dislocation dipole is treated as a single point located at the average spatial location of the two screw cores in the dipole; they are aligned in ½112 and ½121 directions based on their relative positions. The excess energies normalized by the total length of the dislocation lines, which can be regarded as the depth of the Peierls valleys, are   shifted to set the minimum value equal to zero. Based on these data, the left column of Fig. 5a-c shows 3D contours of the Peierls valleys at different SRO states from s1 to s3 and the right column corresponds to their 2D projection. Note that this is not a minimum energy path contour, since no transition-state data were included in these plots. For a pure element metal, the contour in Fig. 5 would be that of a flat surface since the depth of the Peierls valley has a constant value. However, due to variations in local environment within the RHEA, the dislocation dipole energy in these alloys follows a normal distribution, as shown in Fig. 3, which leads to rugged Peierls valleys contours, as shown in Fig. 5. The maximum variation in the Peierls valleys is 0.9 eV/b in the near-random s1 state; with increasing SRO, this decreases to 0.48 eV/b in s2 and to 0.50 eV/b in s3. It is clearly visible in Fig. 5a-c that the Peierls valley contours contains a rugged feature for the RHEA. Similar to Fig. 3, histograms of the differences in Peierls valley energy for different SRO states are shown in Fig. 6 (see the "Methods" section for further details). As discussed further in the next section, the Peierls valley energy differences considered in Fig. 6 are defined as ΔE = E d1 − E d2 , where E d2 is the excess energy of the supercell for one position of the dislocation dipole, and E d1 is the excess energy when this dipole has shifted by glide to the neighboring Peierls valley in the ½112 direction. If we assume that the dislocation dipole energy follows the same Gaussian distribution shown in Fig. 3, based on the properties of Gaussian distributions, the values of ΔE will also follow a Gaussian distribution but with a different variance: $ Normalð0; 2σ 2 À 2σ cov Þ, where σ cov is the covariance of the excess energy for two neighboring positions of the dipole. In Fig. 6, the average value of the energy difference is zero for all the three SRO states, as expected, and the variance of the fitted distribution from the DFT energy data agrees well with the prediction (details of the calculation of σ cov and ffiffi ffi 2 p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ 2 À σ cov p are given in Supplementary Note 5). Similarly, the values of Peierls valley energy difference (ΔE) defined above, corresponding to glide of the dislocations in the ½112 direction, are plotted in Fig. 7 as a function of the initial position of the dipole, and are represented in both 3D contours and 2D projections. Standard analyses of transitions in complex systems are consistent with the basic trend that the energy difference between the final and initial states correlates with the change in the energy barrier. In the contours of valley energy differences shown in Fig. 7, the values range from −0.30 to 0.30 eV/b in the s1 state; these decrease to −0.12 to 0.15 eV/b in state s2 and to −0.14 to 0.15 eV/b in the state s3. The fraction of these energies with relatively high values decreases with the increasing degree of SRO. These results, along with the change of distribution of dipole energies in Fig. 3, demonstrate that the presence of SRO serves to narrow the distribution of dislocation core energies and decrease the spatial heterogeneity of dislocation core energies in the system. For reference, the Peierls barriers in the pure element constituent metals, Mo, Nb, Ta, and W, calculated through the drag method, which is consistent with DFT study 17 , are also plotted on Fig. 6. A significant amount of the Peierls valley energy difference (ΔE) during glide can be seen to have exceeded the highest value of the Peierls barriers in pure bcc elements. This rugged energy landscape and variance in core energies intrinsic in RHEA is anticipated to have a profound effect on the distribution of Peierls barriers, as explored further below.

Peierls barriers of screw dislocations in bcc RHEAs with local chemical order
In pure bcc transition metals 17,20 , the Peierls barriers for ½〈111〉 screw dislocations can be computed from the energy pathway between two equilibrium samples, in which the dislocation dipole is uniformly translated along the ½112 direction on the {110} plane to the nearest neighboring site using the reaction coordinate method (also termed the "drag method" 70 ) or the "nudged elastic band (NEB) method" 71 . However, in a system with a complex energy surface, such as the RHEA considered here, the NEB method is computationally highly costly and difficult to converge. Alternatively, we have found that the "drag method" converges well.
Based on our tentative estimations of Peierls barriers through the "drag method", the most significant feature in the RHEA system is that the equilibrium energies of the dislocation dipoles are not constant due to the different local environments of the dislocation cores, compared with pure element metals. Thus, the potential energy of the initial configuration (where the reaction coordinate is 0), is generally not equal to that of the final configuration (reaction coordinate of 1). The shape of Peierls potential and the barrier values depend markedly on the relative energy difference between the initial and final configurations and can be divided into two distinct classes that we will refer to as Type-1 and Type-2 barriers, as shown in the schematic plot in the Fig. 8a. Generally, the barrier value is higher than the potential energy difference between the final and initial configurations. When the potential energy difference between the initial and final configurations is small, or when the energy of the final configuration is smaller than that of the initial configuration, the barrier curves are usually Type-1, as shown by the red curve in Fig.  8a. However, if the final configuration has a much higher potential energy than that of the initial configuration, the typical barrier curves under this condition will be like the blue curve shown in Fig. 8a; these are referred to as Type-2 barriers, in which the Peierls barrier is dominated by the difference in potential energy between the initial and final configurations.
For pure element metals, we naturally expect 100% Type-1 shape barriers since the dislocation dipole energies are constant and the Peierls potential curve will be perfectly symmetric. For instance, based on our drag method calculations, the Peierls barriers in pure element metals, that is, Mo, Nb, Ta, and W, range from 0.12 to 0.38 eV (0.03-0.09 eV/b if normalized by the total The Peierls barriers of pure bcc metals calculated through drag method are also plotted for reference, reproduced from previous DFT study 17 . Burgers vector) in the current simulation geometry; they are plotted on Fig. 6. If we take the highest Peierls barrier value in the pure bcc elements as the reference for the RHEA, it is found that when the core energies follow the Gaussian distribution, the rugged energy landscape and variance in RHEA will inevitably lead to another scenario during the calculation of the dislocation Peierls potentials, in which the final configuration has a much higher potential energy than that of the initial configuration, as shown by the right-side histogram in Fig. 6, which is noted as a Type-2 barrier. Based on the histograms and contour of the Peierls valley energy difference shown in Figs 6 and 7, there is a significant degree of neighboring valley energy differences that have already exceeded the highest Peierls barriers found in pure bcc metals (0.37 eV or 0.09 eV/b in W). For the case of the random solid solution sample s1, as indicated by the green dash-dot line in Fig. 3, which displays a relatively broad distribution of dislocation core energies, the probability of a Type-2 barrier will be higher. However, with progressively increasing SRO in samples s2 and s3, the distribution of core energies narrows, as shown by the blue and red histograms in Fig. 3. The lower variance of the core , and the plotted energy difference corresponds to the difference in energy between the final state (after glide) and initial state for each position of the dislocation cores. The left column contains the 3D contours and right column is the corresponding 2D projection. The contours were plotted by interpolating data points on grids through bivariate spline. energies leads to fewer Type-2 barriers. In what follows, we argue that the variance or standard deviation of the core energies will lead to the asymmetric barriers and the variance itself is affected by the degree of SRO in the materials.
The transition from Type-1 to Type-2 barriers is highly dependent on the relative energy difference between the initial and final dislocation configurations. Here, we assume that for an alloy with a certain level of SRO, the dislocation dipole energy will follow a normal distribution: Normalðμ; σ 2 Þ (similar to Fig. 3), as shown in Fig. 8b. If we assume that there are two random neighboring dipoles: dipole-1 and dipole-2, then dipole-2 represents the initial configuration and will have a preference to glide to its final configuration dipole-1. The energy of these two dipoles are written as E d1 and E d2 . For the transition from a Type-1 to a Type-2 barrier, we postulate that there exists a critical energy difference E critical that when E d1 À E d2 > E critical , the Peierls barrier will become a Type-2. Based on our assumptions for the distribution in dipole energies in Fig. 8b, the energy of dipole-1 and dipole-2 are: E d1 $ Normalðμ; σ 2 Þ, E d2 $ Normalðμ; σ 2 Þ. The energy difference between the two dipoles is then E d1 À E d2 $ Normalð0; 2σ 2 À 2σ cov Þ. Since the energies of two neighboring dislocation dipoles are not independent, we need to consider the covariance σ cov between these E d1 and E d2 values (see details in Supplementary Note 5). Thus, the probability of observing a Type-2 barrier for this condition can be written as: (1) where Φ is the standard Normal cumulative distribution function and σ cov is the covariance between the energy of two neighboring dislocation dipoles.
Based on this equation, the probability of a Type-2 barrier is a function of E critical , the standard deviation σ (or variance) of the dipole energy distribution and covariance between energy of two neighboring dislocation dipoles. In Fig. 8c, we plot P TypeÀ2 as a function of σ for two different values of E critical with σ cov 2 À0:8σ 2 ; 0:8σ 2 ½ . The two E critical values were chosen as 0.4 and 0.6 eV, which is slightly higher than the Peierls barrier calculated in W (0.37 eV). These curves clearly demonstrate that the probability of a Type-2 barrier will increase monotonically with the standard deviation σ, which is also correlated with the state of SRO. For a single screw dislocation, rather than the dislocation dipole geometry considered in this study, we can obtain similar results as P TypeÀ2 ¼ 1 À Φð e critical ffiffiffiffiffiffiffiffiffiffiffiffi  can conclude that in simple terms, the unique variance of dislocation core energies in RHEA, which is also influenced by the SRO, enhances the probability of observing Peierls barriers of Type-2, which will finally influence the dislocation morphologies and their motion.

DISCUSSION
Using first-principles calculations of dislocation energies and the differences in Peierls valley energies in bcc RHEAs, our results reveal fundamental differences between behavior in the multiple principal element alloys and a pure metal or dilute solution. The variation in local chemical environments within the RHEAs leads to a distribution of dislocation core energies for different dislocation segments; moreover, the characteristics of this core energy distribution are significantly influenced by the presence of SRO. In contrast, all the local environments are constant in pure metals and would be expected to show much smaller distributions for dilute solutions.
With our present DFT calculations, although we have doubled the thickness of the sample, the dimension of the out-of-plane direction is still limited to only two Burgers vectors. The calculated core energies and Peierls potentials thus represent the local characteristics of a small straight segment of dislocation line. When considering a long dislocation line gliding in the RHEA, due to the Gaussian distribution of local energies of dislocation segments, described in Fig. 3, the dislocation line will prefer to form a wavy shape to reduce the total potential energy. For alloys with multiple principal elements in equal molar ratios, statistically the composition fluctuation always exists even for a random solid solution.
The Peierls potential plays a crucial role in governing dislocation motion. Here, we have identified two types of Peierls barriers in the bcc RHEA, which depend critically on the energy distribution of the dislocation segments. Considering a long dislocation motion associated with kink-pair theory 72 , it is extremely difficult for some segments gliding through the path of the Type-2 barriers due to its high magnitude. Under such circumstances, these segments can become pinned or are forced to glide on alternative planes or in different directions. This will serve to facilitate cross slip, dislocation multiplication, and the formation of wavy dislocation lines, all of which will eventually enhance the strength and ductility of the material at the macroscale due to homogenization of plastic strains 73 . Indeed, such a form of wavy slip and enhanced mechanical properties has been reported for a bcc TiZrNbHf RHEA with short-range ordered (O, Ti, Zr) complexes 73 . Recently, a theory 74 developed for screw dislocation strengthening in RHEAs has been presented based on the assumption that screw dislocations will naturally adopt a kinked configuration. Along with the MD simulations of the NbTaV alloy 74 , our DFT data, as shown in Figs 3-5, strongly supports the idea that dislocation lines in this and related RHEAs would tend to form a kinked structure.
In summary, we have systematically studied the dislocation core energy, DAPB energy, dislocation dipole energy distribution, and Peierls valley energy differences in a bcc MoNbTaW RHEA using DFT calculations, considering the effects of chemical SRO. Similar to the pure bcc transition metals, compact cores were found to dominate in screw dislocations in the bcc MoNbTaW RHEA. The average core energy of a screw dislocation is higher in the current RHEA compared with the pure bcc transition metals; however, SRO is found to have only a negligible effect on the average equilibrium energy of this line defect. However, the DAPB energy is found to correlate strongly with the SRO state that could potentially influence the dislocation mobility. In addition, the dislocation core energies were found to follow a Gaussian distribution with the increasing degree of SRO resulting in a progressively lower variance of the distribution of core energies.
Resulting from the intrinsic fluctuation of core energies in HEAs, two types of Peierls barriers were discovered, which depend on the difference in core energies between initial and final configurations. By comparison with pure bcc transition metals, the Peierls barrier of screw dislocations in bcc RHEAs is expected to be higher due to the formation of Type-2 Peierls barriers from the variance of the core energy distribution. The findings from the present work highlight the effect of the variance in core energy distributions in influencing dislocation Peierls potentials and suggest important consequences on dislocation morphology and activity, which is an intrinsic feature of HEAs. As these characteristics are heavily influenced by SRO, such local ordering may have a significant impact on the mechanical properties of RHEAs.

Lattice constant determination and simulation cell with dislocation dipole
The lattice constant of the equimolar MoNbTaW HEA was determined by relaxing the 64-atom quaternary quasi-random structure (SQS) 75 provided by Gao et al. 76 . The calculated lattice constant was 3.230 Å and was adopted in all simulations. For the simulation cell with dislocation dipole, we first defined e 1 ¼ a 0 ½112; e 2 ¼ a 0 ½110; and e 3 ¼ a 0 =2½111. Then, the supercell with a dislocation dipole was built with three edges, h 1 = 7e 1 , h 2 ¼ 3:5e 1 þ 5:5e 2 þ 0:5e 3 , h 3 ¼ 2e 3 , to contain 462 atoms. The periodic length along the dislocation line direction, h 3 , was twice the magnitude of the Burgers vector.
DFT-based MC simulations MC simulations were performed using the supercell geometry described above. For the initial condition in these simulations, the sample was generated as an SQS model of the random alloys. The temperature employed in the MC simulations was 500 K. Energy calculations were performed using the Projector Augmented Wave (PAW) method 77,78 , as implemented in the VASP [57][58][59] . A plane wave cutoff energy of 400 eV was employed, and the Brillouin zone integrations were performed using Monkhorst-Pack meshes 79 with a 3 × 1 × 1 grid, where the first index corresponds to the direction along the dislocation line. PAW potentials 78 were employed with the Perdew-Burke-Ernzerhof generalized-gradient approximation for the exchange-correlation function 80 . Lattice MC simulations were then conducted similar to the methods utilized by Tamm et al. 47 and Ding et al. 55 , which included swaps of atom types with the acceptance probability based on the Metropolis-Hastings algorithm 81 . In the current MC simulations, a total of 2094 swaps were conducted and 471 swaps were accepted. For the choice of PAW potentials, 6 valence electrons were used for Mo and W, 5 valence electrons for Ta, and 11 valence electrons for Nb.
Core structure and Peierls valley energy differences Following the MC simulations, the dislocation dipole was introduced into the sample at all possible sites. All configurations with the dislocation dipole were then relaxed through a conjugate-gradient algorithm using VASP with the settings described above, but with a denser k-point mesh of 7 × 1 × 1. Atomic positions were relaxed with a convergence criterion on forces of 10 −2 eV/Å. For each relaxed sample selected as the initial configuration, we chose the sample with a nearest dislocation dipole on the same {110} plane and displaced in the ½112 direction as the final configuration to calculate the valley energy differences between these two neighboring dipoles. Further details can be found in Supplementary Fig. 3.

Local chemical SRO parameter
Similar to the definition described by Ding et al. 55 , which was modified from the Warren-Cowley parameter 66 , we defined the nonproportional number of local atomic pairs, Δδ ij , to quantify the chemical ordering around an atomic species for the combined first and second nearestneighbor shells in the bcc structure, for which the corresponding coordination numbers are N = 14. The value of Δδ ij was then calculated as: S. Yin et al. where N = 14 is the coordination number of first and second nearestneighbor shells in the bcc structure, p ij is the actual probability of bonds between atoms of type j and type i in the sample, p ideal ij is the ideal probability of bonds between atoms of type j and type i for the random solid solution case based on the species concentrations. Δδ ij = 0 for the case of a random solution. The overall SRO is represented by the sum of all the Δδ ij for all species ðSRO ¼ P ij Δδ ij Þ.

DATA AVAILABILITY
The data that support the findings of this study are available from Dr. Sheng Yin (email: shengyin@berkeley.edu) upon reasonable request.