The influence of dilute aluminum and molybdenum on stacking fault and twin formation in FeNiCoCr-based high entropy alloys based on density functional theory

Stacking faults, as defects of disordered crystallographic planes, are one of the most important slipping mechanisms in the commonly seen lattice, face-centered cubic (FCC). Such defects can initiate twinning which strengthens mechanical properties, e.g. twinning-induced plasticity (TWIP), of high entropy alloys (HEAs) at cryogenic temperatures. In this work, by using density functional theory (DFT), the twinning initiated from stacking faults is discussed with regard to two different solute elements, Al and Mo, in the FeNiCoCr HEAs. Our results show that adding aluminum (Al) has noticeable enhancement of twinnability while molybdenum (Mo) only induces more stacking faults in the FeNiCoCr-based HEAs.

faults in such HEAs, by the quantification of the lattice distortion in each structure using the mean-square atomic displacement (MSAD).

Methods
Density function theory (DFT) is used to calculate the formation energies of various crystal structures in this study 17 . The plane-wave based Vienna Ab initio Simulation Package (VASP) is employed for the simulation with the projector augmented wave (PAW) method 17 as the pseudo-potential treatment. The generalized gradient approximation (GGA) parameterized by Perdew, Burke and Ernzerhof (PBE) 18 is chosen as the exchange-correlation function. Full valence electron of the composited elements (Fe, Cr, Ni, Co, Al/Mo) for these prototype alloys are considered and have valence electron number of 8, 6, 10, 9, 3/6, respectively 19 . Although usually the cutoff radius of plane-wave energy chosen as 350 eV, which is the maximum cutoff radius of Ni among considered components (269 eV) with a pre-factor of 1.3 is large enough for alloy studies, our test result in Supplementary Fig. 2 shows that if the cancellation effect is taken out of consideration, a plane-wave energy cutoff radius of 400 eV is necessary in order to eliminate the uncertainties. Monkhorst-Pack k-point spacing of 0.2 Å −1 is used in our study. The tolerance of the energy and force convergence were 1.0 × 10 −5 and 1.0 × 10 −4 , respectively. In order to represent the chemical disorder of HEAs and solid solution state of Al and Mo in the FeNiCoCr matrix, the Alloy Theoretic Automated Toolkit (ATAT) is used to build the supercells containing 96 atoms (Fe 22 Ni 22 Co 22 Cr 22 X 8 , X stand for Al/Mo) with SQS 20,21 . Supercell models were used in this study 22 , the slab models are built with vacuum added in the 111 direction. The thickness of vacuum is tested and chosen as 26 Å so that the pseudo-interactions between the repeated slabs in the 111 direction due to the periodic boundary condition (PBC) can be eliminated. Detailed vacuum thickness tests and selection can be seen in the Supplementary  Fig. 3.
The size of systems would also have effects on the calculation results. In Zhao et al. 's work 23 , they carried out the molecular dynamics (MD) simulations of NiFe solid solution alloys with various stacking fault areas. In their work, they presented the distributions of SFE in NiFe with stacking fault areas of 96.57, 536.52, 2146.08, 8584. 32,19314.72 and 34337.28 Å 2 . The standard deviation decreases rapidly with increasing stacking fault area. That means, especially for the multicomponent solid solution, larger size systems should be considered for eliminating the uncertainties of SFEs caused by local chemical environments. However, the number of atoms in DFT calculations are usually constrained at tens or hundreds. Therefore, instead of using systems as large as they used in the MD simulation, we conducted multiple samples to get the averaged values of SFEs and we believe the standard deviation would be smaller if very large model could be used in the DFT calculations. Figure 1 shows the stacking sequence of atomic planes in each structure and the atom mapping. The SQSs with FCC lattice were firstly generated, then the Shockley partial dislocations with Burger's vector (1) p were assigned to the atoms in 3 rd , 4 th and 5 th layers to create the desired structures.

Results and Discussion
In Table 1, the lattice parameters of each structure were presented to determine the effect of stacking faults on the lattices. In general, adding Mo or Al atoms into FeNiCoCr will both increase the lattice constant. The Mo-containing HEAs would have larger lattice constant than the Al-containing ones because of the larger atomic radius of Mo. When introducing partial dislocations into the FCC structures, lattice parameters are slightly elevated due to the reduction of geometrical symmetry. Compared to the theoretical determination of FeNiCoCr lattice parameters in Zhang et al. 's work 23 , adding Al atoms slightly decreases the lattice parameters while adding Mo would have the opposite effect. Molybdenum, as one of the transition metals, has similar atomic and ionic radius to the Fe, Ni, Co, Cr atoms in the matrix. However, aluminum has similar atomic radius and much smaller ionic radius. The distinct ionic size of aluminum atoms could serve as one of the reasons of decreased lattice constant. On the other hand, Mo has larger atomic and ion radius than the matrix atoms. These could result in different ways of lattice distortion. www.nature.com/scientificreports www.nature.com/scientificreports/ The cohesive energies in Table 1 are calculated using where E atom represents the energies of each species at single-atom states and E bulk is the energies of each structure at its alloying state. The energy differences are averaged by N the total number of atoms in the model. The stacking fault energies are calculated using Where E isf and E 0 the total energies of faulted and pristine supercells, respectively, A 0 is the area of stacking fault. As shown in Table 2, the surprisingly negative stacking fault energies of both samples suggest that at cryogenic temperature, faulted FCC is more stable than pristine FCC for both Al-and Mo-containing HEAs. In Zhao et al. 's work 22 , the results of the axial interaction model (AIM) and the one dimensional axial next-nearest neighbor Ising (ANNNI) model both showed very low even negative stacking fault energies for the face-centered cubic HEAs. The results from their ANNNI models and our results indicate that HCP is more favorable than FCC energetically and the total energies for both samples experience rebounds when introducing ESFs into them. Apparently, the cohesive energies of Al-containing HEAs increase dramatically with ESF structures (go to positive value). The HCP phase turns out to be the universally lowest energy state of both HEAs. However, the DFT calculations in this study were carried out with cryogenic temperatures, which are not able to predict the phase stability of elevated temperatures. In the Supplementary Fig. 4, the temperature dependent phase stability of the studied HEAs is presented. It can be seen that with the temperature increases, the FCC phase becomes more stable than the HCP phase. Notably, for all the calculated structures (ISF, ESF and HCP), the Al-containing ones always cause greater changes in cohesive energies. These evidences all could point to the fact that the phase stability of FeNiCoCr-based HEAs would be decreased by adding Al atoms.
In order to examine the statistical uncertainties caused by the local chemical environments, the SFEs by the slipping of different atomic layers were calculated. As shown in Table 3, the stacking fault energies at different fault layers were carried out. The averaged SFEs for Al 0.36 FeNiCoCr and Mo 0.36 FeNiCoCr are −2.46 meV/Å 2 and −0.28 meV/Å 2 , respectively. The individual SFEs diverse greatly and the averaged value still show that the SFE of Al 0.36 FeNiCoCr is much lower. Local chemical environments would have significant effects on the results of stacking fault energies, especially in the Al-containing HEAs. Furthermore, in order to study the uncertainties caused by atomic configurations, 10 different SQSs were taken into consideration for each proposed structure. Figure 2(c) shows the averaged value for each structure. The mean values show a clear trend that Al 0.36 FeNiCoCr has lower energy of ISF, HCP when Mo0.36FeNiCoCr has lower energy of ESF. When compare these results to the work of others, e.g. CoCrNi and FeNiCoCr alloys [24][25][26] , it seems that aluminum tends to decrease the SFEs while molybdenum tends to increase the SFEs.  www.nature.com/scientificreports www.nature.com/scientificreports/ In addition to the cohesive energy calculations, the minimum energy paths and the energy barriers during the shifting of atomic planes for each slipping mechanism were investigated. The activation energies are presented in Fig. 2(a,b). The nudged-band elastic (NEB) method is used in the energy barrier calculation. The combination of the climbing image method and the better tangent definition gives more accurate saddle points with fewer images and larger displacement intervals, compared to the originally compiled elastic band (EB) method in VASP 19,[27][28][29][30][31][32] .
As shown in Fig. 2(a), the first saddle of creating ISF in Al 0.36 FeNiCoCr is around 0.023 eV/Å 2 . The Mo-containing samples have slightly lower barrier than the Al-containing samples. From the second saddles for ESFs, the energy barriers basically keep the same value as those for ISFs. Either the formation energies or the activation energies of Mo-containing HEAs have small changes from one structure to another. This validates what is stated in the previous discussion, the Al-containing HEAs are less stable than that the Mo-containing ones, as per se, adding Al enhances the meta-stability and propensity of forming stacking faults and nano-twins in the www.nature.com/scientificreports www.nature.com/scientificreports/ FeNiCoCr-based HEAs 32 . The second saddle in Fig. 2(b) shows great propensity of both samples of forming HCP phases in FCC through generating stacking faults although in experiments the HCP phases are not usually probed because the continuous deformation would eventually induce large amount of stacking faults or nano-twins.
Lattice distortion could be the origin why Al atoms decrease the phase stability of FeNiCoCr-based HEAs. The mean square atomic displacement (MSAD) is extracted from the DFT calculations for each atom. The MSAD in this work was carried out using where the x and x 0 represent the coordinates of atoms at the initial and terminal states, respectively. For the faulted configurations discussed in this work, the cell-shape optimized models with rigid displacements inserted into the optimized FCC structure were used as the initial states and then in the same supercells ionic positions optimization were carried out as the terminal states to ensure the coordinates in the MSAD determination are in the same space (supercell). From Fig. 3, it is obvious that all the atoms except nickel and chromium atoms have larger MSAD in the Al 0. 36    www.nature.com/scientificreports www.nature.com/scientificreports/ the two faulted atomic planes (nearest neighbors in the same close-packed plane are not counted). While in the ISF structures, the Al/Mo-Fe neighbor numbers are reduced to one which reduces the atomic offset of iron atoms.
In Wen et al. 's publication of generalized stacking fault energy, they found that in the nickel based alloys, the atomic radius difference ΔR as well as valence electron difference ΔVEC between solute atoms and solvent atoms all have influence on stacking fault energy in FCC alloys 34 . In their work, larger ΔR result in lower stacking fault energy and larger ΔVEC also tend to drag down the stacking fault energy. As shown in Table 4, in our current study, when we treat solvent atoms valence electron count as the mean value of each element in the quaternary matrix alloy, the ΔVEC are 5.25 and 2.25 for Al 0.36 FeNiCoCr and Mo 0.36 FeNiCoCr, respectively. And the behaviors of stacking fault energies with respect to the ΔVEC agree with their proposition. The ΔR of the Al and Mo are 17.5 pm and 13.5 pm, respectively, are also consistent with their findings. For the species considered in this work, if geometry sizes of atoms are considered, we think the difference between the radius of ionic radius which are ΔR ionic = −21 pm and −3 pm for Al and Mo, respectively, are also one of the factors that can affect the stacking fault energy of alloy systems. The smaller size of Al nuclei compared to Mo shows ease and tendency of forming stacking faults. In Ma et al. 's report 11 , they showed that when the mole fraction of Al in FeNiCoCr-based HEA is above 0.45, the material would tend to form BCC/B2 phase in it. That suggests that although Al itself is naturally FCC metal, adding Al to an FCC alloy matrix will induce BCC/B2 phase separation, says, adding Al will stabilize the BCC structure in FCC HEAs. In King et al. 's work of Al x FeNiCoCr alloys at atomic scale, they give a detailed study of Al content effects on the phase stability of HEAs 35 . In their work, the HEAs tend to remain disordered FCC when the Al content stays below 0.5. And they predict a transition between and FCC-based structure to the partially ordered BCC structure with increasing Al content which has great agreement with the experiments. They also investigated the ordering of the FeNiCoCr FCC phase in a Cu3Au (L1 2 )-type structure which is devoid of Al. Note that if take the partially ordered structures in the FCC HEAs, the results of SFEs of our systems would differ 36 . However, we here use a lower concentration of Al contents and discuss the influence of Al and Mo atoms in a completely random solid solution.
In addition, molybdenum, as a naturally BCC metal, could also stabilize the BCC structure in FCC HEAs. We think that one of the reasons for the tendency of forming twins or stacking faults caused by Al or Mo is the ratio of ionic radii to van der waals radii, k = 0.350 and 0.489 for Al and Mo, respectively.

Conclusion
In this work, the ISFEs, ESFEs and PTEs and their activation energies in single-phase FCC FeNiCoCr-based HEAs with solid solution components (Al/Mo) were calculated using DFT. The MSAD, ionic radius and VEC were used in our study to find the effects on stacking fault formation by adding different types of solute atoms into the quaternary HEA matrix. The temperature-dependent phase stability of FCC and HCP phase in this HEA was also calculated. The negative ISFE shows the great propensity of forming stacking faults in the HEAs used in our work. The even lower ISFE caused by Al solutes took responsibility for the high propensity of twinning in Al 0.36 FeNiCoCr which agrees with the experimental observation from the others' work 37 . The higher ESFE (obvious in Al 0.36 FeNiCoCr while slight in Mo 0.36 FeNiCoCr) implies that forming twins from growth of stacking faults may not be the energetically favorable way for deformation. On the contrary, the surprisingly low PTE and energy barriers of both samples show strong evidence that parallel stacking faults which tends to form one after another by skipping one atomic plane is more likely to be the deformation path for twin growth. We also showed that the lattice distortion is released (especially Fe atoms) after introducing stacking faults which partially contributes to the negativity of SFE. From the aspect of atomic electron structures and geometry size of atoms, higher value of differences in nucleus radius and valence electron numbers between solute atoms and solvent atoms tends to form more stacking faults in FCC HEAs. To sum up, the solid solution HEAs with Al or Mo are both likely to form annealing and deformation stacking faults while Al 0.36 FeNiCoCr would reserve more stacking faults and nano-twins during annealing and form more nano-twins during deformation from the perspective of ionic to atomic radius ratio and valence electron count.

Data Availability
The authors declare that the data supporting the findings of this study are available within the paper and supplement. Also, the data that support the plots within this paper are available from the corresponding author upon reasonable request.