Multi-scale Investigation of Chemical Short-Range Order and Dislocation Glide in the MoNbTi and TaNbTi Refractory Multi-Principal Element Alloys

Refractory multi-principal element alloys (RMPEAs) are promising materials for high-temperature structural applications. Here, we investigate the role of chemical short-range ordering (CSRO) on dislocation glide in two model RMPEAs - TaNbTi and MoNbTi - using a multi-scale modeling approach. A highly accurate machine learning interatomic potential was developed for the Mo-Ta-Nb-Ti system and used to demonstrate that MoNbTi exhibits a much greater degree of SRO than TaNbTi and the local composition has a direct effect on the unstable stacking fault energies (USFE). From mesoscale phase-field dislocation dynamics simulations, we find that increasing SRO leads to higher mean USFEs, thereby increasing the stress required for dislocation glide. The gliding dislocations experience significant hardening due to pinning and depinning caused by random compositional fluctuations, with higher SRO decreasing the degree of USFE dispersion and hence, amount of hardening. Finally, we show how the morphology of an expanding dislocation loop is affected by the applied stress, with higher SRO requiring higher applied stresses to achieve smooth screw dislocation glide.


Introduction
Refractory multi-principal element alloys (RMPEAs), which are composed of refractory metals in highly concentrated proportions, have garnered intense interest as promising candidate materials for high-temperature structural applications [1,2].Most RMPEAs studied to date possess superior yield strengths and several RMPEAs, such as TiZrHfNbTa, TaNbTi, and HfNbTa, have reported good tensile ductility as well [3][4][5][6].To accelerate design of structural RMPEAs within the vast composition space, it is necessary to understand the dislocation-based mechanisms driving their response to applied load.RMPEAs are chemically disordered at the atomic scale leading to unusual dislocation behavior.Several molec-ular dynamics (MD) simulations and transmission electron microscopy studies alike have pointed to stochastic glide, slip on higher-order glide planes, screw dislocation cross-kinking, and relatively low edge dislocation mobility at higher frequencies in RMPEAs than expected in pure refractories [7][8][9][10][11][12].
In most studies and analyses, the chemical disorder was treated as being ideally uniformly random exhibiting no thermodynamically driven chemical short-range order (CSRO).While this may be a suitable approximation at high temperatures or for the bulk, some preferential pairing between elements in the alloy, or SRO, can be expected to manifest, especially over length scales corresponding to a dislocation.In recent years, detectable levels of SRO have been reported in both FCC MPEAs such as CoCrNi [13,14] and CoNiV [15] as well as the BCC RMPEAs such as HfNbZr [16].Supporting calculations of SRO in MPEAs have been provided via density functional theory (DFT) and hybrid Monte Carlo (MC)/MD simulations, where temperature is used to adjust the extent and severity of SRO [17][18][19][20][21][22][23][24].DFT calculations have shown that SRO in MoNbTaW RMPEA reduces the variation in the dislocation core properties, while have little effect on their core energies [25].MC/MD simulations of nanocrystalline FCC and RMPEAs predict a preference for greater degrees of SRO near the grain boundaries than the interiors and overall greater strength [21,26].In MD simulations of high-velocity dislocation motion, SRO tended to increase edge dislocation mobility but decrease screw dislocation mobility in RMPEA [27], while the partial dislocation mobilities were decreased in an FCC MPEA [21].
In this work, we seek to elucidate the extent of CSRO affects the motion of dislocations, in both mechanically under-driven and over-driven situations.We introduce a multiscale modeling approach that links a new machine learning potential with quantum fidelity to hybrid MC/MD for temperature-dependent CSRO calculations at the atomic scale and ultimately to dislocation dynamics simulation of long dislocations and dislocation loops at the mesoscale.The interatomic potential presented here is a highly accurate one developed for the Mo-Ta-Nb-Ti system [27][28][29][30][31].With it, two base ternary systems, MoNbTi and TaNbTi, are extracted and studied via hybrid MC/MD for CSRO at two annealing temperatures.We elected these MPEAs since experimental observations find that their equimolar forms exhibit disparate mechanical properties, with the MoNbTi being substantially greater in tensile yield strength, peak strength and strain hardening than TaNbTi.To study the effect of CSRO on the dislocation glide mechanisms, a method is developed to build statistical instantiations of large 3D crystals with spatial mapping of CSRO and associated unstable stacking fault energies (USFEs).These crystals are readily analyzable by the real-space 3D phase-field dislocation dynamics (PFDD) technique, which predict stress-driven pathways taken by individual dislocations.With the multiscale strategy, we show that SRO strengthening manifests in both MPEAs, with the average USFEs and critical stresses to initiate and sustain propagation of dislocations increasing with SRO above those for the ideal random solid solution.We introduce a figure of merit to measure SRO impact across compositions and annealing treatments and with it, reveal that CSRO strengthening contribution scales linearly with degree of SRO.The computations also reveal that gliding dislocations in subcritical conditions experience significant hardening.This glide hardening is strongly correlated to the statistical dispersion in the local USFE, and since CSRO tends to narrow the distribution in USFE, glide hardening decreases with degree CSRO.In studying dislocation loop expansion across stress regimes, a transition between jerky and smooth dislocation glide is identified and related to stress sensitivity of kink-pair nucleation rates of the screw character portions.Finally, analysis reveals that initially screw-and edge-oriented dislocations will become wavy in glide yet move via different mechanisms-kink-pair formation and migration vs. pinning/depinning.Their individual glide mechanisms do not change with composition, amount of SRO, glide distance, or subcritical or overdriven loading conditions.These computations explain why MoNbTi is the stronger one and has the greater strain hardening and forecasts that it is more amenable to SRO strengthening.

Moment tensor potential
Studies of the RMPEAs here are enabled by the development of a highly accurate machine learning interatomic potential for the Mo-Ta-Nb-Ti system based on the MTP formalism [28][29][30][31].Fig. 1a provides an overview of the MTP fitting procedure, which is based on a well-established workflow developed by some of the current authors in previous works [26,27].To best investigate the effect of composition variations on SRO and dislocation glide, the training data were carefully selected to encompass all known unary, binary, ternary and quaternary phases in the Mo-Ta-Nb-Ti system (see Methods for details).Fig. 1(b,c) show that extremely low test mean absolute errors (MAEs) were achieved for energies (4.1 meV•atom −1 ) and forces (0.067 eV• Å−1 ), comparable to that achieved previously for the NbMoTaW RMPEA [26,27].The MTP also reproduces very well the DFT elastic constants for the constituent elemental systems, as shown in Fig. ??.The shear moduli µ are 29.6 and 32.3 GPa for TaNbTi and MoNbTi, respectively, and the Young's moduli are 82.7 and 90.7 GPa, respectively.

Chemical short-range order
To calculate the chemical SRO, bulk BCC supercells with equimolar MoNbTi and TaNbTi as well as non-equimolar ternaries with elemental ratio of 3:4:4 and 3:1:1, i.e., X 4 Nb 3 Ti 4 , X 4 Nb 4 Ti 3 , X 3 Nb 4 Ti 4 , X 3 NbTi, XNbTi 3 , XNb 3 Ti, where X = Mo or Ta, were constructed.For each composition, three levels of SRO are achieved by studying the as-constructed random solid solution (RSS) and annealing at 300 K and 1673 K using MC/MD simulations with the MTP.The SRO for the final equilibrium structures is characterized using Warren-Cowley parameters, which are defined for each pair type ij as where p ij is the probability of having an atom type j in the nearest neighbor shell of an atom type i and c j is the concentration of type j and δ ij is the Kronecker delta function [32].By definition, for a RSS, α ij ≈ 0 and for greater degrees of SRO, the absolute value of α ij increases.Fig. 2 shows the Warren-Cowley parameters for the annealed RMPEAs.For all compositions, greater levels of SRO can be achieved with the lower annealing temperature.For the same set of elements, greater degrees of SRO can be accomplished with off-equimolar stoichiometry, i.e., the 3:1:1 and 3:4:4 compositions, than equimolar.In materials annealed at 300 K, the SRO exhibited by RMPEAs containing Mo (Mo x Nb y Ti z ) are much higher than those containing Ta (Ta x Nb y Ti z ).These two types of RMPEAs would be expected to respond differently to the same processing condition or heat treatment, with MoNbTi being much more susceptible than TaNbTi.

Unstable stacking fault energy
Fig. 3 plots the calculated unstable stacking fault energy (USFE) on the {110} plane, shifting along with 111 directions for as-constructed RSS and samples equilibrated at two different annealing temperatures.Fig. ?? A-D plots the USFE distributions for all the RMPEAs with different compositions.In all cases, a higher concentration of Ti reduces the USFE.In the Mo-Nb-Ti system, greater concentrations of Mo increases the USFE, in agreement with a prior work using another interatomic potential [33].On average, for the same composition, annealing at 300 K raises the USFE indicating that SRO increases USFE.To correlate the USFE with its local composition, the local composition is determined based on the composition of the first nearest neighbor planes surrounding the cleaving plane.The relationships between USFE and the local composition are plotted in Fig. ?? E. The correlations between the USFE of a plane and its local composition are consistent with those observed for average USFE of different bulk concentrations.Higher local fractions of Mo significantly increase the USFE, while higher fractions of Ta also increase the USFE but not as significantly as Mo.Higher fractions of Ti significantly decrease the USFE.The trend of USFE with local composition is consistent with the trends observed for the bulk composition as discussed above.

Phase-field framework
PFDD was employed to simulate the behavior of dislocations in both alloys (see Methods, [34]).All inputs into PFDD are physical, such as the elastic constants, lattice parameters, and USFE, and are obtained from the MTP calculations in the preceding sections.In an MPEA, the resistance to dislocation glide is not homogeneous within the material due to variations in local composition [7,8,[35][36][37].A dislocation will only be affected by the type and arrangement of atomic elements within a radius of a few Burgers vectors [38].As such, at any position and point in time, the dislocation only samples a small nanoscale region of the crystal with a local composition differing from that of the bulk.In PFDD, the local composition of every nanometer region is associated with its corresponding USFE.We account for these spatial variation in composition by creating a crystal with spatial variation USFE, simulating the variable resistance to dislocation glide [39,40].
To generate position-dependent USFE crystals for MoNbTi and TaNbTi, random lattices with target SRO levels were created using a binary alloy swapping method adapted from Gehlen et al [41] (see the Methods section).Each lattice point in a random lattice is then assigned a local composition, which is defined as the composition of the atom and its neighbors up to the 5th nearest neighbor within two adjacent (110) planes.The 5th nearest neighbor was chosen based on atomistic simulations of solute-dislocation interaction energies [38], and only atoms within the two (110) planes that would be sheared by a dislocation were considered.A local USFE value was then determined using the USFE-composition maps in Fig. 3.In total, 420 lattices containing several million atoms were created for each level of SRO for both alloys, with examples shown in Fig. ??.To quantify SRO based on the Warren-Cowley parameters, we introduce ω as the Euclidean norm of the three like-pair Warren-Cowley parameters ω = α 2 11 + α 2 22 + α 2 33 as a single figure of merit (FOM) ω that can be used to compare the relative extent of SRO in an equimolar alloy.With it, we find that the mean USFE scales directly with ω (Fig. ??A).The MoNbTi alloy reaches a higher ω for the same set of annealing treatment due to the preference for Mo-Ti and Mo-Nb bonds.The coefficient of variation (CV) in USFE tends to decrease with ω, a consequence of the reduction in random atomic pairings (Fig. ??B).

Dislocation glide mechanisms
We first study the role of SRO on the glide behavior of initially straight edge or screw dislocations.Due to the randomness in underlying fault energies, twenty independent realizations are performed for each alloy and each level of SRO.To study critical behavior, the applied shear stress is gradually increased in increments of 0.001µ until the dislocation glides and is held constant until it fully arrests.
Fig. 4 shows snapshots in time of edge dislocation glide.When the stress is initially applied and raised, the dislocation remains straight.Once the applied stress exceeds the first threshold, the edge dislocation becomes slightly wavy, as small portions of the dislocation line bow out into low USFE regions and are held back at the higher USFE regions.The stress must be increased further for the dislocation to glide, and small bowed out segments of the dislocation will glide independently through the lower USFE regions, dragging the neighboring (non-edge) segments through the higher USFE regions.The dislocation arrests several times during the simulation, each time requiring the stress to be raised to restart glide.The arrested dislocation morphologies are wavy, unlike the original pure edge orientation.The stop/start behavior leads to glide plane hardening, a continual increase in applied stress with increasing plastic strain, as seen in the stress-strain curve in Fig. 4.This is in contrast to PFDD simulations of pure metals, in which dislocations gliding have a single critical stress and remain straight during glide.E-H) dislocation glide and their associated stress-strain curve.The edge example is in TaNbTi at 300 K SRO, and screw example is in MoNbTi at 300 K SRO.
Screw dislocation glide proceeds in a different manner from edge dislocation glide.As the stress increases, the dislocation remains completely straight with pure screw character until a kink-pair only a few Burgers vectors wide is nucleated into a low USFE region.Unlike the variable wavy bow out in the edge dislocation, these kink-pairs always have a height of just 1b, and the kinks will usually, but not always, glide along the length of the screw dislocation to advance the full dislocation line forward.Like the edge dislocations, the screw dislocations may become arrested under stress, but unlike the edge dislocations, the arrested dislocation morphologies are nearly straight, apart from a few metastable kinks, recovering the original pure screw orientation.The start/stop mechanism of glide of the screw dislocation also leads to glide plane hardening, although at a lower level than edge dislocation glide plane hardening.

Dislocation glide stresses
Fig. 5A plots the stresses to initiate glide σ i and the final stresses for runaway glide σ f of edge and screw dislocations based on twenty independent initializations.Regardless of SRO, both σ i and σ f for the TaNbTi alloy are lower than those for the MoNbTi alloy.From Fig. 5B, we find that mean glide resistance across the plane scales directly with the average USFE across the plane.Thus, changes in USFE caused by SRO have a direct influence on the stress to initiate and propagate dislocations, as seen in Fig. 5C.The higher the ω, the higher the USFE is increased relative to the RSS case, which translates directly to increased glide resistance for both screw and edge dislocations.
We also find that the hardening in glide resistance is related to the degree of dispersion of the USFE values in the glide plane as opposed to the mean.In Fig. 5D, we analyze the role of composition and its fluctuations by adopting the fractional increase from σ i to σ f as a measure of glide-plane hardening.Importantly, for both alloys and screw/edge character dislocations, hardening scales directly with the USFE CV.The strikingly linear relationship even when considering both alloys implies that it transcends composition.Thus, apparent differences in the hardening seen in these alloys can be explained.Compared to TaNbTi, MoNbTi achieves, on average, greater hardening in the ideal random case and lower hardening in the highest ω case.
Further, while hardening for both screw and edge dislocations increase as the USFE CV increases, the edge dislocations experience greater hardening than the screw dislocations for the same statistically sampled glide plane length (Fig. 5D).The edge dislocations glide by depinning of the segments at the relatively harder regions, segments which have reoriented to non-edge character due to bow out.Continued glide, therefore, relies on overcoming those local regions of higher resistance.Encountering a region ahead of the dislocation of even greater resistance than in the wake more likely occur when the dispersion in USFE is greater.The screw dislocation moves by producing short and narrow atomic advances of screw-oriented segments, i.e., kink-pairs, in the weaker regions and relying on the long advances of easier-to-move edge segments along the length of the dislocation.Those local domains of higher resistance that are more likely encountered when the dispersion in USFE is higher, can be easily overcome by migrating edge dislocations.By virtue of their differing glide mechanisms, edge dislocations experience greater sensitivity to the dispersion in USFE and hence greater hardening than screw dislocations.Fig. 5E examines the influence of ω on hardening.It reveals that the role of SRO on hardening corresponds to the extent to which SRO affects the CV in USFE.As increased ω tends to narrow the dispersion in USFE across the glide plane, it reduces glide-plane hardening.Since the TaNbTi alloy achieves lower ω than MoNbTi for the same annealing treatment, hardening, like its strength, is weakly affected by SRO compared to MoNbTi.

Effect of chemical fluctuations on glide mobility
Next, we study the effect of chemical fluctuations on screw/edge glide mobility by examining loop expansion on the (110) slip plane at constant stress.For given applied stress, thirty dislocation loop expansion simulations are conducted representing different locations in a given alloy and SRO.We find that for all cases, the anisotropy in screw/edge behavior reduces with increases applied stress.At low stresses, the difference in screw/edge behavior is large, causing the loop to expand into an oblong shape (Fig. 6A-D).The edge segments move continuously, constantly changing their wavy appearance.The screw dislocations advance very slowly, nucleating only a few kink-pairs at a time and recovering the nearly straight orientation with each advancement.As the applied stress increases, both screw and edge velocities increase and their ratio decreases towards unity.The loop expands more isotropically, almost FCC-like (Fig. 6E-H).Figure 6H shows the loop at both low and high stress overlaid together.The two loops have similar widths (103b and 101b, respectively) but different heights (38b and 51b, respectively), demonstrating that higher stresses decrease the loop aspect ratio.The screw dislocation has clearly changed its mode of glide, as a result of a higher kink-pair nucleation rate.The screw portions move continuously and adopt a wavy appearance, superficially much like the edge dislocations.Wavy screw glide, however, occurs as many kink-pairs nucleate simultaneously along the same dislocation.Newly advanced portions can nucleate further kink pairs, causing different portions of the dislocations to advance at different rates.
We separate the two extremes of screw dislocation behavior into "jerky" glide at low stresses and "smooth" glide at high stresses.To link the transition between jerky and smooth dislocation the rate of kink-pair nucleation, we calculate the waiting time between kink-pair nucleation events for all simulations.The waiting times from all thirty instantiations of loops are combined into a single distribution for a given alloy and applied stress, and the means are plotted in Fig. 7.Over 22,000 and 52,000 total waiting times were recorded for TaNbTi and MoNbTi, respectively, and each individual distribution contains at least 200 values.For both TaNbTi and MoNbTi, higher levels of SRO correspond to longer average waiting times and thus more jerky dislocation glide at all applied stresses.Normalizing the applied stresses by the average σ i for screw glide, the three distinct SRO curves collapse into one.Thus, the critical stress to transition from jerky to smooth scales with Peierls strength or with static strength.SRO affects the transition stress in dynamic glide indirectly via its strengthening effect on static glide resistance.

Discussion
There are limited experimental measurements of the mechanical properties of TaNbTi and MoNbTi.The tensile yield stresses have been reported as 620 and 950 MPa for TaNbTi and MoNbTi, respectively [8,42,43], which are consistent with the dislocation dynamics   6: The same dislocation loop expanding in MoNbTi under different applied stresses.The initial loop shape is shown by the dotted lines in (A) and (D).When a lower stress is applied, the screw dislocation nucleates kink-pairs infrequently, causing jerky dislocation glide and remaining largely pure screw.At higher applied stresses, the screw dislocation nucleates many kink-pairs at once resulting in smoother glide and a wavy morphology.The final loop from (C) is reproduced by the dashed line in (F) to highlight the difference in aspect ratio between the two loops.predictions of higher glide stresses for MoNbTi.The ultimate tensile strengths are 683 and 1500 MPa, respectively, so MoNbTi exhibits significant strain hardening while TaNbTi does not.Our simulations show higher hardening in MoNbTi than TaNbTi due to the increased coefficient of variation in USFE.The dislocation dynamics simulations also revealed that the amount of dislocation hardening is decreased by the presence of SRO, especially for MoNbTi.As edge dislocations undergo more hardening than screw dislocations, we predict the differences in macro-scale strain hardening in these alloys are largely controlled by edge dislocation behavior.
From the Warren-Cowley parameter calculations, it is clear that MoNbTi has a higher propensity for SRO and will be more affected by processing conditions.There has been interest in tuning the SRO parameters through heat treatment, although experimentally this is difficult to achieve [14].Findings indicate that MoNbTi is a better candidate for exploring SRO strengthening than TaNbTi since SRO promotes two relatively stronger Mo-Nb and Mo-Ti bonds.However, the relative increases in the average dislocation glide resistance due to SRO amount to less than 20% even in the most extreme cases, so large changes in the mechanical properties must be accompanied by changes in the chemical composition, not SRO alone.
Via MD simulations, a few studies have shown SRO-enhanced dislocation glide resistance [14,27,44] and one study on CoFeNiTi alloy reported a slight SRO softening [45].Strengthening or softening was related to the formation of immobile dislocation segments via cross slip.Here, we demonstrate SRO strengthening in dislocation glide without cross slip.While the current dislocation model permits cross slip [46], it was not observed in the present calculations since the influence of thermal fluctuations is not taken into account.Including temperature would undoubtedly increase the chance for cross-slip or cross-kinking, adding another mechanism for SRO strengthening.

Machine-learning interatomic potential
MTP, which has been shown to yield an excellent balance between accuracy and computational cost [30], was utilized in this work.Details about the MTP formalism can be found in previous works for interested readers [47][48][49].Briefly, in MTP model, the energy E is partitioned into the contributions (V ) of individual atomic neighborhood (n i ).V (n) is linearly expanded through a set of basis functions B α : where the basis functions B α (n) are given as: where r ij is the position of the neighbor jth atom relative to the ith atom, and f µ are the radial functions, which only depend on the interatomic distance |r ij | and atomic types z i and z j .The terms r ij ⊗ • • • ⊗ r ij include angular information about the neighborhood n i and are tensors of rank ν.The summation is carried out over all the atoms in the neighborhood n j around ith atom.M µ,ν , i.e., the moment tensors, are invariant with respect to translation of the system and permutation of equivalent atoms.The cutoff distance R cut and maximum level lev max are the two important hyper-parameters that determine the accuracy and computational cost of MTP.Their values in this work are set to 4.8 Å and 20, respectively.The loss function used for optimizing the MTP is given as: where θ are the set of parameters to be optimized; x j are the configurations in the training data; E(θ, x j ) is the MTP predicted energy; E DF T (x (j) ) is DFT energy; f (θ, x j ) is the MTP predicted force; and f qm (x (j) ) is DFT force.w e , w f are non-negative weights expressing the importance of energies, forces in the optimization problem.

Training data generation
The comprehensive training data are essential to the robust potential development.The training data include the elemental, binary, ternary, and quaternary chemical systems.From the structure type perspective, ground state bulk structures, distorted bulk structure, surfaces, stacking faults, and grain boundaries are also considered to improve the accuracy for potential in calculating planar defects.The detailed structure generation for each system is provided as follows.
1. Elemental systems (BCC Mo, Ta, Nb, Ti and HCP Ti) (a) Ground state structures for the elements, i.e., BCC Mo, Ta, Nb, and HCP Ti, plus BCC Ti; (b) Distorted structures constructed by applying strains of −10% to 10% at 1% intervals to the bulk conventional cell for BCC Mo, Ta, Nb, Ti and HCP Ti.Detailed can be found in ref. [50]; (c) Surface structures of elemental metals obtained from the Crystalium database [51]; (d) Grain boundary structures of elemental metals obtained from grain boundary database [52]; (e) Snapshots from N V T ab initio molecular dynamics (AIMD) simulations of the bulk 3 × 3 × 3 supercell for BCC and 4 × 4 × 2 supercell for HCP Ti at room temperature, medium temperature (below melting point), high temperature (above melting point).In addition, snapshots were also obtained from N V T AIMD simulations at room temperature at 90% and 110% of the equilibrium 0 K volume.Forty snapshots were extracted from each AIMD simulation at intervals of 0.1 ps; 100 snapshots from high-temperature (above melting point) N P T simulations are also included for each element.2. Binary systems (Mo-Ta, Mo-Nb, Mo-Ti, Ta-Nb, Ta-Ti, Nb-Ti) (a) Solid solution structures constructed by partial substitution of 2 × 2 × 2 bulk supercells of one element with the other element.Compositions of the form A x B 1 − x were generated with x ranging from 0 at% to 100 at% at intervals of 6.25 at%.
With each alloying percentage, we performed structure relaxations for all symmetrically distinct binary solid solution structures.We include both unrelaxed and relaxed structures in the data set.For ternary and quaternary systems, we considered the permutation of the elements in SQS structures.We also included structures from each ionic relaxation step of permutated SQS structure samples in the data set.
All the training data were generated using the same convergence criteria used in our previous work [26].The database includes 17210 static calculations as training data and 1660 as test dataset.

Atomistic simulations
All the atomistic simulations using MTP were performed using LAMMPS [56].The simulation cells are oriented along the [112], [ 110] and [11 1] cubic crystollagraphic directions, which are the main directions of interest for slip in the bcc system.A supercell of dimensions of ∼ 48 Å× 46 Å× 46 Å with 5760 atoms was used.
• Hybrid MC/MD simulations within the isothermal-isobaric (NPT) ensemble were carried out for different compositions (equimolar, element ratio of 3:3:4 and 1:1:3) for ternary alloys Mo x Nb y Ti z and Ta x Nb y Ti z .The supercells were randomly initialized based on the desired stoichiometry, i.e., a RSS.The simulations were then carried out at 300 K and 1673 K to obtain the equilibrium configurations with different degrees of short-range order.The Metropolis algorithm was used to attempt a swap between species type 1 and 2, species type 1 and type 3, and species type 2 and type 3 at every time step.The MD time step was set to 2 fs.The evolutions of energy, temperature, and SROs of different pairs for various compositions from MC/MD calculations can be found in SI Fig. ?? -??.After the converged configurations were reached through MC/MD, all the systems were then relaxed at 300 K using NPT MD.The structures with different composition and SROs were used for stacking fault energy calculations with energy minimization at 0 K.
• USFE calculations.The conjugate gradient algorithm was used to perform stacking fault structure relaxation and energy minimization.The force components within (110) plane are set to zero.The direction perpendicular to (110) plane was allowed to fully relax.The schematic in Fig. 8 shows the procedures of sampling the stacking

Creation of Random Lattices
Our method to generate the random lattices with SRO is based on the method of Gehlen and Cohen, who used a swapping procedure to create binary FCC lattices with SRO [41].We extend this method to multi-component alloys and only consider SRO in the first nearest neighbor (NN) shell.First, a BCC lattice of the desired shape is created, and each lattice point is randomly assigned an element type such that the overall composition is equilimolar.The total number of each bond type (e.g.Nb-Nb, Nb-Ti, etc.) is calculated.Using the Warren-Cowley SRO parameters, a goal number of each type of bond can be determined.Atoms are then swapped until the actual number of each bond type is equal to the goal.At each iteration, two lattice points with unlike element types are randomly selected.If swapping these atoms will move the bond numbers closer to the goal, the atoms are swapped.If not, the swap is rejected and a new pair of lattice points is selected.To increase the likelihood that a swap is accepted, the atoms selected for a swap are not chosen from the entire lattice but instead from a subset of atoms with desirable neighborhoods.For each element type, the expected number of each type of element in its NN shell can be predicted.When randomly selecting atoms to swap, we only choose from lattice points that do not have the expected atoms in their NN shell.For example, if the current number of i-j bonds is higher than the goal number of i-j bonds, we can select an i-atom with more j-atoms than expected in its NN shell to swap with a j-atom with more i-atoms than expected in its NN shell.Using this smarter swapping method for an MoNbTi lattice with the 300 K SRO parameters increases the acceptance percentage from 15% to 63% and decreases the running time by more than a factor of 10.
Once the lattice has the desired SRO parameters, we assign each point a composition based on the type of the atom and up to its 5th nearest neighbor within two adjacent (110) planes for a total of 31 atoms.This composition is identified on the linearly interopolated USFE plots in Fig. 3 in order to assign each lattice point a USFE value.

Phase-Field Dislocation Dynamics
PFDD is a mesoscale model used to study dislocation glide through crystals.The full details of the PFDD method are given elsewhere [34,39,57].PFDD tracks the dislocation configuration through scalar order parameters φ α (r) where α corresponds to a slip system with Burgers vector b α and plane normal n α .When φ α (r) = 0 or 1, point r is unslipped or slipped by a full dislocation of type α, respectively.Intermediate values of φ α (r) at the interface between slipped and unslipped regions correspond to dislocations.The order parameters are used to calculate a total system energy ψ(φ(r)), which is then minimized by the Ginzburg-Landau equation to evolve the dislocations.
where t is time and m 0 is a dislocation mobility coefficient.There are three energy contributions to the total energy: the lattice energy ψ latt (φ(r)), the elastic energy ψ elas (φ(r)), and the external energy ψ ext (φ(r)): The elastic energy accounts for the elastic strain in the crystal and is given by where c ijkl is the elastic stiffness tensor and e ij is the elastic strain which can be calculated from the order parameters [58].The external energy accounts for the work done by an applied stress field and is given by where σ app ij is the applied stress state.The final energy term, lattice energy represents the energy to break bonds at the dislocation core, and its form is specific to the material being studied.Because dislocations in BCC crystals have compact cores, a simple sinusoidal approximation is used [59]        The coefficient of variation in USFE for alloys with SRO normalized by the coefficient of variation in USFE for the no SRO case, plotted against the SRO figure of merit ω.In both plots, the black × represents the RSS case.

Figure 4 :
Figure 4: Representative examples of edge (A-D) and screw (E-H) dislocation glide and their associated stress-strain curve.The edge example is in TaNbTi at 300 K SRO, and screw example is in MoNbTi at 300 K SRO.

Figure 5 :
Figure 5: Critical stresses for dislocation glide.(A) The mean critical stress for dislocations in each of the alloys.The critical stresses for MoNbTi and TaNbTi are normalized by their respective shear moduli µ.The error bars show the resolution of the simulation, 0.001 µ. (B) The mean initial critical stress vs the mean USFE for each alloy.(C) The mean initial critical stress relative to the mean initial critical stress for RSS alloy vs the extent of SRO, represented by the SRO FOM.The black × represents the RSS case.(D) The hardening, represented by the difference in final and initial critical stresses, vs the coefficient of variation of USFE.(E) The hardening relative to the hardening for the RSS alloy vs the extent of SRO.

Figure
Figure6: The same dislocation loop expanding in MoNbTi under different applied stresses.The initial loop shape is shown by the dotted lines in (A) and (D).When a lower stress is applied, the screw dislocation nucleates kink-pairs infrequently, causing jerky dislocation glide and remaining largely pure screw.At higher applied stresses, the screw dislocation nucleates many kink-pairs at once resulting in smoother glide and a wavy morphology.The final loop from (C) is reproduced by the dashed line in (F) to highlight the difference in aspect ratio between the two loops.

Figure 7 :
Figure 7: Average time between kink-pair nucleation events as a function of applied stress.The error bars show the standard error.The transition between jerky and smooth screw dislocation glide is determined by the applied stress relative to the dislocation critical stress.

Figure 8 :
Figure 8: The procedures of sampling the stacking faults with different local compositions.

)Figure S1 :
Figure S1: Comparison of elastic constants of elemental systems calculated using density functional theory (DFT) and the trained moment tensor potential (MTP) for the Mo-Ta-Nb-Ti system.The values computed by MTP reproduce the DFT values for all elements.

Figure S4 :
Figure S4: The evolutions of potential energy, temperature, and SROs (α ij ) from 300K and 1673K for Mo-Nb-Ti systems with a ratio of 3:4:4.With the cell size of 48 Å× 46 Å× 45 Å of 5760 atoms, the energy and SRO converge after 200 ps for all cases.

Figure S5 :
Figure S5:The evolutions of potential energy, temperature, and SROs (α ij ) from 300K and 1673K for Ta-Nb-Ti systems with a ratio of 3:1:1.With the cell size of 48 Å× 46 Å× 45 Å of 5760 atoms, the energy and SRO converge after 200 ps for most cases except for Ta 3 NbTi and TaNbTi 3 , for these two, energy and SRO converged after 800 ps.There is a significant bonding preference for the Ti-Ti pair and Nb-Ta pairs in Ta 3 NbTi and TaNbTi 3 .The energy drops due to the formation of SRO is within 10 meV/atom in Ta 3 NbTi and 10 meV/atom in TaNbTi 3 .

Figure S6 :Figure S7 :
Figure S6:The evolutions of potential energy, temperature, and SROs (α ij ) from 300K and 1673K for Mo-Nb-Ti systems with a ratio of 3:1:1.With the cell size of 48 Å× 46 Å× 45 Å of 5760 atoms, the energy and SRO converge after 200 ps for all cases.The non-equimolar composition significantly affects the preferred bonding pairs reflected in SROs.I.e., the atoms with a higher percentage prefer to surround and bond with those with a lower percentage.For example, in Mo 3 NbTi, Nb-Mo pair and Ti-Mo are the most preferred bonding pairs.In MoNb 3 Ti, Mo-Nb and Ti-Nb pairs are the preferred bonding pairs.Mo-Nb is slightly more preferred than the Ti-Nb pair.The energy drops due to the formation of preferred bonding are within 10 meV/atom for Mo 3 NbTi and MoNb 3 Ti and within 20 meV/atom for MoNbTi 3 .

Figure S8 :
Figure S8: Example (110) planes with colored by local unstable stacking fault energy (USFE) values for both alloys at different levels of short range order (SRO).The mean and standard deviation of the normally distributed USFEs are given in mJ/m 2 .

Figure
Figure S9: (A) The mean USFE for alloys with SRO normalized by the mean USFE for the random solid solution (RSS), plotted against the SRO figure of merit ω. (B)The coefficient of variation in USFE for alloys with SRO normalized by the coefficient of variation in USFE for the no SRO case, plotted against the SRO figure of merit ω.In both plots, the black × represents the RSS case.