Microstructural evolution and mechanical properties of FeCoCrNiCu high entropy alloys: a microstructure-based constitutive model and a molecular dynamics simulation study

High entropy alloys (HEAs) attract remarkable attention due to the excellent mechanical performance. However, the origins of their high strength and toughness compared with those of the traditional alloys are still hardly revealed. Here, using a microstructure-based constitutive model and molecular dynamics (MD) simulation, we investigate the unique mechanical behavior and microstructure evolution of FeCoCrNiCu HEAs during the indentation. Due to the interaction between the dislocation and solution, the high dislocation density in FeCoCrNiCu leads to strong work hardening. Plentiful slip systems are stimulated, leading to the good plasticity of FeCoCrNiCu. The plastic deformation of FeCoCrNiCu is basically affected by the motion of dislocation loops. The prismatic dislocation loops inside FeCoCrNiCu are formed by the dislocations with the Burgers vectors of a6[1¯12¯]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a \over 6}\left[ {\bar 11\bar 2} \right]$$\end{document} and a6[11¯2]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${a \over 6}\left[ {1\bar 12} \right]$$\end{document}, which interact with each other, and then emit along the 〈111〉 slip direction. In addition, the mechanical properties of FeCoCrNiCu HEA can be predicted by constructing the microstructure-based constitutive model, which is identified according to the evolution of the dislocation density and the stress-strain curve. Strong dislocation strengthening and remarkable lattice distortion strengthening occur in the deformation process of FeCoCrNiCu, and improve the strength. Therefore, the origins of high strength and high toughness in FeCoCrNiCu HEAs come from lattice distortion strengthening and the more activable slip systems compared with Cu. These results accelerate the discovery of HEAs with excellent mechanical properties, and provide a valuable reference for the industrial application of HEAs.


Introduction
The concept of high entropy alloys (HEAs) was put forward by two research teams in 2004 [1] with two main ideas, which mainly focused on the unrevealed core zones of the multi-element phase diagram. HEAs with several of principal elements in equiatomic and near-equiatomic ratios have superior properties [2] . These excellent mechanical properties are caused by high mixing entropy and severe lattice distortion in HEAs, which are not available in traditional alloys [3] . Nöhring and Curtin [4] indicated that owing to severe lattice distortion caused by the mixing of different atomic sizes, slow dislocation activity occurs in HEAs, which is not only conducive to the formation and stability of the solid-solution phase, but also brings enough strength and ductility to HEAs. Instead of formation of intermetallic compounds, HEAs form single phase solid solutions, including face-centered cubic (FCC) or body-centered cubic (BCC) solid solutions. Due to the effect of lattice distortion, the strong lattice resistance, which hinders dislocation slip occuring in HEAs, results in solid solution strengthening, which plays a vital part in improving the yield strength of HEAs [5] .
As a typical type with single-phase FCC, the FeCoCrNiCu HEA has good thermal stability and excellent mechanical properties [6] . Due to outstanding thermal stability, the FeCoCrNiCu HEA has been widely investigated in high temperature environment [7] . Furthermore, it has good wear resistance and high hardness, caused by the addition of Cu [8] . However, the deformation behaviors of the FeCoCrNiCu HEA were rarely revealed, which causes the lack of theory about the FeCoCrNiCu HEA system, thereby narrowing the scope of its industrial applications. Therefore, it is necessary to reveal the deformation mechanism and mechanical response of the FeCoCrNiCu HEA in detail.
With rapid development of processing equipment and technology, an increasing number of material mechanics researchers are focusing on nanometer or micron processing. As a common technique, nanoindentation is utilized to investigate the mechanical response of workpieces [9] . The research [10] represents that the nano-scale contact is dominated by atomic phenomena, and the plastic deformation of workpiece occurs in the form of dislocation nucleation during the nanoindentation process. Since it is extremely difficult to observe this nano-scale phenomenon on-site during the experiment, the introduction of computational tools will help us understand the deformation mechanism at this scale. Researchers generally employ molecular dynamics (MD) simulation to investigate the deformation behavior of materials [11] . In addition, these studies [12][13] further displayed that MD simulation is an appropriate approach to study the mechanical response of materials during nanoindentation.
Although there have been some experimental studies and atomic simulation studies on the FeCoCrNiCu HEA, the theoretical research on its strengthening mechanism is still lacking. Therefore, it is necessary to propose a microstructure-based predictive model to study the stress-strain response in the FeCoCrNiCu HEA. In the paper, we contrastively analyze the nanoindentation responses of the FeCoCrNiCu HEA and single crystal Cu via MD simulations, to better illustrate the superiority of the FeCoCrNiCu HEA. In addition, a microstructure-based constitutive model is established to further investigate the strengthening mechanism in the HEA.

MD model
The nanoindentation model has been constructed by the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). For reference purposes, the FCC crystal structured Cu is utilized as a comparison. Figures 1(a) and 1(b) show the models of nanoindentation into the single crystal FeCoCrNiCu HEA and Cu, respectively. Here, two MD models include (i) FeCoCrNiCu HEA workpiece + virtual indenter and (ii) Cu workpiece + virtual indenter. Considering the effect of indenter size, the indenter radius is set as 3.5 nm, as represented in Fig. 1. The boundary layer atoms with 1 nm thickness are fixed to prevent the whole substrate from any translational or rotational movement. The thermostat layer atoms with 1 nm thickness perform the function of heat dissipation via the velocity scaling method to keep the temperature remaining at 293 K. The Newtonian layer atoms are in accord with the classic Newton's second law [14] . The simulation parameters of nanoindentation are shown in Table 1.   The potential used is a set of embedded atom method (EAM) model interatomic potential [15] that allows for an explicit characterization to quinary alloy. It contains the potentials of the pure elements, and the lattice constant of Cu is 0.362 nm which is close to the experimental value [16] . The atoms of the HEA are randomly distributed at approximately equal proportions on the FCC lattice points, with the lattice constant of 0.356 nm which is consistent with the research in Ref. [17].
Both relaxation and indentation are performed in canonical ensemble (NVT). The virtual indenter created is a repulsive sphere, which interacts with the sample atoms by the potential [18] where R stands for the indenter radius (3.5 nm), and r denotes the distance of a workpiece atom to the indenter center. According to Ref. [19], the force constant K can be set to be 300 eV/nm 2 for the surface of the indenter. As a visualization software, OVITO [20] is employed to observe the atomistic configurations in the current simulation.

MD results and analysis
3.1 Mechanical properties Figure 2(a) shows the load-displacement curves of the single crystal FeCoCrNiCu HEA and Cu. Based on the curves and the Oliver-Pharr method [21] , the hardness H and Young's modulus E of the FeCoCrNiCu HEA are calculated. The hardness H and Young's modulus E are obtained by the following formula [22] : where S is the contact stiffness of the unloading part of the curve, which can be calculated by S = dp dh h = h max . P max stands for the maximum force during indentation. R represents the radius of the spherical virtual indenter. The related factor of the indenter shape ε can be set to be 0.75 for the spherical indenter [21] . Poisson's ratio v of the HEA can be regarded as 0.3 according to Ref. [21]. The constant β related to the shape of indenter can be regarded as 1 for the spherical indenter [21] .
The continuum model would consider that the contact area is complete. However, the surface at the atomic structure is discontinuous. Here, the total area between the indenter and workpiece, A T , can be obtained by the equation [23] A T = N T A α , where N T is the number of atoms interacting, and A α is the average contact area of an individual atom. According to the calculation, the deviation between A C (36.4 nm 2 ) and A T (35.48 nm 2 ) is greatly slight, which further proves the validity of A C .
The results show that Young's modulus of the FeCoCrNiCu HEA is 201 GPa with a small deviation of 4% from the experiment [24] . In addition, Young's modulus of the Cu crystal is 114.8 GPa, which is very close to the experimental value of 114 GPa [25] . The hardness of the FeCoCrNiCu HEA calculated by the current simulation is 17.2 GPa, which is slightly higher than the previous experimental result [26] . The corresponding hardness of the single crystal Cu is 11.9 GPa, which is marginally higher than that obtained in Ref. [27]. Because of the strong effect of indentation size [28] and the scale difference between the simulation and the experiment, there is a slight deviation between our simulation results and the experimental values. For instance, in the microscale, the precipitation of grain boundary void is the major observable defect, which seriously affects the mechanical properties and weakens the workpiece. However, with the decrease in the dimension, when it is smaller than the scale of typical grain boundary voids, the main visible defects are dislocations instead of voids [29] . Therefore, as the scale gradually increases, the number of defect types in the workpiece will augment. Another reason for the deviation is that the simulated models here are perfect single crystals, while the experimental workpieces [26][27] contain various defects. The above statements further indicate the accuracy of our MD simulation results.
In Fig. 2(a), the loading and unloading curves of the HEA and Cu represent some differences. It can be clearly observed that the fluctuation and the magnitude of the nanoindentation force F z in FeCoCrNiCu are larger than those of Cu. This fluctuation is caused by the nucleation and motion of dislocations in the workpiece [30] . During the loading period, new dislocations are generated in the materials due to the action of the indenter, which leads to the increasing force.
The dislocation density of the workpieces during indentation is calculated, which provides a direct reference to explain the flow stress of deformed materials [30] , as shown in Fig. 2(b). Obviously, the dislocation density inside the HEA is larger, which interprets the first-order effect of the traditional metallurgical alloy system (i.e., the strength of alloy is proportional to the dislocation density) [31] . Due to high density of dislocations in the HEA, the interactions between the dislocation-solute atoms lead to strong work hardening, and plentiful dislocations nucleate and move in the activated slip system, leading to good plasticity. This implies that the dislocation mechanism has a significant effect on the plastic deformation and strengthening mechanism of the HEA, which will be explained later.

Evolution of interior defects
, and 3(e) represent the evolution of interior defects in the workpieces during the penetration of indenter. All atoms are colored by the common neighbor analysis (CNA) value, where red atoms, blue atoms, and green atoms represent hexagonal close-packed (HCP) atoms, BCC atoms, and FCC atoms, respectively. In the early stage of nanoindentation (i.e., at the indentation depth of 0.3 nm), a few defects are generated under the indenter, which represents the initiation of dislocation nucleation. As the indenter deepens, more dislocations slip along the {111} slip plane. It can be observed that the emission of dislocation loop in place is a common feature of the two materials, but the difference is that the dislocation loop and stacking fault (SF) in FeCoCrNiCu are greater compared with Cu, which is attributed to the low stacking fault energy (SFE) of the HEA [30] .
Figures 3(c) and 3(f) show the variation curves of the number and percentage of different atomic types in the workpiece with the indentation depth, where HCP atoms and other atoms represent SF atoms and surface atoms, respectively. It can be observed from Fig. 3(c) that with the increase in the indentation depth, the number of HCP atoms and other atoms increases, accompanied by the decrease in the FCC atoms. It means that more SFs are converted into HCP structures. The number of HCP atoms in the FeCoCrNiCu HEA is larger. Figure 3(f) gives the percentage of atomic types in the HEA and Cu. It can be seen that the proportion of HCP atoms in the HEA is higher during the indentation. Therefore, the work hardening induced by nanoindentation in the FeCoCrNiCu HEA is still closely related to the microstructure evolution.
The formation of prismatic dislocation loops in the FeCoCrNiCu HEA usually goes through three stages, i.e., embryonic prismatic dislocation loop, "lasso"-like dislocation loop, and prismatic dislocation loop, as proposed by Remington et al. [32] . Figures 4(a)-4(c) show snapshots of these dislocation loops. At the initial stage of indentation, the embryonic prismatic dislocation loops are formed on the subsurface of the HEA, as displayed in Fig. 4(a). With the indenter penetrating, the two edges of embryonic prismatic dislocation loop interact with each other to form a "lasso"-like dislocation loop, and finally transform into an independent prismatic dislocation loop, as represented in Figs. 4(b) and 4(c). Interestingly, the cross-slip hardly appears in the formation process of prismatic dislocation loops in the FeCoCrNiCu HEA. This is different from the result of Remington's work [32] , where the prismatic dislocation loops were induced following the cross-slip in BCC Ta. The phenomenon indicates that the screw-type Shockley dislocation hardly has cross-slip in the FeCoCrNiCu HEA. It can be seen from Fig. 4(b) that no cross-slip appears, in conformity with the previous work [33] . Only the interaction between Shockley partial dislocations is observed, and there is no interaction between stair-rods.
The transformation mechanism from "lasso"-like dislocation loops to prismatic dislocation   loops in this work is very different from that mentioned by Remington et al. [32] . Here, the dislocation reaction of the transition from "lasso"-like dislocation loops to prismatic dislocation loops is investigated in detail. Thompson tetrahedron [34] is constructed to track the Burgers vectors, as shown in Fig. 4(d). All atoms are hidden to display dislocation lines. Figure 4(e) manifests "lasso"-like dislocation loops, and the Burgers vectors of a 6 [112] and a 6 [112] are marked by and , respectively. It can be observed from Fig. 4(e) that a 6 [112] and a 6 [112] move towards each other, and finally contact and interact to form the prismatic dislocation loops, which can be expressed by a 6 [112]+ a 6 [112] = 0 ( + =0). Due to the interaction of dislocations and , the prismatic dislocations are finally emitted from the indentation region, as shown in Fig. 4(f). ). The prismatic dislocation loops emit along the 111 slip direction, which is consistent with the traditional metals [33] .

Deformation behaviors
Figures 5(a)-5(d) show the distributions of shear strain in the workpieces with different indentation depths. The deformation characteristics of the FeCoCrNiCu HEA are slightly different from those of metallic crystals, because the mechanical behavior inside the HEA is also affected by lattice distortion [2] . The fact that each element in the FeCoCrNiCu HEA may ignore the chemical order and occupy the lattice position, results in severe lattice distortion and affects the deformation behavior of the HEA during the indentation [2] . We analyze the shear strain generated during the nanoindentation process based on this. All atoms are colored by the magnitude of the partial shear strain. The shear region is nucleated in the contact area between the indenter and the workpiece, and then spreads into the workpiece. The inhomogeneous distribution of shear strain in the workpieces is attributed to the strain gradient effect, resulting in diverse strain hardening along the indentation direction. Comparatively, the shear area inside the FeCoCrNiCu HEA is slightly larger, which is mainly caused by dislocation slip in the workpieces.
In addition, the von Mises stress distribution of Figs. 5(e)-5(h) illustrates the difference in deformation behavior between HEA and Cu during indentation, where red depicts higher stress and blue represents lower stress. It can be observed that the local stress value of the FeCoCrNiCu HEA is higher because of the lattice distortion, which leads to high strength of the HEA [35] . During the nanoindentation, the atomic displacement analysis is carried out in the indentation area. The atomic displacement distributions of the FeCoCrNiCu HEA and Cu are shown in Figs. 5(i)-5(l). As the indenter goes deeper, the atom displacement gradually augments. Here, red atoms represent higher local displacement values and indicate severe deformation. It can be observed that the deformation of FeCoCrNiCu is sightly more severe, which shows better plasticity of the HEA.
In order to further study the plasticity of FeCoCrNiCu in the indentation process, the typical snapshots of the atomic displacement vector are presented in Figs. 6(a)-6(d). It displays that the atoms move in a curve around the indenter because of the influence of the shear stress flow. As the indentation progresses, the atoms on the surface of the HEA shift in the opposite direction of the indentation, causing the pile-ups generated on the surface of the HEA. Furthermore, the lattice rotation can be observed in the pile-up area. For the FeCoCrNiCu HEA, the lattice rotation under the indenter has great influence on the dislocation density [36] . In order to analyze the plastic deformation of the zone around the indentation, the surface morphology of single crystal FeCoCrNiCu HEA and Cu at the depths of 1.5 nm and 3 nm is characterized, as shown in Figs. 6(e)-6(h), where all atoms are colored according to the atomic height.
With the increase in the indentation depth, the pile-ups are generated on the surface of the workpieces. The presence of pile-ups is mainly caused by the movement of dislocations. It can be observed from the yx-plane view of the workpieces that atom pile-ups present strong inhomogeneity and asymmetry. Since there are more dislocation loops generated in the HEA,   the dislocation loops slip further toward the orientation 111 parallel to the surface, forming the rectangular steps, as shown in Figs. 6(e)-6(f). However, since there are more sessile dislocations in single crystal Cu and fewer dislocations for slippage, pile-ups are mainly generated around the pit, as shown in Figs. 6(g)-6(h). Comparative analysis of the results manifests that FeCoCrNiCu has good ductility.

Theoretical modeling
The microstructures and their impact on mechanical properties are exhibited in the MD simulation. However, the intrinsic correlation between the microstructure and mechanical properties is still not revealed. Here, a microstructure-based constitutive model is established to illustrate the difference of the mechanical performance between the single crystal HEA and Cu. The constitutive model is as follows.
The total strain rate consists of the elastic strain rate and the plastic strain rate, According to Hooke's law, the elastic part satisfies the linear relation, where M is the elastic compliance tensor, andσ is the stress rate tensor. The plastic strain rate is indicated asε where σ = σ ij − σ kk δ ij /3,ε p = 2ε pijεpij /3, and σ = 3σ ij σ ij /2 denote the deviatoric stress tensor, the von Mises equivalent plastic strain rate, and the von Mises equivalent stress, respectively. The relationship betweenε p and σ iṡ In the above equation,ε denotes the equivalent strain rate, and m is the sensitivity factor of strain rate, whose value is usually taken as 20. σ flow expresses the flow stress, which is mainly affected by lattice distortion, and dislocation for the FeCoCrNiCu HEA in the present work is shown as where σ d denotes the slip resistance of dislocation interaction, indicated as in which M , α, µ, b, and ρ represent the Taylor factor, the empirical constant, the shear modulus, the magnitude of Burgers vector, and the dislocation density, respectively. Due to the Kocks-Mecking model [37] , the evolution of the dislocation density can be shown as where d denotes the grain size, ψ is the proportionality factor, k 20 represents the dynamic annihilation constant,ε 0 is the reference strain rate, and k e = de d 2 expresses the dynamic annihilation coefficient of dislocations [38] , in which d e is the critical grain size and indicates that dynamic recovery begins to intensify [39] . The related parameters during the process of dislocation evolution of the FeCoCrNiCu HEA and Cu metal are shown in Table 2 [38][39] . For HEAs, the lattice distortion strengthening caused by the shear modulus mismatch and size mismatch between different elements cannot be ignored. According to Vegard's law and the previous work [40] , the contribution of lattice distortion to flow stress consists of the superposition of individual contributions of each element, which can be expressed as where c i is the concentration of element i, and σ i ss denotes the effect of element i on lattice distortion strengthening which can be written as where A is a constant, usually taken as 0.04. µ = n i c i µ i represents the shear modulus of the FeCoCrNiCu HEA, which is constituted by the shear modulus µ i of every single element. δ i is the parameter of mismatch which is related to atom size mismatch δr i and modulus mismatch δµ i and can be written as where ξ = 1 in FCC metals, and ξ = 2.5 in BCC metals [41] . The value of β is connected with the type of dislocations. For screw dislocations, 2 < β < 4, and for edge dislocations, β 16. In addition, because the HEA containing n elements can be regarded as composed of (n − 1)-principal-matrix alloy and the additional element i, the atom size mismatch δr i and modulus mismatch δµ i of element i can be written as where δr ave ijklm and δµ ave ijklm represent the average size mismatch and the average modulus mismatch of the ijklm five-element HEA, and for a similar reason, δr ave jklm and δµ ave jklm represent those of the jklm four-element HEA. The average size mismatch and average modulus mismatch can be expressed as where δr ij and δµ ij are the size mismatch and the modulus mismatch between element i and element j, respectively. They can be expressed as where r i is the atomic radius, and µ i is the shear modulus of the pure metal i. Hence, according to the above deduction, the lattice distortion strengthening of HEA can be obtained. The relevant physical parameters of corresponding elements in the FeCoCrNiCu HEA are shown in Table 3 [42] . Based on the above analysis, the evolution of dislocation density as the strain increases, and the stress-strain curve of the FeCoCrNiCu HEA and Cu metal are shown in Fig. 7. For these two materials, there is slight difference of dislocation density evolution in Fig. 7(a). The variation trend of dislocation density evolution in Fig. 7(a) is consistent with that in Fig. 2 from MD simulation. For single crystal Cu, the flow stress primarily derives from the sliding resistance formed by dislocation interaction. However, because of the difference between the shear moduli of the FeCoCrNiCu HEA and Cu metal, there is huge distinction between the contribution of dislocation strengthening to the overall flow stress according to Eq. (9). From the empirical formula µ = n i c i µ i , it can be found that the shear modulus of the FeCoCrNiCu HEA is greater than that of Cu, which is consistent with the experiments [43] and calculated results in our MD simulation. Compared with Cu, the prominently greater shear modulus of the FeCoCrNiCu HEA leads to the significantly higher dislocation strengthening and strain hardening, as shown in Fig. 7(b). The lattice distortion strengthening is non-negligible and unique compared with traditional alloys. For the FeCoCrNiCu HEA, the calculated lattice distortion strengthening is independent of the strain rate and reaches 460 MPa, which improves the yield strength. In addition, we would integrate the constitutive model into the finite-element analysis in the subsequent work, which could provide an opportunity to visualize the differences of the strain field and stress fields in Cu and HEA [44] .

Summary
In the current work, the nanoindentation models of single crystal FeCoCrNiCu HEA and Cu are constructed using MD simulations. Furthermore, we compare the difference of mechanical properties between Cu and HEA from the perspective of theoretical model. It is identified that the origins of high strength and toughness in the FeCoCrNiCu HEA come from the lattice distortion strengthening and more activable slip systems compared with Cu, respectively. The conclusions are as follows: (i) Compared with Cu, high dislocation density and high loading force occur in the FeCoCrNiCu HEA during indentation. These results indicate that the HEA has high strength.
(ii) The plastic deformations in the FeCoCrNiCu HEA are mainly in the form of emission prismatic dislocation loops, and there are more prismatic dislocation loops in HEA. In addition, it is observed that the Shockley dislocation in the HEA can hardly cross slip during the transition from "lasso"-like dislocation loop to prismatic dislocation loop.
(iii) The shear strain zone in the FeCoCrNiCu HEA is slightly larger than that of Cu, due to the long distance of dislocation slip and more dislocation loops emitted. Furthermore, there are many SFs on the subsurface of the HEA, resulting in a mass of pile-ups on the surface. This further indicates that the HEA has good fatigue resistance. Due to the effect of lattice distortion, the high local stress occurring inside the HEA leads to higher strength.
(iv) Based on the analysis of microstructure-based constitutive model in Section 4, it is found that the dislocation evolution between the FeCoCrNiCu HEA and Cu is not obvious, which is consistent with our MD simulation results. However, the shear modulus of FeCoCrNiCu is larger than that of Cu due to the element composition principle, which induces the greater dislocation strengthening on strain hardening in the HEA. In addition, the lattice distortion of FeCoCrNiCu is significantly higher than that of Cu, resulting in the higher yield stress.
Our investigation focuses on the unique deformation behavior of FeCoCrNiCu compared with Cu. Due to the severe lattice distortion, higher lattice distortion strengthening in FeCoCrNiCu leads to greater strength, while more activated slip systems in FeCoCrNiCu contribute to better plasticity. This study reveals the deformation mechanism and good mechanical properties of the FeCoCrNiCu HEA from the atomic level and theoretical perspective, providing guidance for the future exploration of new materials and the industrial application of HEAs.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.