High entropy alloy strengthening modelling

High entropy alloys (HEAs) have recently drawn attention due to their excellent mechanical properties across wide temperature ranges. This is attributed to phase stability and a wide variety of strengthening mechanisms in operation. Solid solution, precipitation, dislocation, grain-boundary, twin-boundary and phase-transformation strengthening have been reported to play an important role in controlling their mechanical properties. With a focus on yield strength, this paper reviews the different hardening mechanisms reported in the literature. Mathematical formulations and key constant for describing each mechanism are presented and discussed. A strengthening mechanism modelling strategy for HEA design is outlined.


Introduction
From ancient times, the development of metallic materials had great impact on human development. The first industrial revolution incorporated the metallic alloys manufacturing technology that shaped today's industry and enhanced our living standards. A key family of alloys is superalloys, which development began in the 1930s; they demand high-temperature phase stability and excellent mechanical properties, as well as good corrosion resistance. More than 30 superalloy families have been developed until now; each superalloy family is dominated by a principal metallic element, and different alloys per family are generated by elemental additions. Traditional alloys display a principal metallic element and their properties are tuned through processing. HEAs are instead composed by various major principal elements in high concentrations, typically not exceeding 35% [1,2]. Disorder caused by multiple elements tends to inhibit the formation of brittle phases by the random scattering of atoms of each element, leading to improved strength and toughness. So far, several HEAs with outstanding properties beyond those of conventional alloys have been discovered, and new superior HEAs are still expected to be developed in future. HEAs are generally defined as having more than five major elements [1]. The atomic concentration per major element is between 5 and 35%, with minor elements being under 5%. HEAs are expected to have a high mixing entropy at high temperatures. Following the cognition of traditional physical metallurgy and multi-phase diagrams, alloys with multiple elements will produce many phases and intermetallics, which make their microstructure quite complex. The high entropy effect causes the elements to mix into one or more simple solid solution phases. The phase change from high to low temperature is relatively simple. HEAs bring a new concept to alloy design. As 5 to 13 miscible elements can be incorporated in equimolar ratio, more than 7000 equiatomic alloy systems can be obtained [1], which much exceeds the 30 families of traditional superalloys. It is hard and inefficient to design high performance HEAs by traditional trial and error approaches, so the simulation of their strengthening mechanisms is an attractive approach to their conception.
A key advantage of HEAs lies in their combination of properties. As reviewed by Zhang et al [2], there is an excellent improvement in HEAs compared to regular alloy systems. HEAs outperform conventional single-principal element alloys in terms of high temperature strength, fatigue, corrosion resistance, ductility, and irradiation-resistance properties [2]. The hardening mechanisms of HEAs include solid solution, precipitation, grain refinement, and Hall-Petchtype effects, as shown in figure 1 [3]. For example, carbon-containing CoCrFeMnNi displays a good combination of strength and ductility after annealing at 800 • C, which is attributed to solid-solution strengthening by carbon and strong grain boundary strengthening [4]. Dynamic recrystallisation, which is a valid way to control grain size, was studied in CoCrFeMnNi alloy [5].
Compositional design is an effective way to alter the mechanical properties of HEAs. Many relevant studies have been reported. He et al investigated the influence of Al on the CoCr-FeMnNi HEA system. Both fracture and yield strength significantly increased at the expense of plasticity [8]. Al in AlNbTiVZr 0.5 exhibited second-phase and solid-solution strengthening [9][10][11]. By changing the Al content in Al x NbTiVZr (x = 0, 0.5, 1, 1.5), the Laves phase and Zr 2 Al can be controlled. Similarly, Al changes to Al x CoCrFeNiTi resulted in a phase transformation from fcc to bcc structure with sub-grains and nanosized precipitates formation [12]. It is well documented that Al addition stabilises the bcc structure in the Al x CoCrFeMnNi [8] and Al x CoCrCuFeNi [1], transforming the crystal structure from fcc to bcc. Al can be added to refractory HEAs to reduce density, Senkov studied the effect of Al on HfNbTaTiZr and CrMo 0.5 NbTa 0.5 TiZr alloy systems [13]. The yield strength and hardness at room temperature increased significantly at wide temperature ranges with Al replacing Cr and Hf. In this case, there is no phase transformation with Al addition to HfNbTaTiZr and CrMo 0.5 NbTa 0.5 TiZr alloy systems. Here, Senkov [13] suggested that the strengthening mechanism of Al x Hf 1−x NbTaTiZr is related to the formation of stronger interatomic bonds due to the strong bonds of Al with each alloying element.
Stepanov et al investigated the AlCr x NbTiV (x = 0, 0.5, 1, 1.5) HEA system [14]. Variations in Cr content significantly affect properties and microstructures. Strength is increased with an increase in Cr content, but at the expense of room temperature ductility. But ductility increased at high temperatures (600, 800, and 1000 • C). Other studies show that the addition of Cr prevents environmental cracking [15]. Hardening mechanisms of HEAs: (a) solid-solution hardening: the lattice distortion increases the resistance to dislocation movement [3]; (b) grain-boundary strengthening: the slip resistance is increased because the grain boundaries block dislocations [3]; (c) phase-transformation strengthening: the hard phases increase strength while the soft phases provide the deformation ability. Reproduced from [3]. CC BY 4.0.; (d) dislocation strengthening: dislocation forest pins or obstructs another group of dislocations. Reprinted from [6], Copyright (2020), with permission from Elsevier.; (e) precipitation strengthening: the process of cutting through small size precipitates; or large size precipitates that block and bend moving dislocations; (f) twin-boundary strengthening: twin boundaries diminish dislocation motion. Reproduced with permission from [7]. '© Cambridge University Press. Reprinted with permission.' [7].
Here, alloying effects are interpreted in terms of their influence on mechanical properties. The addition of Al [16], C [17], and Cr [15] would result in martensitic transformation, nanoprecipitate formation, stacking fault energy (SFE) changes, and increased twining and ductility of the alloy. Elemental concentration variations have a great impact on the microstructure of HEAs.
A comprehensive study was conducted in relating HEAs with their mechanical properties in [2]. The severe lattice distortion and 'cocktail effect' makes glide difficult due to a wide range of lattice spacings. Microstructural features such as disordered crystal lattices, chemical disorder, stacking faults, and twins contribute to the rich behaviour of HEA systems, particularly at high temperatures. Inhomogeneous lattice distortions can lead to nucleation and retard phase growth. Stepanov et al [14] show that an increase in Cr content promotes Laves phases formation which precipitation contributes to strength increase. The alloy was homogenized at 1200 • C for 24 h. 13% and 35% of hexagonal C14 Laves phase are found in the bcc matrix of AlCrNbTiV and AlCr 1.5 NbTiV alloys, respectively. CoCrFeMnNi is often associated with twinning [18]. As strain is applied, the number of primary and secondary twins increases and the mechanical properties improve significantly. Precipitation, twinning, and the matrix solute concentration are key mircostructural features for strengthening HEAs.
HEAs display excellent mechanical properties stemming from various strengthening mechanisms in operation, sometimes simultaneously. Due to the nearly infinite combination of elements and concentrations, it is impractical to find optimal HEAs by trial and error. Therefore, further analysis of their strengthening mechanisms and their quantification through modelling may be of great significance in guiding the design of novel alloy compositions. Many efforts have been carried out to describe the strengthening mechanism of HEAs. In the following sections, we discuss the modelling of HEAs according to different strengthening mechanisms, including solid-solution strengthening (Δσ ss ), precipitation strengthening (Δσ ppt ), dislocation strengthening (Δσ dh ), gran-boundary strengthening (Δσ gb ), twin-boundary strengthening (Δσ twin ), and phase-transformation strengthening (Δσ trip ). The overall yield strength of a HEA can be calculated in an additive manner as follows: σ y = σ fr + Δσ ss + Δσ ppt + Δσ dh + Δσ gb + Δσ twin + Δσ trip (1) where, σ fr is the lattice friction stress [19].

Solid-solution strengthening (Δσ ss )
One of the distinctive features of HEAs is their solid-solution strengthening effect being significantly higher than that of conventional alloys. It is widely accepted that solid solution hardening (SSH), especially at high temperatures, is the main cause for the extraordinary mechanical properties of HEAs. SSH mechanism is generally described in terms of atoms randomly dispersed within a solvent; the solute atoms interaction with dislocations will be augmented by the misfit in their atomic sizes and elastic moduli. Stationary solute atoms acting on moving dislocations lead to frictional models, which are divided into dilute and concentrated alloys, based on the solute spacing, the solute-dislocation interaction energy, and the range of interaction. Sriharitha et al [20] showed that half of the yield strength of HEAs stems from SSH and order strengthening, which occurs when the precipitate is in an ordered crystal structure which bond energy is altered as a result of shearing. Compared to conventional alloys, there are a few computational models for HEAs strengthening, but they typically neglect their complexity. SSH models, initially proposed by Fleisher for binary systems [21], assumed a low concentration of solute atoms. The SSH effect Δσ ss was expressed as: where X i is the i solute content and B i is a constant dependent on the mismatch parameter i , the shear modulus G, and a fitting parameter Z: Here, η i is the elastic misfit, δ i is the atomic size mismatch, and α i is a parameter relating the interaction forces between dislocations and solute atoms. A list of these parameters for HEAs is shown in table 1.
Labusch [22] further considered high solute concentrations. Due to the constant interaction of with solute atoms, a gliding dislocation is subjected to frictional rather than a blocking effects. The expression is similar to Fleisher's model:  [43] (continued on next page) where However, neither of these two models considers the effect of temperature on SSH. Butt et al [23][24][25][26][27] studied the change in critical resolved shear stress (CRSS) with composition at low temperatures. The CRSS (τ ) is expressed as: where τ 0 is the CRSS at 0 K, m is a constant, k b is the Boltzmann constant and W 0 is a constant that describes the binding energy of an edge dislocation with the nearby solute atoms [28]. But equation (6) validity is limited as yield strength at high temperature ranges is observed to reach a plateau. Labusch [29] gave an explanation for this phenomenon. A dislocation segment mostly jumps forward to cross obstacles at low temperatures. At high temperature, a large amount of backward jumps take place due to thermal activation. A similar amount of forward and backward jumps results in a τ plateau as a function of temperature as a result of this thermally-activated mechanism. Some calculated lattice parameters of reported HEAs are listed in table 1, including average hardness (H ave ), average atomic size (r ave ), elastic mismatch inside the lattice (η), and atomic size misfit (δ) required to calculate Δσ ss . Gypen et al [69] proposed a methodology to calculate SSH in multicomponent alloys. The main assumption was that solute atoms in close proximity to each other do not interact. Therefore, each hardening parameter B i can be computed as if the solute i were in a binary system with the solvent. Here, the effect of each triplet or pair is treated as a whole atom. The method sets a principal element as predominant. It is consistent for multicomponent alloys, which is expressed as: However, Gypen's model does not consider the varying atomic size misfit generated by a continuously deformed crystal lattice with different solute species, as in HEAs. By adapting the Gypen's approach, Toda-Caraballo and Rivera-Díaz-del-Castillo (TR) [70] firstly proposed a model for calculating the SSH in HEAs, which extended the Vegard's law to multicomponent alloys. They calculated the interatomic spacing between solute i and j by performing algebraic manipulations, as expressed by where Δ s is the mean interatomic spacing change with composition. In this method, there is no reference solvent or solute atom, making it suitable for calculating the unit cell parameter of HEAs. The model divides the interatomic space of one varying element and an equimolar concentration element system. The s value can be obtained by varying one element content and applying equation (8) for each sub-binary system. As shown in figure 2, the different lines represent the interatomic distance value with the varying elements. The large orange bullet is the reference alloy CrFeNbV. The general elastic misfit δs/δX i thereby can be expressed as:  with this, the parameter for B Gi can be calculated, and then the SSH is computed from equation (7).
Based on above work, Toda-Caraballo [71] presented a general formulation for SSH effects in multicomponent alloys. For Toda-Caraballo model, each element can be treated as solute, while other elements are the solvent, similar to a pseudo-binary system. Compared to previous work [70], this approach used the matrix to describe the lattice distortion and then extended it to multicomponent alloys. The matrix allows the calculation of not only the lattice parameter and the variation with composition, but also of other parameters concerning the crystal stability in HEAs. In this work, the SSH effect is accounted based on the lattice distortion due to the presence of different elements rather than the solute-solvent interaction. Toda-Caraballo used the variation of elemental concentration to calculate the change in the lattice parameter. The minimum variation for the system is by changing one single atom. By replacing each atom type by another, there are n 2 possible 'remove-replace' cases. The variation of unit cell parameter with composition is expressed by matrix D: here, a is the unit cell parameter. The SSH effect in multicomponent systems can be expressed as: Although this model describes well the SSH of different solutes and crystal structures, its accuracy still needs to be improved. Its application is shown in figure 3, the model predicts the magnitude of SSH effect and the trends with different elemental contents for a range of alloys. The SSH and grain size (Hall-Petch) effects are the main factors affecting yield strength. However, the grain size is poorly reported in most references, so the SSH model cannot provide an accurate prediction of yield strength.
Wang and Xu [72] applied the classical solid-solution theory to research (TiZrNbTa) 100−x Mo x (0 x 20) series alloys which aims to clarify the influence of Mo on microstructure and mechanical properties. The authors developed Labusch model to calculate the SSH effect. They first simplified the lattice as a fictive lattice, ignoring the overall mismatch in size and modulus. Then, they integrated all the distortion collectively and applied it to self-bonding configuration. The SSH is expressed as: where Z is a constant depending on the fictive solvent and temperature. By investigating 29 single phase refractory HEAs in the literature, Z is statistically determined to be 0.0074 ± 0.0011. The predicted σ y was compared to experiment in 34 single-phase refractory HEAs. The results show good agreement with experiment with a deviation of 27%. The predictive values produced by this model are more accurate compared to Toda-Caraballo's model. A critical factor in SSH modelling is the solute-dislocation interaction energy. Besides being able to compute this parameter analytically using the elasticity theory of dislocations, it can also be determined from first principles, taking into account both elastic and chemical contributions. Leyson et al [73][74][75] have used the solute-dislocation interaction energy calculated from first principles to explain the dependence of the CRSS on solute concentration and temperature for certain dilute Al and Mg alloys. Due to difficulty in calculating the solute-dislocation interaction energy, this method cannot easily be applied to HEAs; it is hard to define solute and solvent atoms in an equiatomic HEA.
Recently, Varvenne et al [76,77] proposed a parameter-free predictive theory for the composition-, temperature-, and strain-rate-dependence of the yield strength for fcc HEAs. Their work considers the elemental component as a solute embedded in an effective matrix of the surrounding alloy. A reference material is defined for the HEA displaying the average properties of the alloy. Each individual solute then interacts with the dislocation in the reference matrix, thus accounting for the solute interactions with the average chemical surroundings. Strengthening is mainly achieved by energy and local elemental fluctuations. Each solute located at any position has an interaction energy along the dislocation line direction. The mechanical interaction energy U(x i , y j , z k ) is between the dislocation pressure field p(x i , y j ) and the mismatch volume ΔV of the solute to the matrix material. The long straight dislocation usually changed to a wavy configuration because part of it is attracted to energetically-favourable fluctuations and repelled by unfavourable counterparts, as shown in figure 4. With the minimum energy configuration, each segment ζ is located in a local potential energy well. It is necessary for the dislocation to overcome ΔE b by thermal activation in order to glide. An asymptotic approximation is adopted to describe how the stress-dependent energy barrier behaves: In order to rationalise solid-solution strengthening in CoCrFeMnNi HEA, Okamoto [78] proposed mean-square atomic displacements (MSADs) based on first-principles total-energy calculations and successfully predicted trends in the strength of other solid solution alloys by determining the appropriate MSAD values, as the resistance to dislocation motion is solely determined by the amount of distortion in the crystal. The MSAD varies remarkably with the choice of elements in quaternary alloys, indicating that atomic displacement in equiatomic alloys is affected not only by the radius of an individual atom but also by the nature of the constituent atoms. Coury et al [79,80] improved the models by TR [70,81] and Varvenne et al [76,77] by adopting a more accurate atomic radii set based on the effective atomic radii for strength (EARS) method [82]. The crystal structure type and choice of elements in solution influence the effective atomic radius of pure elements, thus affecting the strength model predictions based on lattice distortion. Rather than putting atomic radii directly into the TR or Varvenne models, they used strength values from the literature to select the atomic radii that best match the experimental data. Based on this methodology, the Co 27.5 Cr 45 Ni 27.5 alloy was designed and determined to be 50% stronger than the strongest alloy previously reported for this system, that is, Co 33.3 Cr 33.3 Ni 33.3 [82].
Coury proposed an experimental methodology for rapid yield strength estimations of single-phase HEAs [80], which involves the TR-EARS methodology and the analysis of a compositionally-graded sample from a multiple diffusion joint, as shown in figure 5. To explore the vast range of compositional possibilities of HEAs, it may be more efficient to produce, characterise, and test compositionally-graded samples [83,84]. They determined the yield strength by analysing nanohardness results. By assuming full plasticity and an expanding cavity finiteelement-based model, Clausner and Richter [82] derived the following equation for converting nanohardness and elastic modulus into a single yield strength value: where H is the nanohardness, E is the Young's modulus, α is the effective cone angle of the Berkovich indenter, and K is a constant. This equation describes a linear relationship between hardness and yield strength for alloys with the same elastic modulus. In EARS approach, the elastic constants for all the HEAs are calculated by extrapolating from multicomponent alloys measured by ultrasonic methods for pure Cr, Mn, Fe, and Co [85]. Developing highthroughput methods to estimate and measure the mechanical properties of HEAs is crucial for the development of alloys with enhanced properties. In [85], the combination of TR-EARS and nanoidentation method proved to be a high-throughput alloy developmental tool. Using mechanical testing results combined with literature data, Coury et al [79] also developed a solid-solution strengthening model for both components that take into account the characteristics of single-phase bcc materials. They divided the yield stress of single phase bcc refractory HEAs into thermally-activated and athermal components. The athermal component is affected by elastic modulus mismatches and atomic size mismatches, which are generated by calculating the averages from each alloy. This allows for estimating the athermal component at high throughput. The yield stress component of thermally-activated refractory metals does not correlate with any physical parameter that can be calculated from averaged elemental atomic properties. Furthermore, it is greater than the values found for pure and diluted bcc refractory metals. For the thermally activated yield strength model, they used an 'activation energy barrier' which scales with a total energy F for a dislocation segment to slip where τ is the applied shear stress, τ a is the athermal critical resolved stress component, and τ 0 is the shear stress required to move a dislocation at 0 K. The activation energy can be related to dislocation motion through an Arrhenius expression and the thermal yield stress is derived as: where M is the Taylor factor. The total activation barrier (F) is expected to change with temperature, because the lattice parameter increases and the shear modulus softens. The full yield stress of selected refractory high entropy alloys is described, as shown in figure 6.

Precipitation strengthening (Δσ ppt )
Precipitation strengthening makes use of second-phase formation in the matrix whilst increasing temperature. A two-step heat treatment process incorporates solution-treatment at a high single-phase temperature, and further ageing at a low temperature to precipitate out secondphase particles. It is the blocking of dislocation motion by precipitates what induces the strengthening; this is influenced by precipitate particle size and lattice coherence with the matrix material. The particles maybe either cut or bypassed by dislocations. Precipitation strengthening thus takes contributions from shearing (Δσ sh ), where dislocations cut through small-sized precipitates, or from Orowan bypassing (Δσ Orowan ), where dislocations will pass by large-sized precipitates, bowing their dislocation lines, and eventually producing loops around precipitates [86]. As for the shearing hardening mechanism, dislocation motion is resisted by precipitates which are small and coherent with the matrix. Particle shearing is influenced by coherency strain, stacking-fault energy, interfacial energy, morphology, matrix and precipitate moduli, and lattice friction.
For small precipitates, there are three contributions of precipitation to hardening: particlematrix coherency (Δσ cs ), modulus mismatch (Δσ ms ), and atomic ordering (Δσ os ) [87]. The former two (Δσ cs + Δσ ms ) build up prior to shear, while the latter (Δσ os ) makes a direct contribution during the dislocation cutting through a particle. The equations associated to these mechanisms are [87,88]: where, α ε m is a crystal structure-dependent factor; ε m = 2/3 · (Δa/a) is the constrained lattice parameter mismatch, while Δa is the difference in lattice constant between the precipitate and the matrix and a is the lattice parameter of the matrix; r is the average particle radius; f is the precipitate volume fraction; b is the magnitude of the Burgers vector; ΔG is the shear modulus mismatch between the matrix and the precipitates and m 0 is a constant.
As for a dislocation cutting through particles, it is atomic ordering what dictates the strength increment, which is given by [87,88] Δσ os = M · 0.81 γ APB 2b 3π f 8 1/2 (19) where, γ APB is the average value of the anti-phase boundary energy for precipitates. If just one type of particle size is found after processing, the higher value between Δσ cs + Δσ ms and Δσ os represents the critical stress for shearing. Based on equations (17)- (19), precipitation strengthening gains is calculated for each mechanism separately.
Gradually, the resulting strengthening effect varies as precipitates grow with prolonged ageing time. Figure 7 shows the precipitation hardening of Al 4 (CoCrFeNi) 94 Ti 2 alloy aged at different temperatures as a function of time [89]. The peak hardness increased firstly then decreased upon overageing time at all temperatures, as shown in figure 7(a). At the beginning, the particle-shearing mechanism dominates accounting for the nano-sized and coherent precipitate contribution. As mentioned before, strengthening from matrix-particle coherence and modulus mismatch occurs first, while the order strengthening follows during shearing. In fact, the ordered particle size will affect the shear mechanism, but equation (19) does not incorporate particle size. Therefore, it is necessary to improve order strengthening models. For instance, the dislocation pairs cut through the small particles when they are weakly coupled. The shear stress increment is given by [89]: where γ APB is the anti-phase boundary energy of the precipitates, d is the precipitate particle size, φ is the precipitate volume fraction, and Γ = Gb 2 /2 is the dislocation line tension. On the other hand, the dislocation pairs tend to be strongly coupled when shearing large precipitates, and the shear stress increment becomes [90]: where w is a constant describing the elastic repulsion between strongly coupled dislocation pairs, and it can be approximated to 1. The hardness increment (ΔH) can be calculated from the shear stress increment (Δτ ) via the equation: Equations (20)- (22) have been applied to Al 4 (CoCrFeNi) 94 Ti 2 , and resulting hardness increment as a function of precipitate size at different ageing temperatures are plotted in figure 7(b). The calculation predictions fit well with experimental measurements. It should be noted that the hardness increases firstly to a peak, which can be explained by the weak-coupling mechanism (ΔH ∝ d 1/2 in equation (20)). Then, it decreases with further growing precipitate size due to the strong-coupling mechanism (ΔH ∝ d −1/2 in equation (21)). Strengthening occurs at the intersection of both mechanisms. Heat treated alloys contain particles with more than one type of morphology. As an example, Al 4 (CoCrFeNi) 94 Ti 2 HEA processed with the process 1 (P1) and process 2 (P2) methods had distinct particle morphologies in two distinct regions [87]. Due to the different particle spacing, sizes, and distributions, the strength gained from each region needs to be calculated separately. Therefore, the overall strength can be determined from an integrated model: (23) where, f I and f II are the volume fraction of the different particle regions, Δσ I ppt and Δσ II ppt are the precipitation strengthening contributions from two particle morphologies. The strengthening contribution from different hardening mechanism in [87] is calculated and shown in figure 8.
As precipitate particles grow, the incoherence with the matrix does, and the particledislocation Orowan bypass mechanism prevails. Following the bypassing and bowing of the dislocation line between two adjacent particles, a loop of dislocation may be formed around each particle, adding strength to the material by exerting a back stress [86]. Some HEAs, including Al 6 (CoCrFeNi) 91 Ti 3 , Al 9 (CoCrFeNi) 88 Ti 3 , and Al 14 (CoCrFeNi) 85 Ti 1 , may generate larger precipitate L2 1 -(Ni, Co) 2 TiAl particles (150 to 350 nm radius) coexisting with nano-sized, coherent L1 2 -Ni 3 (Ti, Al) particles (4 to 20 nm radius); both their strengthening mechanisms work concomitantly [87]. The shearing mechanism dominates for the nano-sized coherent L1 2 -Ni 3 (Ti, Al) particles, whilst the Orowan bowing may occur for the larger L2 1 -(Ni, Co) 2 TiAl particles. The contribution to strengthening of both species can be calculated [87]. For large precipitates, the contribution of Orowan bowing is calculated as [87]: where υ is the Poisson's ratio,r = 2/3 · r represents the mean particle radius on the slip planes, while r represents the average particle radius, λ is the mean inter-particle spacing.

Dislocation strengthening (Δσ dh )
Deformation in metals depends on dislocation glide; this may occur on various slip planes which interact differently when in motion. It is common that one group of dislocations pins or obstructs another group. Since dislocations are mutually hindering each other's motion, a higher applied stress is needed to initiate plastic flow. A higher dislocation density is associated with a higher yield strength, which results in dislocation strengthening. The Taylor hardening model is used to describe dislocation strengthening [91,92]: where M = 3.06 (for bcc and fcc materials) is applied, α h is a material-specific correction factor, and ρ is the dislocation density. Equation (25) requires a number of parameters to be determined. For an isotropic and homogeneous material, G = E/2(1 + υ), where υ is the Poisson's ratio. The magnitude of Burgers vector can be calculated by |b| = a/ √ 2. If the alloy is a multiphases system, the lattice constant a can be taken as the average value [91].
Different methods can be used to determine the value of α h . For example, α h is estimated as 1 in [91] and 0.2 in [87,92]. George et al reported a method to determine the hardening factor α h for CoCrFeMnNi [93]. If the modulus-normalised increment of the stress, a plot of Taylor hardening can be deduced from equation (25), showing the linear dependency of the normalised work hardening (σ Max − σ y )/(MG) with respect to (bρ 1/2 ), where σ Max is the maximum applied stress, σ y is the yield stress. In this work, α h was determined as 0.4 ± 0.1, as shown it figure 9.
A common route for determining the dislocation density ρ dis is from x-ray diffraction (XRD) [87,91,94]: in the equation (26), ε is the lattice strain present in materials which can be usually determined by the Stokes-Wilson formula [95], d is the average grain size. A more common technique to simultaneously calculate d and ε from experiment is the Williamson-Hall method [96,97]. This method relies on the theory of XRD peak line broadening β, which lies at the intersection of three distinct broadening effects, the crystalline broadening β G , strain broadening β S , and instrumental broadening β 0 The constant K G is usually approximated as being 0.9, λ WV is the wavelength and θ is the Bragg angle at selected diffraction peaks which is usually set to 0.154 05, the wavelength of Cu Kα radiation. A diffraction peak's line broadening is usually measured as the width at half of its maximum intensity. With a focus on the micro strain ε only, equation (27) can be rewritten as: Parameter ε can be determined by the slope of the linear fit of the β cos θ − 4 sin θ. By this method, the dislocation densities of the two differently processed Al 4 (CoCrFeNi) 94 Ti 2 HEAs were determined [87]. As shown in figure 10, The sample processed by P1 (30% cold rolling, annealing at 1273 K for 2 h, then ageing at 1073 K for 18 h, and water quenching) has a value for microstrain ε 1 ≈ 0, therefore the dislocation density of the sample can be calculated by equation (26), which approaches zero. On the another hand, the sample processed by P2 (70% cold rolling, ageing at 923 K or 18 h, and water quenching) results in ε 2 = 0.102 and ρ 2 = 5.02 × 10 14 m 2 . The reason for this difference is that the number of dislocations varies according to the treatment process. The former is because annealing decreases dislocation content, whilst the latter is attributed to the dislocation density increase by heavy cold rolling.

Grain-boundary strengthening (Δσ gb )
In grain-boundary strengthening, grain boundaries serve as pinning points, preventing further dislocation propagation. Dislocations need extra energy to change direction since the adjacent grains have different orientations. In addition, dislocations are unable to move in a continuous slip plane because of their disordered grain boundary. As the atomic mismatch between grains creates a repulsive stress field to oppose dislocation motion, clusters of dislocations pile up when they are unable to move past a boundary. By decreasing grain size, the amount of possible pile up at the grain boundary increases, increasing the amount of applied stress, and the higher the yield stress [98]. This is described by the Hall-Petch equation: where, σ 0 is a constant describing the stress at which dislocations begin to glide, k s is a materialdependent strengthening factor, and d g is the average grain diameter. Term k s d −1/2 represents the yield strength increment with the decreasing grain size. The value of k s is critical for calculating the grain boundary strengthening effect. Based on a first approximation, the k s value can be estimated using a rule of alloy mixtures [20]. The Hall-Petch effect becomes inverted in metallic systems as grain sizes decrease below 100 nm (the exact value depends on the material). For example, k s = 0.11 MPa m 1/2 for microcrystalline Cu, while it decreased to 0.05 MPa m 1/2 for 100 nm crystallite size. In the calculation of spark plasma sintered (SPS) Al x CoCrCuFeNi for Co, k s is calculated by the formula [20]: where β s is a constant and γ is the surface energy [99]. The values of σ 0 and k s can be obtained by combining equation (31) and experimental data, solving the equation for σ y . Grain size and yield strength may be altered through heat treatments. σ 0 and k s may then be determined but only if other mechanisms such as twinning induced plasticity (TWIP), transformation induced plasticity (TRIP) and precipitation are avoided. A dataset gathering literature values for σ 0 and k s for HEAs is shown in table 2. It should be noted that different processing conditions would result in slightly different σ 0 and k s for the same alloy (i.e. Al 0.3 CoCrFeNi). Ultra-fine grained microstructures with mean grain size under 1 μm have an extra-hardening effect, which will affect the fitting quality of Hall-Petch slope. As shown in figure 11 [100], the yield strength values for CoCrNi and Ni-40Co alloy deviate from linearity for grain sizes under 0.3 μm for CoCrNi and 1.5 μm for Ni-40Co, which is expressed by the dotted line, while the plot of pure Ni shows a standard linear relationship. According to Yoshida [100], CoCrNi has a slow grain growth speed due to sluggish diffusion, which results in an ultra-fine microstructure. Therefore, the grain size order of magnitude for calculating HEA fitting parameter is quite important.
It is widely accepted that the values for σ 0 and k s of HEAs are higher than those of their pure metal counterparts, even with the same crystal structure. In multicomponent alloy systems, the crystal lattice is distorted due to atomic size mismatch. σ 0 increases by solid-solution strengthening due to severe atomic size and modulus mismatch [100][101][102]. The higher k s in HEAs  indicates harder slip transfer between grains and stronger grain boundaries, because the dislocation lines in HEAs are twisted by severe lattice distortion rather than the straight dislocation lines observed in pure metals [103].

Twin-boundary strengthening (Δσ twin )
Dislocation glide can be altered by twin boundaries. Deformation twins can increase both the strength and the strain hardening of alloys [110][111][112][113]. Twin boundaries can serve as a stable interface for strengthening because of their extremely low excess energy [114]. Early works including computer simulations and experiment verified that twin-boundary strengthening becomes dominant by reducing the twin-boundary spacing to the nanometre range [115,116]. By gradually introducing new interfaces during deformation, twins reduce the mean free path of dislocations, resulting in stronger structures. On the other hand, deformation twins containing a high density of sessile dislocations impede dislocation glide, causing increased strain hardening [117][118][119][120][121]. This is commonly called 'dynamic Hall-Petch effect'. The phenomenon of nano-twin strengthening is reported in many HEAs, such as CoCr-FeMnNi [93,122,123] and AlCoCrCuFeNi [91]. In some cases, the twin-boundary strengthening is more prominent than that of grain refinement [124]. Further strengthening can be achieved by nanotwins for a high density of stacking faults [125].
Yield strength and twin spacing can be expressed by a formulation similar to the Hall-Petch relation:  [124].
where f t is the volume fraction of grains containing twins, λ t is the average twin spacing and k t is the strengthening coefficient by twin spacing. Equation (33) is also known as dynamic Hall-Petch relation.
The parameters f t , k t , and λ t are necessary for the calculation of the contribution from twins. f t and λ t can be extracted from statistical analyses using scanning electron microscope, or transmission electron microscope micrographs containing twins [126,127]. The value accuracy of f t and λ t increases with increasing the number of micrographs. It should be noted that only the area size of the grains containing twins can be obtained by two-dimensional micrographs. Therefore, the f t is usually assumed to a first approximation [91,127], or converted by certain geometric relationship [128]. If the number of twin micrographs is small, the average f t and λ t can be calculated from limited number of points [91]. From k t value is approximated as the slope of a σ y /λ −1/2 plot. The intercept σ 1 represents all mechanisms besides twins that contribute to strength. In order to get different data pairs, the processing parameters are usually controlled to prepare samples with differing twin spacing. For example, Su et al [124] used this method to obtain k t for an interstitial carbon alloyed HEA (Fe 50 Mn 30 Co 10 Cr 10 C 0.5 ) as shown in figure 12 [124]. In this case, the varied twin spacing is obtained by changing the annealing temperature and time, as shown in the figures 12(a)-(c). Figure 12(d) shows the dependence of yield strength σ y with average twin spacing λ The strengthening coefficient of nano-twins k t (195 MPa μm 1/2 ) was derived by linear fitting of equation (34). A summary of the contributions of nano-twins, dislocations, and grain refinement to the yield strength is provided in figure 12(e). In analysing the yield strength of bimodal microstructures created by tempering, it is evident that nanostructures within the shear bands and nano-twins within the parent grains contributed primarily to this performance improvement.
To obtain the k t value, approximate methods may be utilised when curve fitting is difficult. For example, the grain size in the nanocrystalline regime can be considered in nano-size to approximate the classical Hall-Petch coefficient.
This method has been applied to calculate twinning strength for some HEAs [91,127]. The values of k t , f t and λ t obtained from the above method in different HEAs are listed in table 3.

Phase-transformation strengthening (Δσ trip )
TRIP effect has been applied to design HEAs to improve both strength and ductility as a novel strategy [129][130][131]. The TRIP effect makes use of mechanical stability of the parent phase. The metastable phase would partly or totally transform into martensite under the action of mechanical loading or thermal processing. The SFE is a key factor determining the TRIP effect [131]; it is mainly determined by deformation temperature and chemical composition. Li et al sought out HEAs with SFEs and FCC-HCP phase stability in Co 20 Cr 20 Fe 40−x Mn 20 Ni x whose deformation-driven phase transformation is promoted by employing density functional theory calculations [131]. On the other hand, tailoring chemical composition to control phase metastability in HEAs; for example, Co 10 Cr 10 Fe 80−x Mn x system, shows TRIP effect (fcc to hcp transformation) when Mn content decrease to 30 at.% [129].
In fact, there are mainly two types of phase transformations that have been reported in the research of TRIP HEAs. One is the fcc → hcp type, found in the Co 20 Cr 20 Fe 34 Mn 20 Ni 6 [131,132] and Co 10 Cr 10 Fe 50 Mn 30 [129]. Fcc → bcc, which was reported in Co 15 Cr 10 Fe 60 Ni 15 [133].
The improved strength of TRIP HEAs is mainly due to the continuously changing volumes of all phases during deformation, whilst the increase in ductility is enabled by the enhanced strain hardening capacity granted by transformation-induced hardening of the metastable phase and dislocation hardening of the stable phase [129,131]. The strengthening effect caused by TRIP effect can be expressed by the linear combination of the two combined phases: where the σ par and σ mar represent the strength of parent phase and martensite. f par and f mar are the volume fractions of each phase. σ 1 is the strength contributed by other mechanisms except TRIP effect. Equation (35) can be used to calculate the strength of alloys containing metastable phases under any stain conditions when a satisfactory volume fraction and strength can be determined for each phase. It should be noted that the σ par and σ mar is the total strength from all the strengthening mechanism contributions, such as solid-solution strengthening, grain boundary strengthening, and dislocation strengthening. The strength of σ par and σ mar correspond to the parent phase and the martensite phases, respectively. The volume fraction of each phase can be detected by XRD, neutron diffraction [130,133,134], and EBSD maps [129,131]. However, the volume fraction changes constantly with the increase in applied stress or strain. Therefore, the volume fraction of each phase needs to be quantified at different stages. For the convenience of using equation (35), there is a function of the volume fraction evolution of martensite with strain: f mar = f (ε). And f par = 1 − f mar . When deformation-induced shear bands at their interface promote the growth of the martensitic phase, the volume evolution of martensite can be expressed as follows: where β M represents the martensite nucleating probability at shear bands intersections which depends on the SFE; f sb is the original phase shear band volume fraction; ε is the strain quantity; α sb is a parameter which depends on strain rate and represents the shear bands formation rate [135]; n is the probability of shear bands intersection and the value is 4.5 suggested by Olson and Cohen [135].

Alloy design
There have been several alloy design attempts adopting statistical techniques such as machine learning [136][137][138][139]. They are useful for the discovery of alloys within relatively well known compositional spaces [140], and for relatively simple properties such as hardness [141].

Variables for design
The present view focuses on a range of hardening mechanisms which common denominator is compositional dependence. Equation (1) can be considered as the function to maximise for, showing the compromises between the hardening terms on the right-hand side of equation (1). This is outlined next.
• σ ss directly depends the amount of solute in solid solution in the matrix, but such solute is decreased with the development of precipitation. Thus, σ ss and σ ppt compete for solute for hardening. • σ dh directly depends on the amount of dislocations per unit area, however, dislocation density triggers recrystallization at high temperatures, which is the most common mechanism to control grain size and thus define σ gb . It follows that σ dh and σ gb compete for dislocations for hardening. • σ twin and σ trip depend on alloy metastability, which in turn depend on the Gibbs free energy difference from the retained high temperature phase, and that of the lower phase. This has been recently illustrated for the titanium alloys in the work of Bignon et al [142,143].
Ignoring σ fr , which is the lowest contribution term, there are then threes contributions to determine HEA hardening: (a) The thermodynamics of phase stability. (b) The evolution of dislocation density. (c) The metastability of high-and low-temperature phases.

Design methodology
The authors are not aware of design calculations considering above three factors simultaneously. And an overall design strategy for HEA would be highly HEA-system dependant. An issue to be born in mind is that several constants in above hardening items are composition dependant, and not known for new systems; however, hardening depends on the equilibrium between those three factors. In spite of the above, it is possible to suggest an overall Figure 13. Flow diagram of HEA design based on different hardening mechanism. HTP, LTP means high-temperature phase and low-temperature phase, respectively. PTMC means phenomenological theory of martensitic crystallography (PTMC). methodology for HEA design. This is shown in the flow diagram shown in figure 13, which methodology can be summarised as follows: • A new composition has to be suggested first, along with the possible deformation conditions. Composition is intrinsically related to phase stability whereas deformation conditions trigger deformation structures and metastable phases. Some deformation structures are highly dependent upon phase metastability. • Next to the composition input is a thermodynamic stability calculations, this can be carried out with standard thermodynamic software such as thermocalc [144]. From this, the temperature dependence of the high-temperature and low-temperature phase (HTP and LTP, respectively). • As shown on the left branch of the flow diagram, the stability of a single-phase high entropy can be estimated by ensuring its sole occurrence within a temperature range. Previous work in high-temperature single-phase HEAs has shown this to be 700 to 1300 • C. • If a single phase takes place then σ ss is calculated. If not, then when a desired precipitate is found then σ ppt can be calculated. This captures the compromise between solid solution strengthening and precipitation. • The right branch of the flow diagram shows the deformation-induced structures, which are harder to predict. Their occurrence depends on thermodynamic and crystallographic factors. The former stem from the metastability of the LTP with respect to the HTP. When this is within a threshold value (e.g. 1000 J mol −1 [145]), then crystallographic coherence, i.e. the ability to accommodate the product phase within the parent phase has to be assessed. The PTMC can then be adopted [145]. Bignon et al have done this for titanium alloys [142,143]. • Once the occurrence of TRIP and TWIP is established based on the deformation-induced occurrence of martensite or twinning, then σ TRIP and σ TWIP can be calculated. • Kocks-Mecking relationships [146,147] can be employed to calculate the dislocation forest development, and from it σ dis . • The algorithm can be repeated until the alloy composition and deformation conditions match the desired output.

Conclusions
This paper summarises the strengthening mechanisms and corresponding modelling approach for HEAs. Compared to traditional dilute alloys, HEAs are distinctive for the lattice resistance and solid-solution strengthening. The severe lattice distortion of HEAs makes them more resistant to dislocation motion than their dilute counterparts, partly resulting in their outstanding properties. Solid-solution strengthening in HEAs is contributed by the multicomponent, action of several elements, where solute and solvent become indistinguishable. Fleischer's model has been used to simulate solid-solution strengthening effect with dilute element, while Labusch's model is suitable for a higher concentration solute, especially for equiatomic alloys. The other strengthening mechanism found in HEAs are quite similar to conventional dilute alloys. Therefore, many classical models used in conventional alloys can be adapted to strengthening simulation of HEAs. HEAs also display other strengthening mechanisms presents in conventional alloys, but it is the combination of them what makes them of particular interest in HEAs. Precipitation strengthening is of particular interest. The compositional complexity of HEAs can be exploited to produce bimodel, and eventually multi-model precipitate distributions. This can promote hard and soft regions improving strength and ductility. Dislocation strengthening can provide a relevant strength baseline, which combined with grain boundary strengthening can readily be tailored through thermomechanical processing.
There are two prominent strengthening mechanisms in HEAs amenable to further exploitation. Twining and TRIP. Owing to the ability to tailor metastability, and by controlling SFE, HEAs can be designed to include those effects.
The present review summarises mathematical formulations to describe the prevalent strengthening mechanisms in HEAs. The key constants for a range of compositions are also presented.