Phase field simulations of FCC to BCC phase transformation in (Al)CrFeNi medium entropy alloys

Microstructure simulations for quaternary alloys are still a challenge, although it is of high importance for alloy development. This work presents a Phase field (PF) approach capable of resolving phase transformation in a multicomponent system with a simple and effective way to include the thermodynamic and kinetic information for such a complex system. The microstructure evolution during diffusional transformation between FCC and BCC phase at 700 °C for AlCrFeNi alloys was simulated, accounting for composition dependence and off-diagonal terms in the diffusion tensor. The reliability of the presented PF method is validated by comparing the 1-D simulation results with simulations by Diffusion Module (DICTRA) of Thermo-Calc Software. Additionally, 2-D PF simulations of precipitate growth and Ostwald ripening are performed for different alloy systems, and the coarsening behavior is compared. Results showed that thermodynamic and kinetic information is accurately described in the applied PF method. The simulation results show that the diffusion behavior is influenced evidently by variations in the amounts of the different elements in the system. These findings demonstrate the necessity of applying accurate thermodynamic and kinetic models to fully understand the complex interdiffusion behavior in high and medium entropy alloys.

is a widely studied system, showing promising features and properties, such as excellent compressive properties, good yield strength, and good plastic deformation ability (Wang, et al. 2008;. However, Co is an expensive component. To reduce the cost of production, researchers investigated the Co-free AlCrFeNi MEAs (Chen, et al. 2017;Dong, et al. 2016;Jiang, et al. 2019;Jin, et al. 2019), which also exhibits a good combination of strength and ductility with adjusting the alloy's compositions. For such HEAs/MEAs, there is a wide window for property and composition optimization. However, with the vast amount of influencing parameters, the study of phase transformations and microstructure evolution can become highly complex. Therefore, besides experiments, efficient and accurate microstructure evolution simulations are prerequisites to reduce the design time and cost for new multicomponent alloys.
The Phase field (PF) method is a powerful tool to investigate the microstructure evolution of alloys. However, the application of the PF method to multicomponent alloys, such as Al-Cr-Fe-Ni, is challenging and still rare because of the enormous thermodynamic and kinetic information required in the simulations. To date, several ways have been applied to introduce composition-dependent thermodynamic information into PF simulations. For example, Taylor's second-order expansion about the equilibrium mole fractions (Jokisaari, et al. 2017;Choudhury, et al. 2015) was used to fit the Gibbs free energy densities of the different phases. However, while the mathematical treatment is efficient and straightforward for binary systems, a truncated Taylor series is only valid near the expansion point, and the second-order expansion often leads to unphysical molar fraction values in the PF simulations when applied for ternary or even higher-order systems. As an alternative, composition and temperature-dependent functions constructed based on the CALculation of PHAse Diagrams (CALPHAD) approach (Kitashima 2008) were used to calculate the Gibbs free energy densities and diffusion potentials in the PF model for multiphase, multicomponent systems. The CALPHAD Gibbs energies are available for most binary and ternary systems and many multicomponent alloys (Kitashima 2008). However, for some CALPHAD models, such as those based on a sublattice or orderdisorder model, an explicit description as a function of the local molar fractions of the different elements (as generally required in PF simulations) is only developed recently (Schwen, et al. 2021). The proposed sublattice phase field model permits direct use of the CALPHAD sublattice model expressions for an arbitrary number of constituents and sublattices; however, such an approach requires that the coefficients of the CAL-PHAD Gibbs energy expressions are publicly available, while for many MEAs/HEAs systems there is only a description available in the commercial databases which do not give immediate access to the Gibbs energy coefficients. Moreover, extra equations are needed compared with the traditional KKS model. Thus for quaternary or quinary alloys, the computational cost may become high. Another alternative is to precompute the thermodynamic and kinetic data for discrete compositions and store them in tables (Heulens, et al. 2011;Chatterjee and Moelans 2021) for use in the PF model. However, a disadvantage for such method is that for multicomponent alloys (four, five, or more elements), the number of data points generated will be massive, which consequently increases the computational cost. Moreover, getting an accurate discrete thermodynamic and kinetic description for HEAs/MEAs is becoming challenging when more than three elements are considered in the simulation, as the size of the dataset increases dramatically with the number of elements. PF simulations for HEAs/MEAs are thus still non-trivial and even close to impossible. To our knowledge, only limited work ) has so far reported using PF numerical modeling for HEAs/MEAs is available. However, the PF model used in that paper was simplified to a pseudo -binary alloy PF model, and simplified Gibbs free energy expressions were used. Recently, (Coutinho, et al. 2020) introduced the tensor completion method for efficient use of thermodynamic information into the PF model. However, the proposed method of sampling data into thermodynamic functions is still quite complicated and may not always be necessary.
This paper aims to illustrate the application of an intuitive and efficient coupling approach for quaternary systems. We show that if the composition domain is appropriately chosen, simple polynomial fitting of the Gibbs energies allows capturing the composition dependence of the thermodynamic and diffusion data accurately for simulations of diffusional transformation in 2-phase AlCrFeNi systems. In "Phase field method and parameters determination" section, the PF model and parameters are introduced. The procedure to include the thermodynamic and kinetic information from CALPHAD databases in the PF simulations is given in section "Fitting compotion-dependent Gibbs free energy densities data". Finally, to illustrate and validate the use of the approach, the diffusional controlled phase transformation between the FCC and BCC phases for AlCr-FeNi alloys at 700 °C is simulated using the PF approach, and the results are compared with results obtained using the Diffusion Module (DICTRA) of Thermo-Calc Software with the same thermodynamic and kinetic database. The DICTRA software (J- O Andersson, et al. 2002) can be used to perform one-dimensonal simulations of multicomponent diffusion-controlled transformations. Moreover, 2D simulations of growth and Ostwald-ripening of spherical BCC precipitates from an FCC matrix were conducted, and the effect of compositional changes upon the growth rate of the precipitates was investigated. Our approach can potentially be extended to arbitrary alloy systems if the corresponding thermodynamic database is available, paving the way for simulationbased microstructure optimization for HEAs/MEAs.

Phase field method and parameters determination
Modeling microstructure evolution requires tracking the migration of the phase or grain boundaries between different regions. The PF method tackles this process by treating interfaces as diffuse and having a finite width (Chen 2002;Moelans et al., 2008a, b, c), which has been proven to be practical and powerful in simulating grain growth (Moelans et al., 2008a, Moelans et al., 2008b, phase transformations (Mohanty, et al. 2008;Wu, et al. 2004) and solidifications (Mi, et al. 2019;Kobayashi 1993). For this reason, the PF method has been selected to model microstructure evolution in the AlCrFeNi system in this work. (Kim, et al. 1999) introduced a general type of phase PF model for binary alloys, usually called the KKS model. Later (Kim 2007;Kim, et al. 2004) further extended the PF method into multicomponent multiphase systems. (Moelans et al., 2008a;Moelans 2011) introduced a new type of interpolation function allowing for quantitative PF model of multiphase systems, which is not only in grains and at interfaces thermodynamic consistent but also guarantees thermodynamic consistency in multi-junctions. Based on these models, a PF model was implemented for the AlCrFeNi system in the present work. The temperature and pressure in the system are kept constant in this study. In the AlCrFeNi PF model, the local mole fractions of Al, Cr, and Fe are represented by the variables x Al , x Cr and x Fe . The mole fraction of Ni is obtained as x Ni = 1-x Al -x Cr -x Fe . Furthermore, two non-conserved order parameters are introduced to represent the FCC and BCC phases, namely η BCC and η FCC . Within the BCC phase, η BCC =1 and η FCC =0, while in FCC phase, η BCC =0 and η FCC =1. At the interface between the two phases, the order parameters show a diffuse transition between the values 0 and 1. For a two-phase system, one order parameter would have been sufficient to distinguish between the two phases. However, since in future research, we aim to extend the model towards multiphase systems, the given representation was chosen as it can easily be extended to include more phases.
The total free energy F of the system is formulated as a function of the order parameters representing the different phases η BCC and η FCC and the molar fractions as: Where f FCC and f BCC are the Gibbs free energy densities for FCC and BCC phases. The first two terms in the integral represent the interfacial free energy, with f 0 a fourth-order Landau polynomial of the order parameters: h FCC and h BCC in Eq. (1) is the interpolation function for BCC and FCC phases. For multiphase systems, it is introduced in the following format, according to (Moelans 2011): The evolution equations for the mole fractions are derived starting from the following equations, which ensure mass conservation: With x i the mole fraction of component i. As done in many PF models, we make the assumption that the molar volume is independent of the composition and assume V m = 10 − 5 m 3 /mol. It should be noted that considering volume effects in a PF model for solid-solid phase transformations would also require including elastic effects.
According to the KKS model, phase concentration variables x p i are introduced for each solute component in each phase, with p = BCC or FCC, and i = Al, Cr and Fe. The phase specific variables are determined at each position, such that for each of the components, the diffusion potential ∼ µ is equal in the coexisting phases, giving the following set of equations: Where G FCC and G BCC are the Gibbs free energy densities of the FCC and BCC phases, respectively, which have the unit: J/mol. Moreover, the overall mole fractions (1) (2) V m M l for l = i and L li = 0 for l ≠ i, M l the atomic mobilities of the 4 elements and where N represents the dependent component (in this case, Ni). In this work, the atomic mobilities M l are assumed to be independent of composition. Their values were calculated using Thermo-Calc software with the MOBNI4 database and assuming the equilibrium composition of each phase. It was verified that the atomic mobilities are indeed almost constant over the considered simulation concentration region.
To improve the numerical stability of the simulations, the PF equations above were nondimensionalized by introducing the following dimensionless quantities: Where x c , t c , m c, and G c are characteristic scales of the variables, of which the values are given in Table 1. We assume the characteristic energy G c = RT, and the dimensionless diffusion mobility is obtained from Eq. (8), Using Eq.(5), Eq. (4) becomes: The time evolution of the non-conserved order parameters is calculated by the Allen-Cahn equations, which assumes ∂η ∂t = −L η δF δη , giving for the Gibbs free energy of the form of Eq. (1), the equation for the order parameter representing for FCC and BCC phase can be obtained: With the substitution of the dimensionless variables in Eq. (7), Eq. (10) becomes: The relation between simulation parameters and dimensional values are: In this model, boundaries between phases are considered as diffuse interfaces with finite width. Therefore, the model parameters m and k are related to the interfacial free energy γ and width of the diffuse interface profile l as: from which the model parameters m and k are then obtained as: In the considered model, the width of the diffuse interface profile is considered as a numerical parameter and not related to the physical interface width. It is chosen such that there are 8-10 grid points (Moelans et al., 2008, b) to resolve the diffuse transition in values of the order parameters at interfaces. In this paper, the interface thickness is chosen to be ten grids (10 nm), and this thickness will not influence the simulation (11) results. The kinetic coefficient L η should be taken according to Eq. (15) to obtain a diffusion-controlled interface migration in the PF model, as derived in (Moelans 2011) where numerically evaluated values for I η and g are chosen according to Ref. (Moelans 2008), the values used in our simulations are listed in Table 1.
The molar factions x i and order parameters η p are obtained at every time step as the solutions of the diffusion and PF equations. The parameters used for 1-Dimensional (1-D) simulations and the characteristic values considered for the non-dimensional parameters are listed in Table 1 The parameters used for the 2D simulations were chosen in the same way as for the 1D simulation and the corresponding parameters ∼ Lη used for 2D simulations are shown in Table 2. The kinetic Eqs. (9) and (11) and the auxiliary Eqs. (5) and (6) were solved using the finite element method (FEM) in Multiphysics Object-Oriented Simulation Environment (MOOSE) framework (Tonks, et al. 2012). In our simulation, the time step is adjusted based on the number of iterations. The optimal iteration controls the number of nonlinear iterations per time step is chosen as 7. The growth_factor and cutback_factor for adjusting the time step are set as 1.5 and 0.8, respectively.

Fitting composition-dependent Gibbs free energy densities data
As stated in the introduction, this paper aims to show that polynomial fitting of the Gibbs energies can be accurate and effective for MEAs, when only a limited composition region has to be considered, as is for example often the case for a diffusion-controlled phase-transformation and for coarsening. The Gibbs energies for BCC and FCC phases in this paper are fitted as a function of composition using polynomial functions, and they can be evaluated efficiently in the PF method. We found that it is impossible to accurately fit the composition-dependent Gibbs free energy densities over the full composition range using simple polynomials for the quaternary system; however, accurate simulation results are obtained for two-phase systems when the composition-dependent Gibbs free energy densities data are fitted considering a limited composition region in which the two-phase region is included. In this work, the data of Gibbs free energy is sampled with Thermo-Calc TC-Toolbox for MATLAB and using the TCHEA2 database. First, the relevant two-phase region has to be selected as described in section "Determination of FCC/BCC two-phase regiono", and then the Gibbs free energy densities data are fitted over the selected composition domain, as discussed in section "Fitting of Gibbs free energy densities".

Determination of FCC/BCC two-phase region
This paper aims to show that accurate and efficient PF simulation results can be performed for the quaternary system using a polynomial fit of the composition-dependent Gibbs free energy densities data within the two-phase region of interest. The first step in this procedure is to determine the two-phase region based on the phase diagrams computed using the TCHEA2 database. Since it is difficult to visualize and analyze the phase diagram of a quaternary system, the pseudo-ternary isothermal sections are calculated for different Al concentrations. Figure 1 shows the phase diagram sections of (Al) x CrFeNi system at 700 °C with xAl= 0.01, 0.03, 0.05, 0.07(mole fraction). From these diagrams, the relevant range can be determined for each element. The two-phase BCC/ FCC region is shown in green in Fig. 1(a)-(d). When the Al concentration is gradually increased from 0.01 to 0.07, the BCC/FCC two-phase region decreases. From these diagrams, it is decided to fit the Gibbs free energy densities polynomials over the composition range 0.01 < xCr < 0.97, 0.01 < xAl < 0.08 and 0.01 < xFe < 0.5, which is slightly larger than the region in the phase diagram to ensure the two-phase region is fully included.
The corresponding composition range over which the Gibbs free energy densities data will be fitted and is listed in Table 3. 2D simulation results for the quaternary system are also compared with simulation results for the ternary CrFeNi system to test the capability of the PF model and study the phase transformation behavior with different alloy systems. A similar fitting procedure was used for this ternary CrFeNi system, although it is much easier because it is straightforward to determine the relevant two-phase region. The fitting details for the ternary system are given in the Additional file 1.

Fitting of Gibbs free energy densities
First, the Gibbs free energy densities were sampled over the two-phase composition range (given in Table 3) using Thermo-Calc toolbox in combination with the TCHEA2 database. A step size δ xCr = δ xNi = δ xFe = 0.001 is used. For the fitting process, the polynomial's optimal order is determined, as the one for which the deviation between the sampled data and those obtained by evaluation of the polynomial fit is minimized. For example, it is checked that the R-square value, which represents the square of the correlation between the response values and the predicted response values, is close to one. Fig. 2 shows the scatter plot of the absolute error △G as a function of G, for the case (b) and (d) with the ideal mixing term ( RT i x i lnx i ) included in G and the case (a) and (c) for which the ideal mixing contribution is subtracted before the fitting for the fourth-order polynomial. △G is the difference between data extracted and fitted. Results for BCC phase show that the absolute error between the sampled and fitted data is within the range of [− 30,50] when the ideal mixing term is not considered and [− 150, 200] when the ideal mixing term considered in the data. For FCC phase, the absolute error between the sampled and fitted data is within the range of [− 60,20] without the ideal mixing term and [− 150, 200] with the ideal mixing term. The ideal mixing term Zuo et al. Materials Theory (2022)  increased the fitting error because of the difficulty in fitting the logarithmic behavior using polynomial functions. Therefore, since the form of the ideal mixing term is known, the ideal mixing contribution is subtracted in the preprocessing step and added back to the Gibbs free energy when conducting the simulation since this contribution is not system or phase-specific. As shown in Fig. 2 (a) and (c), △G is near zero; the relative error(△G/G*100%) is small (less than 0.12% for FCC and 0.0656% for BCC). The probability density estimate for Fig. 2 Fig. 2 (e) and (f ), respectively. Results show that most data are near zero and located in the range [− 20, 20]. Figure 3 shows the comparison of △G as a function of G between different polynomial orders for both BCC and FCC phases when the ideal mixing term is subtracted. It is shown that the △G for the third-order polynomial is within the range of [− 100, 100] for BCC and [− 200, 200] for FCC, which is significantly larger than that for the fourthorder and fifth-order polynomials. The majority of data between the fourth-order and fifth-order polynomials for the FCC phase are located within the range of [− 20, 20]. Such difference is not significant. The △G for the fifth-order polynomial of the BCC phase is decreased compared with the fourth-order. The majority of data for the fifthorder is located within the range [− 10, 10] for BCC phase. A fourth-order polynomial is chosen to fit the Gibbs free energy densities data for the considered system and composition range. The details on the fitting accuracy are validated by comparing the simulation results between the PF model and DICTRA. But it is important to note that this is not the only option (e.g., the fifth-order for BCC and fourth-order for FCC is also a suggested choice). For other systems or composition ranges, the optimal order of the polynomial can differ, or a different step size in the data collection may be more appropriate.

Validation by comparison with 1-D sharp Interface calculations (DICTRA)
To validate the obtained thermodynamic and kinetic information and corresponding parameters for PF simulation. Sharp interface simulations using DICTRA and the same Gibbs free energy densities and diffusion mobility databases were carried out for the same conditions as the PF simulations. Also in DICTRA, a constant molar volume of 10 − 5 [m 3 /mol] is assumed for both phases. The results are compared and discussed in this section for the quaternary AlCrFeNi alloy. A similar validation was performed for the ternary CrFeNi alloy and is discussed in the Additional file 1. The 1-D simulations are performed at 700 °C. The overall composition of the system was arbitrarily taken within the BCC/FCC two-phase region. The corresponding equilibrium phase fractions and compositions of the BCC and FCC phases for the given overall composition were calculated using Thermo-Calc software and are listed in Table 4.
The initial phase compositions and volume fractions set for the calculations for each phase are listed in Table 5. In the simulation, the composition of the FCC phase is always taken differently from the equilibrium composition, while the composition of the BCC phase is kept at the equilibrium value (from Table 4). The initial phase fraction for BCC is set as 0.2 and FCC as 0.8, such that the overall composition is still the same as the one given in Table 4. In the simulations, the composition and phase fractions will evolve towards the equilibrium values listed in Table 4. The simulation results could be compared with values in Table 4 to see whether the simulation can finally reach the same equilibrium state. = 4 (f) order = 5 Fig. 4 (d). The interdiffusion calculations with the same initial conditions is also carried out using the PF method. The concentration profiles at t = 1, 20, 30 and 100 days for PF method are shown in Fig. 4(a1),(b1),(c1). The evolution of the BCC/FCC interface position is shown in Fig. 4 (d1). Results show that, in general, the PF method calculations reproduce the results from DICTRA simulations. For a more quantitative determination of the error, the relative errors (absolute error/ measurement being taken:[|a-b|/a]) for the local mole fractions of the elements between the PF method and DICTRA simulation during the diffusion process at different times are shown in Fig. 5 (a1,b1,c1) for 1 day, (a2,b2,c2) for 20 days, (a3,b3,c3) for 30 days and (a4,b4,c4) for 100 days. In the PF model, the interface is diffuse, while in DICTRA a sharp interface is assumed. The relative errors over several grid points in the interfacial area are therefore large; however, this deviation at the interface has no considerable effect on the bulk composition and position of the interface. In the bulk, the relative errors on the mole fractions are around zero. These findings show that the proposed method to incorporate the fitted Gibbs free energy densities in PF simulations for quaternary systems gives accurate results that are consistent with results obtained using DICTRA simulation software using the full Gibbs free energy densities description.
In order to quantify the error between PF and DICTRA simulations along the diffusion time, the average relative error for each composition in the bulk phase is calculated and shown in Fig. 6(a). The relative errors on the results for the mole fraction of Al, Cr, and Fe at 1 day are less than 0.2, 0.15, and 0.05, respectively. Moreover, the relative error decreases with simulation time. For the simulations of 100 days of diffusion, the relative error for all elements drops below 0.05 at 100 days. The comparison of results for the interface position is shown in Fig. 6(b). We find that the interface position shows good agreement between PF and DICTRA simulations. The insert in Fig. 6(b) shows the maximum difference for the interface position is less than one grid point (The size of a grid point Δx was taken equal to 1 nm, as listed in Table 1).

Growth and coarsening of BCC precipitates
In order to compare the results between different alloy systems, validation for ternary CrFeNi alloy is also conducted. The details for the validation process of CrFeNi alloy are shown in the Additional file 1. After the PF model was validated by comparison with 1-D DICTRA simulations, it was applied to study the growth and coarsening phenomena of BCC precipitates inside the FCC matrix. Three alloy systems were analyzed to study the effect of different concentrations on the BCC/FCC transformations, and the initial concentration settings for different alloys are listed in Table 6. For alloy 1, there is no Al in the system. For alloy 2, an Al mole fraction of 0.04 is added to the system of alloy 1. Compared with alloy 2, alloy 3 has a larger amount of Fe. For each alloy, two cases are considered. In Case 1, the initial compositions of the FCC and BCC phases are taken equal to the equilibrium compositions of the phases (thus, only Ostwald ripening is expected to occur). For Case 2, the initial composition in the FCC phase deviates from the equilibrium composition (therefore, both precipitate growth and Ostwald ripening can be expected). All simulations have the same initial geometrical configuration (Dimension: 100 nm × 100 nm with 100 × 100 grid points), as shown in Fig. 7(a).
For all cases, the small-sized precipitates gradually dissolve into the matrix, and the large-sized precipitate grows with time. The Ostwald ripening is well represented in the simulations. Since all alloys have a similar growth behavior, only the quaternary alloy 2 -case 2 is further discussed in detail. Figure 7 shows how the precipitates evolve with time for alloy 2 -case 2 with an initially supersaturated FCC matrix. The These results show that the mole fraction of Cr in both BCC and FCC decreases during the diffusion process while the mole fraction of Fe and Al increases in both phases. At the early stage, the concentration profile is not symmetrical during the diffusion process due to the influence of neighboring precipitates, and it will become symmetrical again when all smaller precipitates are dissolved. Figure 8 shows the evolution of the BCC precipitates area as a function of time for the three alloy systems and Cases 1 (alloy 1-case 1, alloy 2 -case 1, and alloy 3 -case 1) and 2 (alloy 1 -case 2, alloy 2 -case 2, and alloy 3 -case 2), for which the FCC phase is respectively initially in equilibrium and initially supersaturated. It is evident that the precipitates grow faster and the final size of precipitate area is larger when the FCC phase is initially supersaturated (case 2) than for the equilibrium condition (case 1). For all cases, due to the curvature effect, the precipitates shrink at the beginning and then grow at different speeds due to different initial conditions. The PF model is thus able to study coarsening processes for this quaternary alloy system.
Results also show that without Al in the system (alloy 1), the precipitates in the ternary system grow faster. Compared with alloy 2, alloy 3 has a higher amount of Fe, but a lower amount of Ni in the system and the growth of precipitates in alloy 3 is the slowest. The effect of changes in the amounts of different components in the system on the growth rate of BCC precipitates is evident, indicating the importance of having full thermodynamic and kinetic information for simulating real HEAs/MEAs material.
The growth rate in the multicomponent diffusion process depends on the supersaturation conditions and thermodynamic quantities in the system. The gradient of diffusion potential as the driving force for the movement of atoms in the alloy system Table 6 Initial composition for BCC and FCC phases for the three different alloys considered in the simulations and assuming equilibrium(EQ) and nonequilibrium(NEQ) composition of the matrix as initial setting of the simulation  The diffusion potentials of case 2 among alloy 1, alloy 2, and alloy 3 are shown in the Additional file 1. The corresponding gradient of diffusion potentials and mobilities are shown in Fig. 9 and Fig. 10, respectively. For the different alloys, results at time steps close to each other were taken to make a comparison between these cases possible. At the chosen time, coarsening, i.e., the growth of the larger particles at the expense of the smaller ones, is occurring. The color bar is set as the same range to compare the gradient of diffusion potentials for alloy 1, alloy 2, and alloy 3. The gradient of diffusion potentials is more significant at the phase boundaries. The difference becomes smaller with time, indicating that the microstructure evolves towards equilibrium. The red color shown in Fig. 9 represents the maximum gradient of diffusion potential, which shows that such a value (around 0.3) for Cr, Al, and Fe appears in alloy 2 and alloy 3.
The number of independent interdiffusion mobilities necessary to describe multicomponent diffusion for the n component system is n(n-1)/2. The color bar is set the same for each alloy. Since the diffusion in the matrix phase determines the kinetics of the growth of the precipitates, the average values for each mobility in the FCC phase are compared.  The different growth behavior between alloy 2 and alloy 3 (same alloy system) shows that increasing the amount of Fe (with reducing the dependent amount of Ni) reduces the growth rate of BCC. The value of M CrFe is negative and the absolute value is twice as large as in alloy 3, slowing down the diffusion process. The M CrCr and M FeFe are larger in alloy2 while M AlAl is slightly larger in alloy3, which indicates that the precipitate growth rate of each alloy is the combined result of both gradient diffusion potential and all the dependent mobilities.

Discussions
Multicomponent PF models were developed for many years; however, there have been few reports on PF models for actual multicomponent material diffusion. One of the reasons for the paucity of the phase-field modeling in a multicomponent alloy is the complexity of coupling the chemical free energy into the multicomponent PF model. For the alloys studied in the paper, the thermodynamic database is not yet publicly available. This paper introduces a relatively straightforward way to couple the thermodynamic and kinetic data into the PF model to study the phase transformation in MEAs/HEAs.
-As this work aims to deal with the coupling of thermodynamic and kinetic information into the PF model, for the sake of simplicity, the elastic effect on the kinetics and morphology of the precipitates is neglected, the partial molar volumes of the components are set as 10 − 5 [m 3 /mol]. The stresses arising from compositional strains due to the differences in the partial molar volumes of the components could change the stress and strain tensors, thus changing the microstructure of the alloys. Such elastic effects will be discussed in our future work. -There is no generally applicable optimal polynomial order for Gibbs energy fitting.
The optimal polynomial orders can vary among systems and depending on the desired accuracy. The fitting accuracy is affected by the data collection steps, fitting regions, and the polynomial orders. The accuracy can be evaluated in two ways: (1) the deviation between the sampled data and those obtained from the polynomial fitting△G is small, as shown in Fig. 2 (the R-square value is reaching or equal to one). Thus, the relative error is small (less than 0.12% for FCC and 0.0656% for BCC).
(2) the PF simulations could reproduce the equilibrium phase fractions and compositions as calculated through Thermo-calc software. (If the initial condition for the system is set as equilibrium, then the phase transformation would not happen in PF simulation. If the initial condition is set as nonequilibrium, the system will finally reach the equilibrium condition, as shown in Fig. 4). In our simulation, the fourthorder polynomials are chosen for the simulation and found it gives accurate results. Thus it is considered the fourth-order polynomials are appropriate for our system. -The kinetic data is considered to be accurately coupled if the PF model can reproduce the diffusion path in DICTRA simulations with the same initial conditions. The partial molar volume of each constituent in our PF model is assumed to be 10 − 5 [m 3 / mol], equal to that in DICTRA simulations.

Conclusions
A PF approach is presented to simulate the diffusion-controlled phase transformation between FCC and BCC in (Al)CrFeNi alloys. Comparison of 1-D PF simulations with results from DICTRA calculations shows that the composition dependence of the Gibbs free energy densities and diffusion mobilities for this multicomponent system is considered accurate and effective with the proposed methodology. The current approach has thus the potential to be applied for multicomponent alloy design, such as in the optimization of HEAs and MEAs.