Fluctuations in local shear-fault energy produce unique and dominating strengthening in metastable complex concentrated alloys

Significance Complex concentrated alloys (CCAs) with multi elements in equiatomic ratios are exotic alloys but what governs their strength is currently not clearly understood. This work has proven the prime origin of strength in CCAs as arising from the fluctuations (not the averaged value) of the shear-fault energy, which is an important feature of CCAs unrecognized beforehand. Through identification of the local hard atomic motifs far beyond the traditional concept of the globally averaged short-range ordering parameters, this work has unfolded the physical basis of strength in CCAs and provides important clues for designing high-strength CCAs by focusing on the alloy composition to produce large fault-energy fluctuations.

High-entropy alloys (HEAs) were originally defined as alloys composed of five or more alloying elements with concentrations at 5 to 35% (1) according to the assumption that a high configurational entropy of mixing could stabilize the formation of single-phase solid solutions. However, mixing entropy has significant effects on the Gibbs free energy of an alloy system only at elevated temperatures, and practical alloys are usually homogenized through thermal processing such as annealing or quenching before being used at ambient temperatures. The configurational entropy is also known to distinctly decrease with lowering of annealing temperature (2,3) rather than preserving the ideal random mixing as at high temperatures. Specifically, complex concentrated alloys (CCAs) exhibit local chemical ordering (2,4,5), phase decomposition (6,7), and fluctuations of staking fault energy (2,8) after given thermal processing due to the intrinsic metastable nature in CCAs (9).
Numerous efforts have been made to understand how the complex microstructures mentioned above give rise to the physical basis of strength in HEAs and CCAs (10)(11)(12). Destruction of SRO by dislocation motion has been stipulated as an important contribution to strength in these alloys (13), while current strengthening models are based primarily on classical solid-solution strengthening due to the elastic modulus and size misfit of alloying elements (11,12). Although SRO effects have been considered through taking the energy change of the system as a strengthening contribution when specific chemical bonds are broken (4), the bond-breaking contribution is actually on the basis of mean binding energy and globally averaged SRO parameters, with local fluctuations ignored. In fact, a distinctive feature of HEAs and CCAs is that their dislocations are highly wavy when at rest or traversing, as shown in previous simulations (2,4,14,15) and experiments (16,17). Yet, a comprehensive understanding of such wavy configurations of dislocations and their jumpy motion, and what forms the key basis of strength in CCAs, is still lacking.
Here, we study the above issues by molecular dynamics (MD) simulations on a representative metastable NiCoCr CCA, which is a face-centered cubic (FCC) single-phase solid solution with excellent mechanical properties, especially better than those of the CrMnFeCoNi Cantor alloy (18)(19)(20)(21). SRO in NiCoCr has been confirmed by experiments and atomistic simulations: Local Co-Cr ordering pairs were identified from molecular dynamics (MD) simulations (2) and density-functional theory (DFT) calculations (8), as well as experiments using energy-filtered transmission electron microscopy (5) and extended X-ray absorption fine structure spectroscopy techniques (22). In order to understand the underlying strengthening mechanism in the NiCoCr CCA, here in this work, MD simulations were carried out on species of this alloy with different SRO states generated by Monte Carlo MD at different annealing temperatures. The pinning sites of the dislocations and SRO effects on their waviness and jumpy glide were studied. By accounting for the line tension effects of the wavy dislocation configurations, a form of dislocation resistance unique to CCAs is identified and analyzed.

Jumpy Glide, Pinning, and Waviness of Dislocations in NiCoCr
The MD simulations carried out are summarized in Fig. 1. Here, hybrid Monte Carlo/molecular dynamics (MC/MD) was employed to obtain different alloy species for mimicking different annealing conditions, comprising a random state where the three elements of Ni, Co, and Cr in equal ratio are randomly distributed in the simulation cell, and three other states relaxed from the random state at 1,350 K, 950 K, and 650 K (details in Materials and Methods). The calculated Warren-Cowley parameters characterizing the degree of SRO (23) for different annealed states is given in Fig. 1A, which are comparable to previous results (2). For dislocation simulations, the simulation cells used comprised 20.3 million atoms with the three orthogonal directions as X-110 , Y- [112] , and Z-111 as shown in Fig. 1B, with periodic boundary conditions (PBCs) employed along X and Y, while free boundary conditions were used on Z surfaces for stress application. As the NiCoCr alloy has the FCC crystal structure, its unit ½<110> dislocation is dissociated into two 1/6<112> Shockley partial dislocations separated by a stacking fault with energy SFE or SF , as shown in Fig. 1C. While the FCC crystal structure is recovered by the advancement of the trailing 1/6<112> Shockley partial, for CCAs with SRO, certain chemical disorder will remain in the wake of the trailing partial in comparison with the original SRO ahead of the leading partial; the residual chemical disorder is thus manifested in a trailing antiphase boundary (APB) with specific energy APBE or APB . The MD-calculated energy landscapes at 0 K for the random and the SRO samples are given in Fig. 1C. The first valleys, corresponding to the passage of the leading Shockley partial, correspond to the intrinsic SFEs (24) and complex stacking fault energies (CSFEs) for random and SRO samples, respectively. On decreasing annealing temperature from 1,350 K to 650 K, i.e., increasing the SRO degree, the SFE of −14.4 mJm −2 for the random state increases to 19.2, 33.4, and 39.2 mJm −2 . After the trailing partial passes through the cell, a nonzero APBE remains in all states as in ordered alloys, and as the annealing temperature decreases or SRO increases, the APBE increases from 11.5 mJm −2 to 72.5, 92.7, and 101.2 mJm −2 . However, successive shears of the ½<110> dislocation system will gradually destroy the SRO and hence the APB, as indicated by the fall of the APBE on successive dislocation passes in Fig. 1D.
In previous simulations of dislocations in FCC HEAs (4,15), the entire dissociated ½<110> system comprising the Shockley partials and the stacking fault glides as a whole on the {111} slip plane. However, in the present NiCoCr alloy with negative or weakly positive SFE, the two Shockley partials are pinned and move rather independently when the applied stress is below the critical resolved shear stress (CRSS) of the full dislocation, resulting in a spectrum of metastable dissociation distances as analyzed by Vaid et al. (25). Fig. 1E shows the glide distances of the two partials at different applied stresses at 5 K for the sample annealed at 950 K. In this annealed state and others, both partials show jerky movements under a smaller applied stress, and as the applied stress increases, the leading partial is first unpinned and continuously glides through the periodic image of the simulation cell (the slip distance of the leading partial in the periodic image is not plotted in this figure and in SI Appendix, Figs. S2 and S3) before the trailing partial starts to move. This implies that the stable separation between partials could be spanning nearly half of the simulation cell due to the low or negative SFE. As the applied stress increases up to the CRSS of the full dislocation, the entire extended dislocation begins to move. The CRSSs of the two partials at 5 K in different annealed states are given in Fig. 1F, indicating that the CRSS of the trailing partial is larger than that of the leading partial. On the first pass of the ½<110> dislocation system, the predicted CRSSs of both partial dislocations show significant strengthening attributed to the increase in SRO with decreasing of the annealing temperature, in consistence with the trends of the CSFEs and APBEs in Fig. 1D and previous simulations (2). However, after a few passes of the ½<110> dislocation system, the strengthening due to the SRO gained at lower annealing temperatures diminishes, as shown in Fig. 1F. In particular, there is a softening of ~100 MPa, compared with the random state, for the leading partial.
The CSFEs and APBEs shown in Fig. 1C are global values averaged over the slip plane, and they serve to illustrate the contribution of ordering to strengthening. However, the planar fault energies of CCAs are locally varying (2,8,16) over the slip plane due to the spatial variations of the alloying elements and chemical ordering. Fig. 2 A-C show the statistical distributions of the CSFE and APBE and Δ = APBE − CSFE for the different annealed states studied, respectively. With increasing SRO degree, the mean CSFE and APBE increase in consistence with the results in Fig. 1C, and correspondingly, the mean Δ also increases with the SRO degree while its fluctuations decrease. The local values of the CSFE and APBE and Δ , respectively, on the {111} plane for the example of the alloy state annealed at 650 K are given in Fig. 2 D-F, where the cross-sectional area used for local fault-energy calculations is 15 × 17 Å 2 . Successive snapshots of the leading and trailing partials in time intervals of 6 ps and 14 ps, respectively, are superimposed on the local fault-energy maps, as shown in Fig. 2 D and F, as they traverse from the left of the simulation cell to the right. The motion of the partials is evidently jerky, which corresponds well to the jumpy flow in Fig. 1E, and between the jumps, the dislocations are pinned at locations with high values of the corresponding local fault energy (examples marked by double-arrowed lines).  Fig. 2I, indicate that the pinning points that correspond to high local fault energies are predominantly Co-Co-Cr-Cr, Co-Co-Co-Cr and Ni-Cr-Co-Co tetrahedrons, which we call "hard atomic motifs" (HAMs) here.
Therefore, to summarize, the present MD simulations indicate that dislocation motion in the NiCoCr CCA is jerky, with the dislocations adopting highly wavy configurations due to pinning at sites with high local fault energies. This aspect will be analyzed in detail in the next section. Fig. 3A shows the various Peach-Koehler forces (forces on unit length of dislocation) acting on the two 1/6<112> Shockley partial dislocations in the present alloy system with SRO contributions. Let us first consider the case of mean-field mobility, i.e., the average mobility of a longenough dislocation with local fluctuations along the dislocation length ignored. The motion equations are given as follows:

Mean-Field Mobility of Dislocations.
, m and v are the inverse of mobility and velocity of the dislocation, and all quantities refer to the averages along the leading (L) or trailing (T) partial. When v L,T > 0 , F L,T f reaches a critical value F c f which is the kinematic lattice friction. When v L,T = 0 , i.e., the dislocation is not moving, the lattice friction F L,T f can be any value between ±F c f to balance the other forces acting on the dislocation. Two cases are imminent: Case I-when the leading partial just starts to move and the trailing partial is still stuck. Then, Eqs. 1 and 2, which is the CRSS to move the leading partial. Case II-when the trailing partial starts to move, and the leading partial is already moving. Then, In Eq. 3, the friction F T f is a changeable quantity called to balance the forces acting on the trailing partial before it begins to move, and as such, it is dependent on the initial separation of the two partials. However, a T from Eq. 4 is a well-defined value as APB and F c f are constants for a given alloy. Since Eq. 3, the lower and upper bounds for a L are as follows: We first assume F c f as a constant for the random and annealed samples. The MD-calculated CRSSs and average APBEs from  6 for the Shockley partial, where a = 3.557 Å is the lattice constant from the present MD simulation which is also comparable to the experimental value of 3.575 Å (26). The intercept of the fitting equation is 675 MPa, which, for a mean-field straight dislocation as described by Eq. 4, would be the Peierls stress F c f ∕b of the trailing partial in the random state of the alloy. However, it is important to note that for the CRSS data predicted by MD in Fig. 3B, the dislocations are highly wavy as seen in Fig. 2, and so, instead of the mean-field Peierls stress, a significant contribution to the 675 MPa actually arises from the line tension of the wavy dislocations, as discussed in the next section. In any case, the linear increasing trend of the MD-predicted CRSS with the APBE in Fig. 3B is in good agreement with Eq. 4, and this confirms that the mean-field dislocation mobility is hampered by APBE strengthening due to the SRO in the annealed states compared to the random state. However, as can be seen from Fig. 3B, the amount of APBE strengthening is only about 300 MPa for the state with the highest SRO (the 650 K annealed state) compared to the much higher basic strength of ~680 MPa for the random state. Moreover, as shown in Fig. 1D ∆γ (x,y) quickly on successful dislocation passes, and so, the APBE strengthening demonstrated here is actually unimportant for macroscopic slip. The key understanding of strength in CCAs would still require an understanding of the basic strength of the random alloy state, which is pursued below. Fig. 2, dislocations are highly wavy in the present alloy due to interactions with sites of high local fault energies and special HAMs. Furthermore, as can be seen from the Insets of Fig. 2 A and C, the random state of the alloy has the largest fluctuations in the local fault energy. Here, the contribution to the CRSS based on the local fault-energy fluctuations is deduced. Fig. 3C shows the jump displacement w y of a wavy dislocation between two dwelling states and the various forces acting. Balancing the work done as the dislocation undergoes the jump gives

Mobility of Wavy Dislocations. As indicated by the results in
Here, the first term is the work done by the net force F a − F f which is equal to b , where is the resistance due to the stacking fault; the second term is the change in self-energy (line tension) of the dislocation due to its change in line length during the jump, and the third term is the work done against the fault energy which is in general a field x, y over the slip plane. Line tension T may be written as T = b 2 , where is shear modulus, and depends on the dislocation character inner and outer cutoff radii of the dislocation.
The dislocation waviness is due to the fluctuating x, y field where local regions of high pin the dislocation (resulting in zero or small w ), as further advancement of the dislocation would produce faulting with a high energy (Fig. 2). Therefore, the line-tension term is also a consequence of the nonuniform x, y field, so that both terms on the right side of Eq. 6 are due to x, y . The x, y field is spatially fluctuating as schematically shown in Fig. 3D, and can be written as follows: [7] x, y = + Δ x, y , where is the average fault energy over the slip plane (hence it is no longer a field of x and y), and Δ x, y is the fluctuating component of x, y that averages to zero. The fault-induced resistance above may therefore be considered to be composed of two corresponding parts: = − + Δ , where, from Eq. 6, is the old term in Eq. 3 for a "mean-field" dislocation, and Δ is given by As the diminishes on successive dislocation passes as shown in Fig. 1D, − from Eq. 8 would also disappear gradually on successive dislocation passes as shown in Fig. 1F, and so, this term is unimportant for macroscopic slip. However, the fluctuation Δ x, y should always remain and so should be the Δ from Eq. 9. The fluctuating Δ x, y is a unique feature of HEAs and CCAs, and so, the Δ is also a unique resistance for CCAs.
Of the two terms on the right side of Eq. 9, we note that ∫ L 0 ∫ w(y) 0 Δ (x, y)dx dy ∕ ∫ L 0 w y dy is the average value of Δ (x, y) in the slipped region, which, for a long-enough dislocation, should also approach zero. Hence, Eq. 9 becomes Δ (x, y)dx dy averages to zero for a long-enough dislocation, the line tension term in Eq. 10 is still due to the fluctuating displacement w y of the dislocation and hence the fluctuating Δ x, y field.
A general and versatile representation of the jump displacement w y is a Fourier series of the form: Fig. 3E shows an example of accurate fitting to MD-computed dislocation shape). Substituting Eq. 11 into Eq. 10, as the domain length L → ∞, Here, ∫ L 0 cos 2 2πy Hence as L → ∞ , only the factor L/2 prevails in the first summation in the second step of Eq. 12, which leads to the third step. Similarly, Therefore, Eq. 10 reduces to is the aspect ratio, with i and i being the bowing amplitude and length of the i th harmonics of the dislocation jump displacement.

Analysis of Alloy Strength
As mentioned in section 2, the partial dislocations glide separately rather than as a whole initially when the required stresses for their movement reached. As the applied stress increases, the leading partial continuously glides before the trailing partial starts to move, therefore, for the leading partial, a repulsive force arising from the unslipped image trailing partial gradually increases. Correspondingly, a pulling force arising from the gliding image leading partial acts on the unslipped trailing partial. This is actually an artificial strengthening arising from the periodic image dislocations due to the periodic boundary conditions applied. Besides, the plasticity in the present NiCoCr alloy is accommodated through the glide of partial dislocations, the CRSS of which is xz cos(30 • ) , where xz is the applied stress in the <110> direction on {111} planes. Therefore, to allow comparison with experiments, all MD and CRSSs from the simulations should be corrected as, where "+" applies to trailing partials and "−" to leading partials.
is the distance between the gliding leading partial and the unslipped image trailing partial, shear modulus G = 59 GPa, Poisson's ratio = 0. 3 , and thus, image ≈ 61MPa. Next, we consider the lattice friction or Peierls stress. Zhang et al. (27) incorporated lattice distortion effects into the classical Peierls-Nabarro model for the lattice friction by scaling the interplanar potential (or generalized stacking fault energy) by a factor that is stochastically distributed according to the normal distribution with mean 1 and SD Δ. Their model predicts the Peierls stress to be [16] Fig. 3B. The large gap between the MD-predicted CRSS and the Peierls stress can be filled by the strengthening Δ arising from the fault-energy fluctuations as proposed in Eq. 14. To obtain the i and i in the Fourier series in Eq. 11, the atomic positions of a wavy partial dislocation were first fitted into a curved line using a smoothing algorithm (28,29), and then, the computed w i (y), being the difference between two successive dislocation shapes X i (y) and X i+1 (y) over the segment length L as shown in Fig. 3C, was fitted to a Fourier series as shown in Fig. 3E. As shown in SI Appendix, Fig. S6, the calculated fault-energy fluctuation resistance is still a distribution, which arises from the spatial variance of the hard/soft zones on the slip plane. For both the leading and trailing partials, there is remarkable contribution to the CRSS, namely, the Δ values are on the order of 300 to 400 MPa, respectively, as shown in Fig. 4A. It can be seen that the Δ is pretty constant with a slight decreasing trend with increased annealing temperature that tallies with the trend of the fluctuations of the fault energies in Fig. 4B. Furthermore, sample calculations of Δ were also performed using Eq. 10 without involving the approximations leading to Eq. 14, but the means and SDs are very similar to those calculated directly using Eq. 14 (c.f. SI Appendix, Figs. S6 A and E and S7 A and B). These confirm that fluctuations of fault energies play a crucial role in the strength of such CCAs.
Up to this point, as shown in Fig. 4C, the MD-predicted CRSS for the annealed alloy states can be interpreted as the summation of the various strengthening contributions discussed above, namely, Peierls stress ( p ) and fault-energy fluctuation resistance ( Δ ), together with the SRO strengthening ( APB = APB ∕(2b) , from Eq. 5) which is present for the annealed alloys. However, in terms of magnitude, it can be seen from Fig. 4C that the most significant contribution is Δ , which makes up 49 to 55% of the total CRSS. Moreover, as mentioned above and illustrated in Fig. 1 D and F, the SRO strengthening ( APB ) which appears to be significant is only important for the first dislocation pass and will subdue on further dislocation passes as in macroscopic slip. Therefore, even for an annealed alloy state with initially large SRO strengthening, the real strength in macroscopic deformation is still dominated by the Δ contribution which remains on successive dislocation passes (Fig. 1F).
The various strengthening contributions are also compared with experiments to understand their importance. The reported experimental CRSS of a NiCoCr single crystal at 14 K is about 170 MPa (30). Considering that the samples in the experiments conducted by Li et al. (30) were annealed at 1,473 K and were close to the random solid-solution state, we took the random state in the MD simulations for comparison. Combining the mean Δ values of the random state (393 MPa and 479 MPa for the leading and trailing partials, respectively) with the Peierls stress ( p = 98 MPa), the CRSSs are 491 MPa and 617 MPa (including APB = 40 MPa due to the small APB of 11.5 mJm −2 for the random state) for the leading and trailing partials, respectively, which are in good agreement with that obtained from the MD simulations as shown in Fig. 4D. A caveat needs to be noted when comparing the present MD-simulated CRSS to the experimental CRSS. The MD CRSS is calculated for a single dislocation in the simulation cell, but an experimental CRSS measured from a macroscopic sample is for a slip system comprising many dislocations. In the latter case, the group interactions of dislocations will cause a portion of the dislocations to start moving at an applied stress smaller than the Peierls stress of the individual dislocations, and then, this will bring in stress concentrations causing dislocation avalanches and sudden strain bursts (31)(32)(33)(34). This effect is thought to explain why MD-calculated CRSSs at 0 K are about three times larger than the experimental values (35,36), which is also confirmed by discrete dislocation dynamics simulations (37) with a large range of initial dislocation densities, i.e., 10 12 to 10 16 m −2 , that are close to general experimental conditions. In the light of this effect, we apply a correction factor of 1/3 to the MD-simulated CRSSs for single dislocations; the results are 166 MPa for the leading partial dislocation and 232 MPa for the trailing partial dislocation, which are well in line with the experimental CRSS of NiCoCr single crystal at 14 K, i.e., ~170 MPa (30), as shown in Fig. 4E. Since p is small, either the Δ alone or the superimposed ( p + Δ ) compares rather favorably with the MD results in Fig. 4D, from which the CRSS obtained by applying the 1/3 scaling factor is in accordance with the experimental results as shown in Fig. 4E, thus confirming the importance of the Δ strengthening. We note in passing that a solid-solution model (11,38,39) also exists in the literature that serves to predict strength in HEAs by taking into account dislocation line-tension effects and size mismatch as follows: where the CRSS at 5 K was approximated as that at 0 K, = 0.125, = 0.02063, and = 0.268 (38) and G V = 107 GPa is the Voigt shear modulus as listed in SI Appendix, Table S1. The calculated macroscopic critical stress for a slip system through the above 1/3 scaling correction is 42.6 MPa as shown in Fig. 4E, which is much smaller than the corresponding experimental results.

Discussion and Conclusion
The simulations and analysis presented above have identified a dominating strength contribution due to fault-energy fluctuations using an EAM interatomic potential for the NiCoCr CCA. Previous understanding of the interactions between dislocations and pinning sites in CCAs was based primarily on solid-solution strengthening models such as that in Eq. 17, which is inadequate for the present studied alloy as indicated by the comparison in Fig. 4E. This probably is because such classical models involve a number of hard-tojustify assumptions, such as equal energy of dislocation segments (hence neglecting irregular positioning of pinning sites on the global slip plane) and neglecting chemical interactions between solutes and dislocation cores. Unlike such models which give a local, characteristic unpinning energy barrier and stress (based primarily on a representative pinned segment and bow-out amplitude), the approach in this work does not rely on any strong assumptions on the energy and spatial positioning of the pinning sites; instead, it simply uses the line tension description of dislocations in conjunction with very general dislocation shapes representable by the Fourier series to fit to MD results to obtain the dislocation resistance (c.f. Eq. 14). In the present approach, all elastic and chemical nature of the interactions is embodied in the MD simulation data for the dislocation shapes and does not appear as input parameters as in the classical solid-solution strengthening models. The result is a global resistance Δ based on a statistical ensemble of dislocation waviness in terms of i and i obtained directly from MD simulations without any strong assumptions.
The waviness of dislocations was experimentally observed by weak-beam transmission electron microscopy in CCAs including the NiCoCrFeMn alloy 17 (FCC) and Ti 38 V 15 Nb 23 Hf 24 alloy (BCC) (40). An important question to answer, at this juncture, is whether the strength contribution Δ is unique to the interatomic potential used in the present study. To address this, we have also employed other potentials to simulate the same NiCoCr ternary alloy and the quinary NiCoCrFeMn alloy. For NiCoCr, we have also used a recently developed machine-learning potential (41) to calculate the SFE distribution for the random alloy state, and as shown in SI Appendix, Fig. S8, the SFE distribution is similar to that obtained from the present EAM potential, i.e., there is a distribution varying from −170 to 155 mJm −2 (c.f. −120 to 100 mJm −2 from the EAM potential) with mean SFE −32 mJm −2 (c.f. −16 mJm −2 from the EAM potential). However, the machine-learning potential is very costly to use-its use has only been demonstrated for systems of several thousand atoms (41,42) which are too small for dislocation simulations. In addition, we also employed an effective pair potential developed by Groger et al. (43) to perform similar dislocation mobility simulations for the random NiCoCr alloy. The predicted CRSS at 5 K after correction by multiplying the factor of 1/3 is around 90 MPa which is much lower than the experimental values given in Fig. 4E, and as shown in SI Appendix, Fig. S9A, the pair potential predicts a positive mean value for the SFE with very small SD, and the dislocations are very straight during slip. The results suggest that the pair potential is not satisfactory for the NiCoCr alloy, probably due to the fact that it does not account for multiatomic interactions which are important to capture the fluctuation of SFE as confirmed by multiple DFT reports (8,44,45). Without the SFE fluctuation properly captured, the strengthening contribution from SFE fluctuations or wavy dislocations would be artificially suppressed, which may explain for the significant underprediction of the CRSS by the pair potential. Employing a similar pair potential to the quinary alloy NiCoCrFeMn (19,46), the predicted CRSS after the similar correction above is 103 MPa (the critical stress for a single dislocation is about 310 MPa as shown in SI Appendix, Fig. S9B), which is lower than the experimental single-crystal CRSS at 0 K, i.e., 168 MPa as reported by Kawanura et al. (47). As shown in SI Appendix, Fig. S9B, the SFE fluctuation is rather weak (at ~ 8 mJm −2 ) compared with that in the NiCoCr alloy (at ~35 mJm −2 as given in Fig. 4B), with also much less wavy dislocation shapes as simulated by the pair potential that takes no consideration of multiatomic interactions. Therefore, to round up, the SFE fluctuations calculated from the EAM potential for the NiCoCr alloy are consistently reproduced by another, machine-learning, potential (SI Appendix, Fig. S8). However, the choice of the potential does affect the prediction-as in the case of ternary NiCoCr, the use of a pair potential would significantly underpredict the CRSS as it artificially suppresses the SFE fluctuation and dislocation waviness (SI Appendix, Fig. S9A). In a more general context, while the strength contributor of Δ is predicted by the present EAM potential to be dominating in the ternary NiCoCr, its significance in other alloys and any governing rule behind will be interesting topics to investigate in the future. In particular, how this strength contributor depends on the alloy complexity and chemistry can be investigated as a data-mining exercise after reliable interatomic potentials have become available for a wider range of CCAs.
Recently, Ostesky and Morris (48) have shown that partial dislocation cores in NiFe alloy tend to favor certain chemical environments, and this effect may manifest itself as local obstacles against dislocation glide. In the present simulations, there was no atom swapping during motion of the dislocations thus, the composition of dislocation cores depends only on the instantaneous location of the dislocation, and the resistance arising from the HAMs depends completely on the spatial distribution of the SFE. From both Eqs. 1 and 2, one can see that the pinning sites or HAMs correspond to regions of high fault energies, which is confirmed through the local fault-energy statistics as shown in Fig. 2 G and H. However, atomic exchange via short-circuit diffusion near the dislocation cores may add to the SFE fluctuations, and further work may investigate this type of dynamic obstacles using additional Monte Carlo atom swapping during dislocation motion.
The motivation of the present work is to study single phases of solid solutions with different degrees of SRO, and for this reason, the NiCoCr alloy was annealed at different temperatures to obtain different SRO states. The strongest SRO state was obtained from annealing at 650K, and after MC exchange, the obtained phase was still a single phase. At this low annealing temperature, it is possible that the truly equilibrium alloy state is two phase, although further increasing the MC steps to 8.87 million did not result in significant changes of the potential energy, alloy appearance, and APBE, as shown in SI Appendix, Fig. S10 A and B, and in agreement with previous simulations using the same potential (2), Co-Cr ordering and Ni-Ni clustering are obvious. Furthermore, the computed APBEs after five passes of the <110>a/2 unit dislocation do not exhibit a dramatic jump as reported in ref. 2 (from −50 to 50 mJm −2 ), indicating the stability of the alloy state. The truly equilibrium state of a strongly ordered alloy may be difficult to achieve using MC/MD simulations from small systems, and it is likely that the present alloy state at 650 K is metastable rather than truly stable. In general, if ordering is strong enough to result in a two-phase microstructure, resistance to dislocation motion will be more complicated, and this aspect deserves further investigation.
In conclusion, in the present work, we conducted MD simulations and deduced a model to unfold the physical basis of strength in the metastable NiCoCr alloy at 5 K. The Shockley partials glide in a characteristic jumpy and wavy manner due to successive pinning at sites with high local shear-fault energies associated with special, hard atomic motifs (HAMs) of predominantly Co-Co-Cr-Cr, Co-Co-Co-Cr, or Ni-Cr-Co-Co tetrahedrons. For different SRO states of the alloy achieved by different annealing temperatures, there is an initial strengthening due to the initial fault energy (i.e., CSFE and APBE), which increases with the SRO. However, the APBE and the associated strengthening subdue quickly after multiple passages of dislocations and so are unimportant in macroscopic slip. The basic CRSS which is significant in the random solid-solution state is dominated by the dislocation resistance due to the fluctuations in the local shear-fault energy, which will always remain in the alloy. This strengthening mechanism is unique in CCAs as it is caused by the local fluctuations in the shear-fault energy, corresponding to the HAMs mentioned above, that result in wavy shapes of dislocations of high line energies. The strengthening provided by this mechanism is nearly the entirety of the random alloy strength as the Peierls stress contribution due to lattice misfit distortions is very small. The present identification of the local HAMs as locations of high local shearfault extends the understanding of dislocation-lattice interactions in CCAs far beyond the traditional concept of the globally averaged SRO parameters. The findings from the present study highlight the effect of fault-energy fluctuations on dislocation waviness and confirm that the fault-energy resistance to dislocation slip is an intrinsic feature of CCAs. The present study not only reveals the unique strengthening mechanism in CCAs but also provides important clues for designing high strength CCAs through tuning the alloy composition to produce large fault-energy fluctuations.

Materials and Methods
The MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (49) within the three-element Ni-Co-Cr interatomic potential developed by Li et al. (2). Local lattice and defects were identified through common neighbor analysis (CNA) (50) and the centrosymmetric parameter, respectively. The Open Visualization Tool (OVITO) was adopted to visualize dislocation structures (51). √ 3a (where a is the lattice constant of 3.557 Å) along X-110 , Y-112 , and Z-111 respectively. Periodic boundary conditions (PBCs) were employed along all three orthogonal directions of X, Y and Z. The simulation cells were then heated up to annealing temperatures of 650, 950, and 1,350 K using the Verlet algorithm and then relaxed under the desired temperatures with the NPT (the number of particles N and pressure P and temperature T are constants) ensemble for 50 ps to reach thermal equilibrium (14). After that, the Monte Carlo (MC) scheme was used to swap atoms inside the cells as per an acceptance criterion of energy reduction under the NPT ensemble at the desired temperatures. To do this, the MC swapping was attempted up to 100 times for every 50 MD steps with a time step of 2.5 fs under the NPT ensemble (4). As shown in SI Appendix, Fig. S1 A and B, both the acceptance rate of MC swapping and potential energy of the annealing samples achieve complete convergence after 6 million MC steps. The calculated Warren-Cowley parameters characterizing the degree of SRO (23) for different annealing temperatures are given in Fig. 1A, which are comparable to previous results (2).

Simulations of Dislocation Slip.
The equilibrated configurations after exhaustive MC/MD swapping were then replicated nine times along the Z-1, 1, 1 direction to generate a series of 30-nm-long simulation cells each comprising 20.3 million atoms as shown in Fig. 1B for use in the subsequent simulations for the CRSS of the ½<110> edge dislocation in the NiCoCr alloy. The edge dislocation was introduced parallel to the Y-112 direction as shown in Fig. 1B according to its anisotropic elastic displacement field, and the simulation cells were subsequently relaxed via a conjugate-gradient scheme to obtain equilibrated dislocation structures comprising two 1/6<112> Shockley partials separated by a stacking fault (14). From this point onward, the hitherto PBCs applied on the Z surfaces (for MC annealing) were changed to free boundary conditions for stress application, while the PBCs were maintained in the other two directions X and Y. To minimize the effects of thermal fluctuations on the dislocation glide and mimic the release of elastic energy during dislocation motion (52), the cells were heated up to 5 K in a quasi-static manner with the NVT ensemble for 50 ps. In order to determine the CRSS, external constant forces were gradually imposed on two thin slabs of thickness 1.2 nm at the lower and upper Z surfaces to produce a desired shear stress xz through an incremental loading approach (53). Initially, a larger stress increment of 100 MPa was used to quickly identify a rough CRSS, and then, a smaller stress increment of 20 MPa was applied to figure out the more exact CRSS. To minimize shock wave effects, the whole simulation cell was linearly strained under the applied stress, and then, the strained cell was relaxed except for the atoms in two surface thin slabs being fixed. After that, the corresponding external stress was directly applied on the lower and upper thin slabs in which atoms were able to move to respond to the applied stress, and at the same time, relaxation was carried out using an NVT thermostat at 5 K for 200 ps (14,15).
As indicated in Fig. 1B, the length of the simulation cell in the Y direction (the dislocation line direction) was about 600 Å. To ensure that this is sufficiently long to not affect the calculated CRSS (54), convergence tests with shorter lengths of ~400 and ~500 Å were carried out, and it was found that the CRSS results were very similar. The simulations presented in this work were therefore based on the length of about 600Å (~ 240b) along the Y direction. Data, Materials, and Software Availability. All study data are included in the article and/or SI Appendix.