A first-principles phase field method for quantitatively predicting multi-composition phase separation without thermodynamic empirical parameter

To design tailored materials, it is highly desirable to predict microstructures of alloys without empirical parameter. Phase field models (PFMs) rely on parameters adjusted to match experimental information, while first-principles methods cannot directly treat the typical length scale of 10 μm. Combining density functional theory, cluster expansion theory and potential renormalization theory, we derive the free energy as a function of compositions and construct a parameter-free PFM, which can predict microstructures in high-temperature regions of alloy phase diagrams. Applying this method to Ni-Al alloys at 1027 °C, we succeed in reproducing evolution of microstructures as a function of only compositions without thermodynamic empirical parameter. The resulting patterns including cuboidal shaped precipitations are in excellent agreement with the experimental microstructures in each region of the Ni-Al phase diagram. Our method is in principle applicable to any kind of alloys as a reliable theoretical tool to predict microstructures of new materials.

M icrostructures involving precipitations, grain boundaries, dislocations, and other defects play a decisive role in many important properties of alloys such as ductility, plasticity, toughness, magnetism, oxidation-and heat-resistances. The quest of new materials with desirable properties requires microstructure engineering of materials by changing composition, morphology, pressure, and temperature as well as doping, casting, and forging. Since there are many degrees of freedom, it is highly desirable to find fundamental breakthrough toward the design of tailored materials. For this purpose, powerful computational techniques are required to predict microstructures. Since the length scale of microstructures is typically 1-100 μm and the number of atoms involved is 10 11 -10 17 , first-principles methods such as density functional theory (DFT) 1 are not applicable. Phase field models (PFMs) 2-4 offer a promising computational tool to study such phenomena, where microstructures are described by order parameters. However, PFMs are purely empirical and adjustable parameters must be used so that the computational output matches the experimental one. Therefore, they are not widely accepted in industries 5 .
In this paper we first introduce how to combine first-principles method and PFM, and solve the previous problems of PFMs relying on parameters. Although, a number of theoretical studies have been conducted using atomistic methods [6][7][8] , they are limited, for example, to the level of the cluster variation method (CVM) 9 , which relies on a particular crystal lattice and a particular alloy composition. A CVM-based phase field method has been also developed to study microstructure evolution in alloys 10,11 . In this method the alloy composition was kept fixed and the order-disorder phase transformation was treated at different temperatures. However, this method cannot discuss the entire phase diagram. Here, we propose a first-principles and non-parameter based PFM. The method can clearly discriminate the local composition and evolution mechanism of microstructures of alloys without any material parameter.
Although our method is in principle applicable to any kind of alloys, we demonstrate its ability by treating Ni-Al binary alloys as an example, which have attracted considerable attention for their excellent mechanical properties; very hard and good oxidation-and heat-resistances suitable for turbine disks and blades 12,13 . With varying Ni and Al compositions, they undergo many phase transformations. We will reproduce the experimental phase diagram 14 and show the time evolution of microstructures. We will focus on the ordered phases only, because the description of liquid and disordered phases requires to introduce another non-conserved order parameter.
In conventional PFMs, the order parameters are described by continuous functions of space and time. The free energy is defined as a polynomial of these order parameters and therefore are continuous too. These kind of representations do not take into account of the local structures and compositions, which are very important to find phases and microstructures in materials. By decomposing the space into a fine regular 3D grid of unit cells containing tetrahedron or other polygon (Fig. 1) and then taking the continuous limit where the lattice constant goes to zero, one can consider number density as the order parameter, which will include the atomic arrangement and local composition and eventually can construct a more effective free energy. In case of Ni-Al alloys, the four atomic sites of each tetrahedron can be filled with either Ni, Al, or vacancy, as shown in Fig. 1b. Thus, instead of continuous order parameter, one will obtain integer functions defined as the local composition of Ni n Al m : φ Ni , φ Al = 0-1, 1-2, 2-3, 3-4, 4-5 for n, m = 0, 1, 2, 3, 4, and use cluster expansion theory [15][16][17] to determine the local energy from ab initio DFT. To include off-lattice effects, the potential renormalization theory can be implemented, where the atomic displacement is renormalized in a length scale shorter than the lattice constant 18 . This theory is applicable at high temperatures where Einstein model is valid. Using this discretized free energy definition, together with the potential renormalization theory for the temperature effect, we simulate the evolution of microstructures and phases in Ni-Al alloy systems at various compositions. We fix the temperature at T = 1027°C (=1300 K), which is typical in jet-engine turbines. Our results are of great match with the experimental and conventional phase field findings.

Results
The first-principles free energy and the diffusion equation. The resulting local free energies F(φ Ni , φ Al ) are shown in Figs. 2 and 3, and are plotted in 1D and 2D in Fig. 4a, b. Each block (plateau) in the 2D (1D) plot corresponds to an integer (n, m) composition of the Ni n Al m alloy, including vacancies for 0 ≤ (n + m) ≤ 4 and interstitial atoms for (n + m) > 4. For n + m = 5, the most stable trigonal unit cell 19 is chosen with the same volume as the cubic unit cell. The explicit functional form for the free energy can be represented as: ð1Þ and E Ni n Al m is the energy/ renormalized energy of the Ni n Al m cluster. The free energy is non-negative when (n + m) > 6 and increases sharply, which is described by a polynomial for (n + m) < 0 and (n + m) > 6 so as to avoid unrealistic compositions. The chemical potentials are given by where ε X (X = Ni or Al) denotes the gradient energy coefficient of φ X , and ε X ∇ 2 φ X represents the interface energy contribution. From the continuity equation ∂φ X /∂t = −∇ · J X and the flux introduced by J X = −M X ∇μ X , the generalized Cahn-Hilliard equation is derived as Cahn-Hilliard equation is applicable to conserved order parameters such as the composition and the atomic density, while for non-conserved order parameters, Allen-Cahn equation should be used in the phase field model. We assume that M X and ε X are independent of the species X. Then, M X and ε X can be set arbitrary, e.g., at 0.00125 and 0.5, respectively, by rescaling time with M X and length with ε X . Unlike atomistic simulations, in our coarse-grained model, there is only the change in the concentration and in this case to determine the time and length scales becomes particularly difficult. In order to avoid this difficulty, we assume simulation time and simulation cell size to be arbitrary units, which can be scaled to the experimental time and length if required. To obtain an exact relation between the simulation time (cell size) and the experimental time (length), one need to develop a method to calculate the mobility (interface energy) from firstprinciples and give this as an input parameter. But this generally requires a huge computation, and the method we are proposing here is to simply avoid this difficulty. Because the free energy is replaced with its local values, it is basically necessary to include the random force in the phase field equation. However, this does not affect the final pattern much. We have chosen the amplitude of the random force as 0.5 for all the calculations. For simplicity, we use a 2D model in a grid space of 124 × 124. The grid space  Fig. 1 Pictorial representation of cluster expansion and potential renormalization theory. a The real space is divided into a fine regular 3D grid. One grid cube represents one unit cell and different colors represent different alloy composition within the unit cell. One unit cell of the Ni-Al alloy containing a tetrahedron formed by four nearest neighbor atomic positions for b fcc lattice, c bcc lattice. d Unit cell of trigonal structure. These atomic positions are occupied by n Ni atoms (yellow) and m Al atoms (orange) representing the composition Ni n Al m , in the cluster expansion theory. The 3D grids around each of the atoms are the points, where the atom was displaced during potential renormalization calculation  Microstructure evolution at various alloy compositions. Here, we restrict ourselves to the solid-state phases of Ni n Al 4−n with n > 0.88 (Ni > 22%) as shown in Fig. 5a. The initial structure contains Al-and Ni-rich seeds in a uniform matrix as shown in Fig. 5b. The resulting structures (φ Ni ) are plotted in Fig. 5c Fig. 5c-f. For small Ni concentrations, the Ni-rich seeds start dissolving followed by Al-rich seeds forming a homogeneous solid solution with very small variation (Fig. 5c). In the middle of this region, Al 4 precipitates are obtained within a matrix formed by mixed Ni 2 Al 3 and NiAl 3 . These precipitates were random with dendritic signatures for Ni 25.5% and rectangular for Ni 30% as shown in Fig. 5d, e. For snapshots of the simulation, see Fig. 6a. The Al-rich seeds grow and Ni-rich seeds dissolve with time, showing Ostwald ripening. The particle size decreases with increasing Ni concentration and the microstructure becomes almost homogeneous for Ni 35.5% ( Fig. 5f) as we enter region II.
In the next region II ranging from Ni 38 to 40%, the initial structure disappears (Fig. 5g for Ni 40%) by Ni-rich seeds first getting dissolved followed by the Al-rich seeds. Figure 6b shows the snapshots. This region is the single phase region of Ni 2 Al 3 , where the total composition is n + m = 5.
In region III also, we obtained homogeneous pattern with some slight variation as shown in Fig. 5h for Ni 41%. The snapshots in Fig. 6c shows similar time evolution as Ni 40%.
Region IV is the widest among all solid phase regions. We observe uniform solid solution (Fig. 5i, j for Ni 47 and 52.5%) in this region except near the right boundary. Near the right boundary, Ni-rich seeds precipitate. First, circular shaped Ni 3 Al particles are formed in a matrix of Ni 2 Al 2 , which transformed  into square shapes at t = 10 4 as shown in Fig. 6d for Ni 60%. The particle size grows and coalesced particles are formed by some particles at t = 6 × 10 4 along with the formation of some smaller particles in the matrix. With increasing Ni concentration, the particle size increases. The microstructure composition changes in region V ranging from Ni 62 to 73% (Fig. 5l, m for Ni 65 and 70%). Here, Al-rich seeds precipitate as shown in Fig. 6e for Ni 65%. Spherical Ni 2 Al 2 particles are formed within a uniform matrix of Ni 3 Al at t = 4 × 10 3 . They transform into rectangular shapes along with the formation of some coalesced ones (t = 4 × 10 4 ). The particle size increases as the time proceeds (t = 10 5 ). The particle size decreases with increasing Ni concentration and distribution becomes homogeneous in the next region. Region VI is a single phase having no microstructure as shown in Fig. 5n for Ni 75%. The corresponding snapshots are shown in Fig. 6f.
In region VII, rectangular Ni particles are formed in the Ni 3 Al matrix from the Ni-rich seeds of the input structure. Internal structures appear in these particles as shown in Fig. 6g for Ni 78% at t = 2 × 10 3 . These particles grow further at later times to form square shapes composed of Ni 4 only (t = 8 × 10 4 ). The particle size further increases and rectangular particles with various sizes are formed at t = 1.4 × 10 5 (Fig. 5o). With increasing Ni concentration, the particle size increases and some coalesced particles are formed (Fig. 5p). Some smaller particles with irregular shapes appear for Ni 80% (Fig. 5q). For Ni concentration higher than 80%, Al-rich seeds grow forming spherical shapes, that transform into rectangular shape as shown in Figs. 5r and 6h for Ni 82%. The resulting microstructure consists of Ni 3 Al particles embedded in pure Ni 4 matrix. For higher Ni concentrations, the particle size reduces and the microstructure disappears forming a uniform phase in the next region; see Figs. 5t and 6i for Ni 92%.

Disscussion
To understand the growth mechanism of the microstructures from the initial pattern, we plot the change in local concentration and hence the local free energy at various time steps, as shown in Supplementary Fig. 1a, b for Ni 60%. We also plot the magnitude of the free energy gradient (|∇F|), which corresponds to the local stress of the system at each grid point in Supplementary Fig. 1c. As expected the local stress is non-zero only at the interface. The interfacial thickness is related to the (visible) non-zero gradient region. These plots show a sharp interface and we can expect bulk region even in the smallest precipitates of the microstructures. The plot of the total stress of the system versus the simulation time step shows an initial increase in stress until a maxima is    0  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100   0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100   0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100   0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100   0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100   0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80 100  0  20  40  60  80  reached. After this increase, the stress decreases slightly and becomes constant after sufficient simulation time (see Supplementary Fig. 1d). This initial increase in the stress corresponds to the diffusion of the individual species. After the system forms a stable microstructure, the stress also becomes constant. In addition to the growth mechanism, this inherent local stress in particular at the phase boundaries influences the shape of the precipitates too and we have obtained cuboidal precipitates, which are peculiar feature of the NiAl alloy without explicitly introducing any external parameter for coherency stresses and/or strain. In atomic scale, the lattice mismatch is the basic origin of the coherency strain. The difference in the lattice parameter can be included in our PFM as vacancies and interstitials as well as the mixture of the two phases. This effect can be included in our PFM in which the free energy is a discrete variable for integer compositions. This is not possible in all previous PFMs, in which the free energy is a continuous variable with respect to compositions. This is a distinct advantage of the present PFM. The most important point is that, in our PFM, the density is automatically fine tuned by the clusters with vacancies and interstitials as well as the mixture of the phases. Therefore, the detailed lattice mismatch can be handled without introducing anisotropic elastic energies. Thus, we can conclude that the main origin of the cuboidal precipitates is that our free energy is a discrete variable. Indeed, if we change the definition of the free energy to be a continuous variable with respect to compositions, we obtain round precipitates. The resulting microstructures and phases from our phase field method are listed in Table 1 for each region in the solid phase diagram and compared with the experimental and/or empirical phase field model results. In region VII, we show inversion of microstructure from γ-Ni precipitates in γ′-Ni 3 Al matrix (Fig. 5o-q), to γ′-Ni 3 Al precipitates in γ matrix (Fig. 5r, s). However, there is not yet any experimental demonstration of microstructures having Ni precipitates in Ni 3 Al matrix, commonly termed as inverse superalloy. However, there is indication of possibility of such microstructure in some experimental papers as these microstructures are expected to have better hardness compared to the normal alloys (γ′ precipitates in γ matrix). Vogel et al. 20 showed γ-Ni particles in the γ′-Ni 3 Al precipitates in hierarchical microstructures of Ni-Al alloy. A very recent paper 21 established an analogy between the inverse alloy and the hieirarchical Ni 86.1 Al 8.5 Ti 5.4 alloy. Therefore, we believe that the prediction of such an inverse alloy will be seen in experiments too.
Our results are everywhere in excellent agreement with the previously reported microstructures. We compare the rectangular particles in Fig. 5r with that of Fig. 1 (heat-treated at 1350°C) of ref. 22 , both for Ni 82%, in terms of its size, and estimated our length scale as Δx ≅ 0.03 μm. With this value we estimate ε X to be 8.9 × 10 −11 Jm −1 (see Supplementary Note 1), which is comparable to the gradient energy coefficient used in the previous phase field calculation on Ni-Al alloys 23 .
We have assumed fcc lattice for various compositions of the Ni-Al alloy system for our phase field calculation. The exceptions were for n + m = 5 cases and Ni 2 Al 2 composition. For all the n + m = 5 compositions, we have used trigonal unit cell (similar to Ni 2 Al 3 structure). Ni 2 Al 2 on the other hand is known to have a bcc (CsCl) structure. For a more accurate calculation, we have modified the free energy by constructing a tetragonal cluster from the bcc lattice point as described by Allen et al. 24 . Unlike fcc lattice, the cluster is formed by two adjacent unit cells as shown in Fig. 3. We kept the volume of the unit cell containing one tetragonal cluster same as that of the fcc cluster such that the number density is preserved. For the free energy calculation, we constructed a 2 × 2 × 2 supercell for each Ni-Al combination. The structures, total energy per tetragonal cluster, renormalization correction and the renormalized energies are shown in Fig. 3 for each no vacancy Ni-Al combination, respectively. The dotted plateaus in Fig. 4a denote the free energies for the clusters in a bcc lattice. As expected, the bcc energy is lower than the fcc energy for the Ni 2 Al 2 composition only. We first repeated the firstprinciples phase field calculation for Ni 47, 52.5, and 60% using this free energy, and confirmed that the resulting microstructures are similar to the ones obtained with fcc lattice. Next, we replaced the fcc energy with the bcc energy only for the Ni 2 Al 2 composition, and repeated more realistic phase field simulations shown here. In all regions of the phase diagram, the resulting microstructures using the bcc Ni 2 Al 2 energy perfectly coincide with those using the fcc Ni 2 Al 2 energy. This suggests the validity of using the fcc energy only. This is due to the fact that the difference between these two energies was occasionally small. Of course, it is better to choose the most favorable lattice having the lowest energy for each Ni n Al m composition as in the present simulation. This example nicely demonstrates the validity of this treatment. The most important point is that the lattice (unit cell) can be different in different composition. Our method is a universal one and can be used for any alloys in any lattice structure.
Since a 3D model with those parameters determined by firstprinciples calculations can be used for constructing a more realistic model, we performed 3D simulations for some of the alloy compositions by taking a smaller system size of 40 × 40 × 40 grid points, using the same simulation parameters as in the 2D simulations. The initial structure consists of two seeds in a uniform matrix of the alloy. The resulting microstructures (shown in Fig. 7) are very similar to the ones obtained in the 2D simulation even for the spatial scale of the resulting patterns. This suggests that the 2D simulations are good enough to reproduce the real microstructures. Using this model, we confirmed that the resulting PFM works well giving the reliable time evolution of microstructures for various compositions with no empirical parameter in the thermodynamic part of the model. Thus, our first-principles phase field method is definitely a successful achievement in determining microstructures of the length scale of 1-100 μm purely from quantum mechanical ab initio theory. We strongly believe that the present method becomes a future theoretical standard for materials, which is fundamentally different from the existing methods limited only within explanation of experimental observations with empirical parameters. We are creating an automatic submission protocol for potential renormalization calculation that can be used to provide all the necessary input data for any system. This method has a potential to predict new useful materials in industries by supplying T-C phase diagram with microstructures, which is crucial to discuss about material properties in realistic applications. Our first-principles PFM is a universal method to perform large scale simulation for variety of materials at less computation time.

Methods
Coarse graining procedure. In this model we have considered the number density φ X of the constituent elements, i.e., X = Ni and Al as the conserved order parameters. We calculate the free energy, F as a functional of these phase field variables by using cluster expansion theory [15][16][17] . It provides a very simple yet powerful approximation to calculate total energy of a system with a large number of substitutional structures. It allows us to calculate thermodynamic properties of a very large system by simplifying it into Ising-type models which deals with much simpler and discretized coarse-grained systems. This method has been widely used in the calculations such as formation energies of random alloys and temperature-composition phase diagrams. For Ni-Al alloys, which are mostly found in an fcc structure, we use the tetrahedral approximation 15,16 . In a tetrahedron cluster (as shown in Figs. 1b, 2, and 3) there are maximum four possible atomic sites to be filled by the two atomic species Ni and Al or left vacant. In cluster expansion theory, the total energy can be expressed as a summation of the product of many-body interaction potentials (J i ) and multisite correlation functions ξ i for the ith order cluster (i = 0, .., 4 for tetrahedron approximation). The sum is over all the ith order clusters for this lattice type. Using DFT, the interaction potentials (J i ) and the total energy for various configurations Ni n Al m (0 ≤ m + n ≤ 4) are calculated by filling the four sites of the tetrahedral unit cell by Ni or Al or vacancy. For Fe-Pt alloys with vacancies, see ref. 25  Effect of temperature. In the total energy calculations, the internal entropy arising by the displacement of the atoms in a length scale shorter than the lattice constant is very important in describing, e.g., nucleation, defects, crack propagation, and crystal growth. Without this effect, lattice models overestimate the order-disorder phase transition temperature. To include the internal entropy effect into our model, we apply potential renormalization theory 18 , in which, the bare interparticle interaction in the off-lattice system is renormalized so as to give the same local partition function at a given temperature. This is done by discretizing the configuration space (Fig. 1b) and then taking the trace of the local movement of the atom from the lattice point for each of the configurations as follows: Here, ΔF denotes the change in the free energy due to the movement of the atom, ΔE jk is the increase in energy when one atom of the species j = Ni or Al is moved to r i , and C j is the concentration ratio. For normalization correction, Misumi's procedure [18] of dividing by 8 is inappropriate and a more suitable procedure is to divide by 4. The details about the potential renormalization calculation can be found in the Supplementary Note 2. Adding ΔF calculated in this way to the on-site energy, we obtain the required local free energy for each of the compositions shown in Figs. 2 and 3. The energies are renormalized at T = 1300 K (1027°C). This procedure is valid at high temperatures, where the Einstein model becomes applicable. To partially incorporate the magnetism of Ni in the free energy, we performed spin polarized energy and renormalization (for non vacancy clusters) calculations for all the clusters containing Ni atom in VASP. The 1D and 2D plots of the free energy calculated in this way as a function of φ Ni and φ Al are shown in Fig. 2.
Free energy functional. In the previous PFMs, the free energy functional is expressed by a very simple polynomial of a continuous field variable 27 . In our PFM, the field variables φ X (X = Ni or Al) have numbers varying from 0 to 5 such that φ X = 1~2 corresponds to one site in the tetrahedron being occupied by X atom, φ X = 2~3 corresponds to two sites in the tetrahedron being occupied by X atom, and so on. Because of this definition of φ X , the condition for all sites of the tetrahedron to be filled by either Ni or Al is given by φ Ni + φ Al = 5. This is equivalent to the condition in terms of atomic fraction, i.e., n + m = 4 for the Ni n Al m alloy. For a uniform definition, we normalize φ X to x of Ni x Al 4−x as . This relation roughly produces x ¼ 4 5 φ X for x < 1.0 and x > 3.8, and x = φ X − 0.5 for 1.0 < x < 3.8. Then our results can be compared with the experimental results, with any mixing composition of Ni and Al.
Phase field simulation. We have used this free energy, to perform phase field simulation at various regions in the phase diagram along the temperature line as shown in Fig. 3a. For the effect of temperature on the microstructures please see Supplementary Fig. 2 (for Ni 60% composition). For each of the global compositions Ni n Al m , we define the initial phase field densities, φ Ni and φ Al by giving a constant value each corresponding to n and m, respectively, calculated using the expression relating x to φ X . In this uniform matrix we introduce some random circular seeds, half of them having composition Ni n+c Al m−c and the remaining half as Ni n−c Al m+c , so that the total composition of the complete system remains the same as Ni n Al m as shown in Fig. 3b for Ni 30% alloy. Here, c is a very small constant parameter in our simulation (typically around 0.3). Defining this initial pattern, we perform the phase field simulation as per the above method. It is very important to input an initial pattern, without which the microstructures will not be formed. This is because a homogeneous pattern is a trivial solution to Eq. (3) corresponding to a stable distribution. To avoid the system going to such local minimum position, we need to assign some initial fluctuation in the input structure for example, by distributing random initial seeds. The amplitude c of these seeds is of more importance than the number and distribution within the matrix. We observed that if c < 0.3, the system becomes homogeneous. This result strongly suggests the nucleation growth mechanism, which coincides with the experimental evidence of NiAl alloys as discussed by R. Moskovic 28 .

Data availability
All data generated and analysed during this study are included in the published article and its Supplementary Information.

Code availability
Code that supports the findings of this study is available in the published article as Supplementary