Critical Evaluation and Thermodynamic Re-Optimization of the Si–P and Si–Fe–P Systems

Abstract

Thermodynamic modeling of the Si–P and Si–Fe–P systems was performed using the CALculation of PHAse Diagram (CALPHAD) method based on critical evaluation of available experimental data in the literature. The liquid and solid solutions were described using the Modified Quasichemical Model accounting for the short-range ordering and Compound Energy Formalism considering the crystallographic structure, respectively. In the present study, the phase boundaries for the liquidus and solid Si phases of the Si–P system were reoptimized. Furthermore, the Gibbs energies of the liquid solution, (Fe)₃(P,Si)₁, (Fe)₂(P,Si)₁, and (Fe)₁(P,Si)₁ solid solutions and FeSi₄P₄ compound were carefully determined to resolve the discrepancies in previously assessed vertical sections, isothermal sections of phase diagrams, and liquid surface projection of the Si–Fe–P system. These thermodynamic data are of great necessity for a sound description of the entire Si–Fe–P system. The optimized model parameters from the present study can be used to predict any unexplored phase diagrams and thermodynamic properties within the Si–Fe–P alloys.

1. Introduction

Si has been commonly recognized as the prime candidate for solar industry applications owing to its excellent photoconductive and electrical properties [1]. However, the photovoltaic conversion efficiency and electrical conductivity of Si solar cells highly depend on the level of impurities. In order to fabricate solar-grade Si with high purity, relatively inexpensive metallurgical-grade Si and ferrosilicon alloy are often selected as raw materials for Si refinement [2]. One of the dominant challenges of such a process is to eliminate P, which needs to be controlled to as low as $1\times {10}^{-5}$ wt.% to meet the requirements for solar cells [3]. Currently, several methods, including vacuum refining [4,5], directional solidification [5,6], slag refining [7,8,9], electron-beam melting [10,11], and solvent refining [12], have been employed to remove P from Si. Implementation of these processes demands a sound thermodynamic description of the Si–P and Si–Fe–P systems.

So far, the Fe–P [13,14,15,16,17,18,19] and Fe–Si [20,21,22,23,24,25,26,27] systems have been thermodynamically modeled in many studies and were recently reoptimized by the present authors [28,29]. A complete thermodynamic modeling of the Si–P system was performed by Jung and Zhang [30] and Liang and Schmid-Fetzer [31]. However, both modeling results exhibited some discrepancies in the phase equilibria of the Si–P system and inconsistency with the ternary Si–Fe–P system. The Si–Fe–P system was thermodynamically assessed by Yan et al. [32] and Miettinen and Vassilev-Urumov [33]. The former study [32] assessed only liquid solutions using the molecular interaction volume model. In the latter assessment [33], phase equilibria and thermodynamic properties of the whole Si–Fe–P system were not well determined using the substitutional model. Therefore, a complete re-optimization of the Si–Fe–P system is necessary to obtain an accurate and self-consistent thermodynamic description of this system.

The aim of this work is to conduct a critical thermodynamic re-optimization of the Si–P and Si–Fe–P systems to develop a more accurate and reliable thermodynamic database, which is of great importance for materials engineering embracing such alloy systems. The recently optimized Fe–P [28] and Fe–Si [29] systems by the present authors were adopted in this study. The phase equilibria and thermodynamic properties of the Si–P and Si–Fe–P systems were determined based on reliable experimental data. In particular, the Gibbs energies of the liquid, solid Si, SiP(s), SiP₂(s), FeSi₄P₄(s), (Fe)₃(P,Si)₁, (Fe)₂(P,Si)₁, and (Fe)₁(P,Si)₁ phases were carefully optimized to resolve the discrepancies left in previous thermodynamic assessments [30,31,33]. The developed thermodynamic database of the Si–Fe–P system was used to predict experimentally unexplored phase diagrams and thermodynamic properties and can be applied to process optimization of Si refining and alloy design. All the calculations were performed using FactSage 8.2 software [34].

2. Thermodynamic Models

2.1. Gas Phase

The gas phase of the Si–Fe–P system is a mixture of Si(g), Si₂(g), Si₃(g), Fe(g), P(g), P₂(g), and P₄(g) species. The Gibbs energy per mole of gas ($G_T^{\mathrm{g}\mathrm{a}\mathrm{s}}$) was calculated using Equation (1):

$$G_T^{\mathrm{gas}}={\displaystyle \sum x_{\mathrm{i}}\left(G_{\mathrm{i}}^{°}+RT\mathrm{In}{x}_{\mathrm{i}}\right)}+RT\mathrm{In}\left(f/P^{\theta }\right)$$

(1)

where $x_{\mathrm{i}}$ is the mole fraction of species i, $G_{\mathrm{i}}^{°}$ is the molar Gibbs energy (J/mol) of species i that can be taken from the FactPS database stored in FactSage 8.2 software [34], R is the gas constant (=8.314 J/(mol·K)), T is the temperature in Kelvin (K), $P^{\theta }$ is the atmospheric pressure (=1 atm), and f is the gas fugacity and is identical to the gas pressure (in atm) at normal pressure.

2.2. Pure Elements and Stoichiometric Compounds

The Gibbs energies of all pure liquid and solid Si, Fe, and P elements were taken from the Scientific Group Thermodata Europe (SGTE) [35] data compilation. The Gibbs energies ($G_T^{°}$) of all intermediate stoichiometric compounds, including Fe₃P, Fe₂P, FeP, and FeP₂ of the Fe–P system, Fe₂Si, Fe₅Si₃, FeSi, FeSi₂, and Fe₃Si₇ of the Fe–Si system, SiP, SiP₂ of the Si–P system, and FeSi₄P₄ of the Si–Fe–P system were determined from their corresponding heat capacity $C_{\mathrm{P}}$ (J/(mol·K)), standard enthalpy of formation ${∆H}_{298.15\mathrm{K}}^{°}$ (J/mol) and standard entropy $S_{298.15\mathrm{K}}^{°}$ (J/(mol·K)) as expressed by Equation (2):

$$G_T^{°}=\left(\mathsf{\Delta }{H}_{298.15\mathrm{K}}^{°}+{\displaystyle {\int }_{298.15\mathrm{K}}^T{C}_P}dT\right)-T\left(S_{298.15\mathrm{K}}^{°}+{\displaystyle {\int }_{298.15\mathrm{K}}^T\frac{C_P}{T}}dT\right)$$

(2)

In the cases of pure elements and stoichiometric compounds exhibiting magnetic behavior, an additional Gibbs energy of magnetic contribution $G^{\mathrm{m}\mathrm{g}}$ (J/mol) will be applied. In this study, $G^{\mathrm{m}\mathrm{g}}$ was applied to Fe (BCC_A2, FCC_A1), Fe₃P, and Fe₅Si₃ and determined from an empirical expression proposed by Inden [36] and modified by Hillert and Jarl [37]:

$$G^{\mathrm{mg}}=RT\mathrm{In}\left(\beta +1\right)g(\tau )$$

(3)

where τ is expressed by $T/T^{\mathrm{*}}$, $T^{\mathrm{*}}$ is the Curie temperature $T_C$ (K) for ferromagnetic ordering or the Néel temperature $T_N$ (K) for anti-ferromagnetic ordering, $g\left(\tau \right)$ is a polynomial function that can be found elsewhere [37], and β is the mean magnetic moment per mole of atoms expressed in Bohr magnetons (μ_B/mol).

2.3. Solid Solutions

It has been well known that Fe and P are soluble in solid Si to form a solution in diamond_A4 cubic structure. Further, Si and P can dissolve into γ-Fe (FCC_A1) to generate a substitutional solution. In the binary Si-Fe system, the BCC phase exhibits a long-range ordering from disordered structure (BCC_A2) to ordered structure (BCC_B2) transition [29], so this order/disorder transition was thus taken into account by the present modeling for the ternary Si–Fe–P system. Furthermore, Si was considered to substitute the P atoms of Fe₃P, Fe₂P, and FeP to form (Fe)₃(P, Si)₁, (Fe)₂(P, Si)₁, and (Fe)₁(P, Si)₁ solid solutions represented by Me₃P, Me₂P, and MeP, correspondingly. The Gibbs energies of all involved solid solutions in the sub-systems of Si–Fe–P were described using the Compound Energy Formalism (CEF) with consideration of their crystal structures [38].

2.3.1. FCC_A1 and Solid Si Solutions

FCC_A1 and solid Si (diamond_A4) solutions were described with the formula (Fe, Si, P)₁(Va)₁, and their Gibbs energies per formula unit were calculated as follows:

$$\begin{array}{c}{G}_{\mathrm{S}}^{\mathrm{disorder}}={\displaystyle \sum _{\mathrm{i}=\mathrm{Fe},\mathrm{Si},\mathrm{P}}{x}_{\mathrm{i}}{G}_{\mathrm{i}}^o}+R{\rm T}{\displaystyle \sum _{\mathrm{i}=\mathrm{Fe},\mathrm{Si},\mathrm{P}}{x}_{\mathrm{i}}\mathrm{In}{x}_{\mathrm{i}}}\hfill \\ +{\displaystyle \sum _{m=0,1,2...}}{x}_{\mathrm{Fe}}{x}_{\mathrm{P}}{L}_{\mathrm{Fe},\mathrm{P}}^m+{\displaystyle \sum _{k=0,1,2...}{x}_{\mathrm{Si}}{x}_{\mathrm{P}}{L}_{\mathrm{Si},\mathrm{P}}^k}+{\displaystyle \sum _{p=0,1,2...}{x}_{\mathrm{Fe}}{x}_{\mathrm{Si}}}{L}_{\mathrm{Fe},\mathrm{Si}}^p+{\displaystyle \sum _{q=0,1,2...}{x}_{\mathrm{Fe}}{x}_{\mathrm{Si}}{x}_{\mathrm{P}}}{L}_{\mathrm{Fe},\mathrm{Si},\mathrm{P}}^q+G^{\mathrm{mg}}\hfill \end{array}$$

(4)

where $x_{\mathrm{i}}$ is the mole fraction of component i and $G_i^{°}$ is the molar Gibbs energy (J/mol) of pure solid i (i = Fe, Si, P), $L_{\mathrm{F}\mathrm{e},\mathrm{P}}^m$, $L_{\mathrm{S}\mathrm{i},\mathrm{P}}^k$, $L_{\mathrm{F}\mathrm{e},\mathrm{S}\mathrm{i}}^p$, and $L_{\mathrm{F}\mathrm{e},\mathrm{S}\mathrm{i},\mathrm{P}}^q$ are adjustable interaction parameters of corresponding binary and ternary systems (J/mol), and $G^{\mathrm{m}\mathrm{g}}$ is the magnetic contribution to the Gibbs energy (J/mol).

2.3.2. Disordered/Ordered BCC Solid Solution

The Gibbs energy of the BCC solid solution was modeled by combining the disordered part described with the formula (Fe, Si, P)₁(Va)₃ and the ordered part with the formula (Fe, Si, P)_0.5(Fe, Si, P)_0.5(Va)₃. The Gibbs energy of the disordered part can be calculated from Equation (4), while that of the ordered part can be calculated using the following equations:

$$\mathsf{\Delta }{G}_{\mathrm{BCC}}^{\mathrm{order}}=G_{\mathrm{BCC}}^{\mathrm{order}}(y_{\mathrm{i}}^{\prime },y_{\mathrm{j}}^{″})-G_{\mathrm{BCC}}^{\mathrm{order}}{(y_{\mathrm{i}}^{\prime },y_{\mathrm{j}}^{″})}_{y_{\mathrm{k}}^{\prime }=y_{\mathrm{k}}^{″}}$$

(5)

$$\begin{array}{cc}\hfill G_{\mathrm{BCC}}^{\mathrm{order}}(y_i^{\prime },y_j^{″})=& y_{\mathrm{Fe}}^{\prime }{y}_{\mathrm{Fe}}^{″}{G}_{\mathrm{Fe}:\mathrm{Fe}}+y_{\mathrm{Si}}^{\prime }{y}_{\mathrm{Si}}^{″}{G}_{\mathrm{Si}:\mathrm{Si}}+y_{\mathrm{P}}^{\prime }{y}_{\mathrm{P}}^{″}{G}_{\mathrm{P}:\mathrm{P}}\hfill \\ \hfill & +y_{\mathrm{Fe}}^{\prime }{y}_{\mathrm{Si}}^{″}{G}_{\mathrm{Fe}:\mathrm{Si}}+y_{\mathrm{Si}}^{\prime }{y}_{\mathrm{Fe}}^{″}{G}_{\mathrm{Si}:\mathrm{Fe}}+y_{\mathrm{Fe}}^{\prime }{y}_{\mathrm{P}}^{″}{G}_{\mathrm{Fe}:\mathrm{P}}\hfill \\ \hfill & +y_{\mathrm{P}}^{\prime }{y}_{\mathrm{Fe}}^{″}{G}_{\mathrm{P}:\mathrm{Fe}}+y_{\mathrm{Si}}^{\prime }{y}_{\mathrm{P}}^{″}{G}_{\mathrm{Si}:\mathrm{P}}+y_{\mathrm{P}}^{\prime }{y}_{\mathrm{Si}}^{″}{G}_{\mathrm{P}:\mathrm{Si}}\hfill \\ \hfill & \hfill +0.5RT\left(y_{\mathrm{Fe}}^{\prime }\mathrm{In}{y}_{\mathrm{Fe}}^{\prime }+y_{\mathrm{Si}}^{\prime }\mathrm{In}{y}_{\mathrm{Si}}^{\prime }+y_{\mathrm{P}}^{\prime }\mathrm{In}{y}_{\mathrm{P}}^{\prime }\right)\\ \hfill & +0.5RT\left(y_{\mathrm{Fe}}^{″}\mathrm{In}{y}_{\mathrm{Fe}}^{″}+y_{\mathrm{Si}}^{″}\ln y_{\mathrm{Si}}^{″}+y_{\mathrm{P}}^{″}\ln y_{\mathrm{P}}^{″}\right)\hfill \\ \hfill & +{\displaystyle \sum _{\mathrm{i},\mathrm{j},\mathrm{k}}{y}_{\mathrm{i}}^{\prime }{y}_{\mathrm{j}}^{\prime }{y}_{\mathrm{k}}^{″}}{L}_{\mathrm{i},\mathrm{j}:\mathrm{k}}+{\displaystyle \sum _{\mathrm{i},\mathrm{j},\mathrm{k}}{y}_{\mathrm{k}}^{\prime }{y}_{\mathrm{i}}^{″}{y}_{\mathrm{j}}^{″}}{L}_{\mathrm{k}:\mathrm{i},\mathrm{j}}+G^{\mathrm{mg}}\hfill \end{array}$$

(6)

where i, j, and k represent Fe, Si, P. $y_{\mathrm{i}}^{\prime }$, $y_{\mathrm{j}}^{\prime }$, $y_{\mathrm{k}}^{\prime }$ and $y_{\mathrm{i}}^{″}$, $y_{\mathrm{j}}^{″}$, $y_{\mathrm{k}}^{″}$ are site fractions of component i, j, k in the first and second lattice of the formula (Fe, Si, P)_0.5(Fe, Si, P)_0.5(Va)₃, respectively. The Gibbs energy of the BCC solid solution combining both ordered and disordered contributions was determined from Equation (7):

$$G_{\mathrm{BCC}}^{\mathrm{sol}.}=G_{\mathrm{S}}^{\mathrm{disorder}}+\mathsf{\Delta }{G}_{\mathrm{BCC}}^{\mathrm{order}}$$

(7)

when the site fractions of component i in the first sublattice are equal to that in the second sublattice ($y_i^{\prime }=y_i^{″}$), then the ordering contribution ${\mathsf{\Delta }G}_{\mathrm{B}\mathrm{C}\mathrm{C}}^{\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}}$ equals nil, and the Gibbs energy of the BCC phase is the same as that of the disordered BCC_A2 ($G_{\mathrm{S}}^{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}}$) calculated by Equation (4). In the case of $y_i^{\prime }\ne y_i^{″}$, then ordering contribution ${\mathsf{\Delta }G}_{\mathrm{B}\mathrm{C}\mathrm{C}}^{\mathrm{o}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{r}}$ becomes negative, and the Gibbs energy of the BCC_B2 phase can be calculated using Equation (7).

2.3.3. Other Solid Solutions (Me₃P, Me₂P, MeP)

The solid solutions Me₃P, Me₂P, and MeP were described with the formula (Fe)_n(P, Si)₁, where $\mathrm{n}=3,2,1$ for Me₃P, Me₂P, and MeP, respectively. Their Gibbs energies per formula unit were calculated based on the CEF [38] as follows:

$$G_{{\mathrm{Me}}_{\mathrm{n}}\mathrm{P}}^{\mathrm{sol}.}=y_{\mathrm{P}}{G}_{{\mathrm{Fe}}_{\mathrm{n}}\mathrm{P}}^{°}+y_{\mathrm{Si}}{G}_{{\mathrm{Fe}}_{\mathrm{n}}\mathrm{Si}}+R\mathrm{{\rm T}}\left(y_{\mathrm{P}}\ln y_{\mathrm{P}}+y_{\mathrm{Si}}\ln y_{\mathrm{Si}}\right)+{\displaystyle \sum _{m=0,1,2...}{y}_{\mathrm{P}}{y}_{\mathrm{Si}}}{L}_{\mathrm{Fe}:\mathrm{P},\mathrm{Si}}^m+G^{\mathrm{mg}}$$

(8)

here, $G_{{\mathrm{F}\mathrm{e}}_3\mathrm{P}}^{°}$, $G_{{\mathrm{F}\mathrm{e}}_2\mathrm{P}}^{°}$, and $G_{\mathrm{F}\mathrm{e}\mathrm{P}}^{°}$ are optimized Gibbs energies (J/mol) of stoichiometric ${\mathrm{F}\mathrm{e}}_3\mathrm{P}$, ${\mathrm{F}\mathrm{e}}_2\mathrm{P},$ and $\mathrm{F}\mathrm{e}\mathrm{P}$ compounds of the Fe–P system [28]; $G_{{\mathrm{F}\mathrm{e}}_3\mathrm{S}\mathrm{i}}$, $G_{{\mathrm{F}\mathrm{e}}_2\mathrm{S}\mathrm{i}}$, and $G_{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}$ are Gibbs energies (J/mol) of Fe₃Si, Fe₂Si, and FeSi combinations respectively, which need to be optimized in this work; $y_{\mathrm{P}}$ and $y_{\mathrm{S}\mathrm{i}}$ are site fractions of P and Si in the second sublattice, respectively; $L_{\mathrm{F}\mathrm{e}:\mathrm{P},\mathrm{S}\mathrm{i}}^m$ is the adjustable interaction parameter (J/mol), and $G^{\mathrm{m}\mathrm{g}}$ is the magnetic contribution to the Gibbs energy (J/mol).

2.4. Liquid Solution

The Gibbs energies of all liquid solutions within the Si–Fe–P system were described by the Modified Quasichemical Model (MQM) [39,40] in pair approximation. The MQM, with consideration of the bond structure in a liquid solution, gives a more realistic thermodynamic description of the liquid solution than the Bragg–Williams Random Mixing Model. In MQM, the pair formation Gibbs energy can be expressed as a polynomial in the pair fraction instead of the component fraction, and the coordination number of each component can be varied with composition to reproduce the short-range ordering more easily.

In the case of the binary A–B liquid phase, A atoms and B atoms are distributed over the quasi-lattice sites, and the following atom pair exchanging reaction is considered:

$$(\mathrm{A}-\mathrm{A})+(\mathrm{B}-\mathrm{B})=2(\mathrm{A}-\mathrm{B});\mathsf{\Delta }{g}_{\mathrm{AB}}$$

(9)

where (A–A), (B–B), and (A–B) represent the first-nearest-neighbor pairs between components A and A, B and B, A and B, and $\mathsf{\Delta }{g}_{\mathrm{A}\mathrm{B}}$ is the Gibbs energy change for the formation of 2 moles (A–B) pairs from 1 mole (A–A) pairs and 1 mole (B–B) pairs. The Gibbs energy of the A–B solution was calculated using Equation (10):

$$G_{\mathrm{AB}}^{\mathrm{L}}=(n_{\mathrm{A}}{G}_{\mathrm{A}}^{°}+n_{\mathrm{B}}{G}_{\mathrm{B}}^{°})-T\mathsf{\Delta }{S}_{\mathrm{AB}}^{\mathrm{conf}.}+n_{\mathrm{AB}}(\mathsf{\Delta }{g}_{\mathrm{AB}}/2)$$

(10)

where $n_{\mathrm{A}}$ and $n_{\mathrm{B}}$ are the mole numbers of A atoms and B atoms (mol), $G_{\mathrm{A}}^{°}$ and $G_{\mathrm{B}}^{°}$ are the molar Gibbs energies of pure A and B in a liquid state (J/mol), and ${\mathsf{\Delta }S}_{\mathrm{A}\mathrm{B}}^{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}.}$ is the configurational entropy of mixing (J/(mol·K)) given by Equation (11):

$$\mathsf{\Delta }{S}_{\mathrm{AB}}^{\mathrm{conf}.}=-R(n_{\mathrm{A}}\ln X_{\mathrm{A}}+n_{\mathrm{B}}\ln X_{\mathrm{B}})-R\left[n_{\mathrm{AA}}\ln \left(\frac{X_{\mathrm{AA}}}{Y_{\mathrm{A}}^2}\right)+n_{\mathrm{BB}}\ln \left(\frac{X_{\mathrm{BB}}}{Y_{\mathrm{B}}^2}\right)+n_{\mathrm{AB}}\ln \left(\frac{X_{\mathrm{AB}}}{2Y_{\mathrm{A}}{Y}_{\mathrm{B}}}\right)\right]$$

(11)

where $n_{\mathrm{A}\mathrm{A}}$, $n_{\mathrm{B}\mathrm{B}}$, and $n_{\mathrm{A}\mathrm{B}}$ represent the mole numbers of (A–A), (B–B) and (A–B) pairs (mol), $X_{\mathrm{A}\mathrm{A}}$, $X_{\mathrm{A}\mathrm{B}}$, and $X_{\mathrm{A}\mathrm{B}}$ are pair fractions of corresponding atom pairs, and $Y_{\mathrm{A}}$ and $Y_{\mathrm{B}}$ are coordination equivalent fractions of A atoms and B atoms. The pair fractions $X_{\mathrm{A}\mathrm{A}}$, $X_{\mathrm{B}\mathrm{B}}$, $X_{\mathrm{A}\mathrm{B}}$ and coordination equivalent fractions $Y_{\mathrm{A}}$, $Y_{\mathrm{B}}$ can be calculated using Equations (12)–(16):

$$X_{\mathrm{AA}}=n_{\mathrm{AA}}/\left(n_{\mathrm{AA}}+n_{\mathrm{AB}}+n_{\mathrm{BB}}\right)$$

(12)

$$X_{\mathrm{AB}}=n_{\mathrm{AB}}/\left(n_{\mathrm{AA}}+n_{\mathrm{AB}}+n_{\mathrm{BB}}\right)$$

(13)

$$X_{\mathrm{BB}}=n_{\mathrm{BB}}/\left(n_{\mathrm{AA}}+n_{\mathrm{AB}}+n_{\mathrm{BB}}\right)$$

(14)

$$Y_{\mathrm{A}}=X_{\mathrm{AA}}+\frac{1}{2}{X}_{\mathrm{AB}}$$

(15)

$$Y_{\mathrm{B}}=X_{\mathrm{BB}}+\frac{1}{2}{X}_{\mathrm{AB}}$$

(16)

${\mathsf{\Delta }g}_{\mathrm{A}\mathrm{B}}$ in Equations (9) and (10) is the model parameter for reproducing the Gibbs energy of the A–B liquid solution (J/mol) and can be expanded as a polynomial in terms of the atomic pair fractions $X_{\mathrm{A}\mathrm{A}}$ and $X_{\mathrm{B}\mathrm{B}}$.

$$\mathsf{\Delta }{g}_{\mathrm{AB}}=\mathsf{\Delta }{g}_{\mathrm{AB}}^{°}+{\displaystyle \sum _{\mathrm{i}\ge 1}{g}_{\mathrm{AB}}^{\mathrm{i}0}}{X}_{\mathrm{AA}}^{\mathrm{i}}+{\displaystyle \sum _{\mathrm{j}\ge 1}{g}_{\mathrm{AB}}^{0\mathrm{j}}}{X}_{\mathrm{BB}}^{\mathrm{j}}$$

(17)

where ${\mathsf{\Delta }g}_{\mathrm{A}\mathrm{B}}^{°}$, $g_{\mathrm{A}\mathrm{B}}^{\mathrm{i}0}$ and $g_{\mathrm{A}\mathrm{B}}^{0\mathrm{j}}$ are the adjustable model parameters (J/mol) that can be functions of the temperature. In the MQM, the coordination numbers of A and B, and $Z_{\mathrm{A}}$ and $Z_{\mathrm{B}}$, can be varied with the composition to reproduce the short-range ordering.

$$\frac{1}{Z_{\mathrm{A}}}=\frac{1}{Z_{\mathrm{AA}}^{\mathrm{A}}}\left(\frac{2n_{\mathrm{AA}}}{2n_{\mathrm{AA}}+n_{\mathrm{AB}}}\right)+\frac{1}{Z_{\mathrm{AB}}^{\mathrm{A}}}\left(\frac{n_{\mathrm{AB}}}{2n_{\mathrm{AA}}+n_{\mathrm{AB}}}\right)$$

(18)

$$\frac{1}{Z_{\mathrm{B}}}=\frac{1}{Z_{\mathrm{BB}}^{\mathrm{B}}}\left(\frac{2n_{\mathrm{BB}}}{2n_{\mathrm{BB}}+n_{\mathrm{AB}}}\right)+\frac{1}{Z_{\mathrm{BA}}^{\mathrm{B}}}\left(\frac{n_{\mathrm{AB}}}{2n_{\mathrm{BB}}+n_{\mathrm{AB}}}\right)$$

(19)

here $Z_{\mathrm{A}\mathrm{A}}^{\mathrm{A}}$ is the value $Z_{\mathrm{A}}$ when all nearest neighbors of an A atom are A atoms, and $Z_{\mathrm{A}\mathrm{B}}^{\mathrm{A}}$ is the value of $Z_{\mathrm{A}}$ when all nearest neighbors of the A atom are B atoms. $Z_{\mathrm{B}\mathrm{B}}^{\mathrm{B}}$ and $Z_{\mathrm{B}\mathrm{A}}^{\mathrm{B}}$ are defined in an analogous manner. In the present study, $Z_{\mathrm{F}\mathrm{e}\mathrm{F}\mathrm{e}}^{\mathrm{F}\mathrm{e}}=Z_{\mathrm{S}\mathrm{i}\mathrm{S}\mathrm{i}}^{\mathrm{S}\mathrm{i}}=Z_{\mathrm{P}\mathrm{P}}^{\mathrm{P}}=6$ [28,29], $Z_{\mathrm{P}\mathrm{F}\mathrm{e}}^{\mathrm{P}}=Z_{\mathrm{P}\mathrm{S}\mathrm{i}}^{\mathrm{P}}=Z_{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}^{\mathrm{F}\mathrm{e}}=Z_{\mathrm{S}\mathrm{i}\mathrm{F}\mathrm{e}}^{\mathrm{S}\mathrm{i}}=Z_{\mathrm{S}\mathrm{i}\mathrm{P}}^{\mathrm{S}\mathrm{i}}=6$ [28,29,30], and $Z_{\mathrm{F}\mathrm{e}\mathrm{P}}^{\mathrm{F}\mathrm{e}}=3$ [28], as listed in

Table 1. Optimized model parameters for the Si–Fe–P system. Heat capacity $C_P$ (J/(mol·K)), standard enthalpy of formation ${∆H}_{298.15K}^{°}$ (J/mol), standard entropy $S_{298.15K}^{°}$ (J/(mol·K), adjustable interaction parameter L (J/mol), Gibbs energy G, g, $∆\mathrm{g}$ (J/mol), Curie temperature $T_C$ (K), magnetic moment $\beta$ (μ_B/mol).

Phase	Model Parameters
Liquid (Si, Fe, P)	$Z_{\mathrm{FeFe}}^{\mathrm{Fe}}=Z_{\mathrm{SiSi}}^{\mathrm{Si}}=Z_{\mathrm{PP}}^{\mathrm{P}}=6$ [28,29]
	$Z_{\mathrm{PFe}}^{\mathrm{P}}=Z_{\mathrm{PSi}}^{\mathrm{P}}=Z_{\mathrm{FeSi}}^{\mathrm{Fe}}=Z_{\mathrm{SiFe}}^{\mathrm{Si}}=Z_{\mathrm{SiP}}^{\mathrm{Si}}=6$ [28,29,30] $Z_{\mathrm{FeP}}^{\mathrm{Fe}}=3$ [28] $\Delta g_{\mathrm{FeP}}=-56902+6.569T+(5481+3.033T)X_{\mathrm{FeFe}}-(11966-2.51T)X_{\mathrm{FeFe}}^2-9623X_{\mathrm{PP}}$ [28] $\Delta g_{\mathrm{SiP}}=-8812+2.092T-2343X_{\mathrm{SiSi}}+6276X_{\mathrm{PP}}$ [*] $\Delta g_{\mathrm{FeSi}}=-33710+2.26T-(12552-5.02T)X_{\mathrm{FeFe}}-(8368-4.81T)X_{\mathrm{FeFe}}^2-(3054-6.49T)X_{\mathrm{SiSi}}$ [29] $g_{\mathrm{FeSi}(\mathrm{P})}^{001}=-3849+5.86T$ [*], $g_{\mathrm{FeSi}(\mathrm{P})}^{011}=14770$ [*], $g_{\mathrm{FeSi}(\mathrm{P})}^{101}=-34095+15.06T$ [*]
	“Toop-like” interpolation with Fe as an asymmetric component [*]
FCC_A1 (Si, Fe, P)1(Va)1	$G_{\mathrm{Fe}:\mathrm{Va}}^{\mathrm{FCC}}=G_{\mathrm{Fe}(\mathrm{FCC})}^{°}$, $G_{\mathrm{Si}:\mathrm{Va}}^{\mathrm{FCC}}=G_{\mathrm{Si}(\mathrm{FCC})}^{°}$, $G_{\mathrm{P}:\mathrm{Va}}^{\mathrm{FCC}}=G_{\mathrm{P}(\mathrm{FCC})}^{°}$ [*] $L_{\mathrm{Fe},\mathrm{P}:\mathrm{Va}}^{\mathrm{FCC}}=-139787+6.49T$ [28] $L_{\mathrm{Si},\mathrm{P}:\mathrm{Va}}^{\mathrm{FCC}}=0$ [*] $L_{\mathrm{Fe},\mathrm{Si}:\mathrm{Va}}^{\mathrm{FCC}}=-115254-2.19T-(84777-44.33T)\left(x_{\mathrm{Fe}}-x_{\mathrm{Si}}\right)+20007{\left(x_{\mathrm{Fe}}-x_{\mathrm{Si}}\right)}^2$ [29] $T_{C\mathrm{Fe}:\mathrm{Va}}=-201$, ${\beta }_{\mathrm{Fe}:\mathrm{Va}}=-2.1$ [41]
BCC_A2 (Si, Fe, P)1(Va)3	$G_{\mathrm{Fe}:\mathrm{Va}}^{\mathrm{BCC}\text{\_}\mathrm{A}2}=G_{\mathrm{Fe}(\mathrm{BCC})}^{°}$, $G_{\mathrm{Si}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{A}2}=G_{\mathrm{Si}(\mathrm{BCC})}^{°}$, $G_{\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{A}2}=G_{\mathrm{P}(\mathrm{BCC})}^{°}$ [*]
	$L_{\mathrm{Fe},\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{A}2}=-203476+15.48T+33472\left(y_{\mathrm{Fe}}-y_{\mathrm{P}}\right)$ [28]
	$L_{\mathrm{Si},\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{A}2}=-52300$ [*]
	$L_{\mathrm{Fe},\mathrm{Si}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{A}2}=-154014+32.29T-(63511-13.25T)(x_{\mathrm{Fe}}-x_{\mathrm{Si}})+35728{(x_{\mathrm{Fe}}-x_{\mathrm{Si}})}^2$ [29]$T_{C\mathrm{Fe},\mathrm{P}:\mathrm{Va}}=-285$ [28], $T_{C\mathrm{Fe},\mathrm{Si}:\mathrm{Va}}=504(y_{\mathrm{Fe}}-y_{\mathrm{Si}})$ [29] $T_{C\mathrm{Fe}:\mathrm{Va}}=1043$, ${\beta }_{\mathrm{Fe}:\mathrm{Va}}=2.22$ [41]
BCC_B2 (Si, Fe, P)0.5(Si, Fe, P)0.5(Va)3	$G_{\mathrm{Fe}:\mathrm{Si}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{Si}:\mathrm{Fe}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=-20930$ [29] $G_{\mathrm{Fe}:\mathrm{Fe}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{Si}:\mathrm{Si}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=0$ [29] $G_{\mathrm{P}:\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{Fe}:\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{P}:\mathrm{Fe}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{Si}:\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=G_{\mathrm{P}:\mathrm{Si}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=0$ [*] $L_{\mathrm{Fe},\mathrm{Si}:\mathrm{Si}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Fe},\mathrm{Si}:\mathrm{Fe}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Fe}:\mathrm{Fe},\mathrm{Si}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Si}:\mathrm{Fe},\mathrm{Si}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=0$ [29]$L_{\mathrm{Fe},\mathrm{P}:\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Fe},\mathrm{P}:\mathrm{Fe}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Fe}:\mathrm{Fe},\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{P}:\mathrm{Fe},\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=0$ [28]
	$L_{\mathrm{Si},\mathrm{P}:\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Si},\mathrm{P}:\mathrm{Fe}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{Si}:\mathrm{Si},\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{P}:\mathrm{Si},\mathrm{P}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=0$ [*] $L_{\mathrm{Fe},\mathrm{Si}\text{:}\mathrm{P}:\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=L_{\mathrm{P}:\mathrm{Fe},\mathrm{Si}\text{:}\mathrm{Va}}^{{\mathrm{BCC}}_{-}\mathrm{B}2}=-3870$ [*]
Diamond_A4 Si (Si, Fe, P)1(Va)1	$G_{\mathrm{Si}:\mathrm{Va}}=G_{\mathrm{Si}(\mathrm{diamond}\text{\_}\mathrm{A}4)}^{°}$ [29], $G_{\mathrm{Fe}:\mathrm{Va}}=G_{\mathrm{Fe}(\mathrm{FCC})}^{°}+1000$ [29] $G_{\mathrm{P}:\mathrm{Va}}=G_{\mathrm{P}(\mathrm{white})}^{°}+67100-2T$ [*]$L_{\mathrm{Si},\mathrm{P}:\mathrm{Va}}^{\mathrm{Diamond}\text{\_}\mathrm{A}4}=-50208$ [*] $L_{\mathrm{Si},\mathrm{Fe}:\mathrm{Va}}^{\mathrm{Diamond}\text{\_}\mathrm{A}4}=113001-0.5T$ [29]
Me3P (Fe)3(P, Si)1	$G_{\mathrm{Fe}:\mathrm{P}}^{{\mathrm{Me}}_3\mathrm{P}}=G_{{\mathrm{Fe}}_3\mathrm{P}}^{°}$ [28] $G_{\mathrm{Fe}:\mathrm{Si}}^{{\mathrm{Me}}_3\mathrm{P}}=3G_{\mathrm{Fe}(\mathrm{BCC})}^{°}+G_{\mathrm{Si}(\mathrm{Diamond}\text{\_}\mathrm{A}4)}^{°}-94500$ [*] $L_{\mathrm{Fe}:\mathrm{P},\mathrm{Si}}^{{\mathrm{Me}}_3\mathrm{P}}=0$ [*]
Me2P (Fe)2(P, Si)1	$G_{\mathrm{Fe}:\mathrm{P}}^{{\mathrm{Me}}_2\mathrm{P}}=G_{{\mathrm{Fe}}_2\mathrm{P}}^{°}$ [28] $G_{\mathrm{Fe}:\mathrm{Si}}^{{\mathrm{Me}}_2\mathrm{P}}=2G_{\mathrm{Fe}(\mathrm{BCC})}^{°}+G_{\mathrm{Si}(\mathrm{Diamond}\_\mathrm{A}4)}^{°}-116480+24T$ [*] $L_{\mathrm{Fe}:\mathrm{P},\mathrm{Si}}^{{\mathrm{Me}}_2\mathrm{P}}=0$ [*]
MeP (Fe)1(P, Si)1	$G_{\mathrm{Fe}:\mathrm{P}}^{\mathrm{MeP}}=G_{\mathrm{FeP}}^{°}$ [28] $G_{\mathrm{Fe}:\mathrm{Si}}^{\mathrm{MeP}}=G_{\mathrm{Fe}(\mathrm{BCC})}^{°}+G_{\mathrm{Si}(\mathrm{Diamond}\text{\_}\mathrm{A}4)}^{°}-87350+17T$ [*] $L_{\mathrm{Fe}:\mathrm{P},\mathrm{Si}}^{\mathrm{MeP}}=14644$ [*]
FeP2 (Fe)1(P)2	$\Delta H_{298.15K}^{°}=-191100$, $S_{298.15K}^{°}=51.05$ [28] $C_{\mathrm{P}}=77.52563+0.009348T-443846T^{-2}-1.1\times {10}^{-6}{T}^2$ [28]
SiP (Si)1(P)1	$\Delta H_{298.15K}^{°}=-64000$ [*], $S_{298.15K}^{°}=34.78$ [*] $C_{\mathrm{P}}=38.343+0.010878T-565000T^{-2}$ [30]
SiP2 (Si)1(P)2	$\Delta H_{298.15K}^{°}=-79950$ [*], $S_{298.15K}^{°}=64$ [*] $C_{\mathrm{P}}=67+0.0171T$ [30]
Fe2Si (Fe)2(Si)	$\Delta H_{298.15K}^{°}=-53889.7$ [29], $S_{298.15K}^{°}=106.39$ [29] $C_{\mathrm{P}}=2C_{\mathrm{P}}({\mathrm{Fe}}_{\mathrm{BCC}})+C_{\mathrm{P}}({\mathrm{Si}}_{\mathrm{Diamond}\_\mathrm{A}4})$ [29]
Fe5Si3 (Fe)5(Si)3	$\Delta H_{298.15K}^{°}=-234740$ [29],$S_{298.15K}^{°}=209.1$ [29] $C_{\mathrm{P}}=180.3069+0.085912T-1060722T^{-2}+2.665\times {10}^{-7}{T}^2$ [29] $\beta =2.32$ [29], $T_{\mathrm{C}}=360$ [29], $P=0.28$ [29]
FeSi (Fe)1(Si)1	$\Delta H_{298.15K}^{°}=-76410$ [29], $S_{298.15K}^{°}=46.024$ [29] $C_{\mathrm{P}}=48.5666+0.01472T-428220T^{-2}-1.7511\times {10}^{-6}{T}^2$ [29]
FeSi2 (Fe)1(Si)2	$\Delta H_{298.15K}^{°}=-96940.44$ [29], $S_{298.15K}^{°}=55.48$ [29] $C_{\mathrm{P}}=79.02985-0.0181469T-999009T^{-2}-1.782\times {10}^{-6}{T}^2$ [29]
Fe3Si7 (Fe)3(Si)7	$\Delta H_{298.15K}^{°}=-247842.42$ [29], $S_{298.15K}^{°}=207.3$ [29] $C_{\mathrm{P}}=214.2176+0.10993T-2345707T^{-2}-2.3033\times {10}^{-6}{T}^2$ [29]
FeSi4P4 (Fe)1(Si)4(P)4	$\Delta H_{298.15K}^{°}=-334600$ [*], $S_{298.15K}^{°}=175$ [*] $C_{\mathrm{P}}=C_{\mathrm{P}}({\mathrm{Fe}}_{\mathrm{BCC}})+4C_{\mathrm{P}}({\mathrm{Si}}_{\mathrm{Diamond}\text{\_}\mathrm{A}4})+4C_{\mathrm{P}}({\mathrm{P}}_{\mathrm{White}})$ [*]

* optimized in the present study.

.

The Gibbs energy of the ternary Si–Fe–P liquid solution can be predicted from the interpolation of Gibbs energies of its sub-binary systems based on their nature. In the present study, an asymmetric “Toop-like” geometric interpolation [40] with Fe as the “asymmetric component” was used for the Si–Fe–P system since the Fe–P and Fe–Si liquid solutions show much more negative deviations from the Si–P liquid solution. Based on this interpolation, the entropy of mixing and Gibbs energy of the ternary Si–Fe–P liquid solution were calculated using Equations (20) and (21):

$$\mathsf{\Delta }{S}_{\mathrm{Si},\mathrm{Fe},\mathrm{P}}^{\mathrm{conf}.}=-R{\displaystyle \sum _{\mathrm{i}=\mathrm{Si},\mathrm{Fe},\mathrm{P}}{n}_{\mathrm{i}}\ln X_{\mathrm{i}}}-R\left[{\displaystyle \sum _{\mathrm{j}=\mathrm{Si},\mathrm{Fe},\mathrm{P}}{n}_{\mathrm{jj}}}\left(\frac{X_{\mathrm{jj}}}{Y_{\mathrm{j}}^2}\right)+{\displaystyle \sum _{\mathrm{k},\mathrm{m}=\mathrm{Si},\mathrm{Fe},\mathrm{P}}^{\mathrm{k}\ne \mathrm{m}}{n}_{\mathrm{km}}\ln \left(\frac{X_{\mathrm{km}}}{2Y_{\mathrm{k}}{Y}_{\mathrm{m}}}\right)}\right]$$

(20)

$$G_{\mathrm{Si},\mathrm{Fe},\mathrm{P}}^{\mathrm{L}}={\displaystyle \sum _{\mathrm{i}=\mathrm{Si},\mathrm{Fe},\mathrm{P}}{n}_{\mathrm{i}}{G}_{\mathrm{i}}^{°}}-T\mathsf{\Delta }{S}_{\mathrm{Si},\mathrm{Fe},\mathrm{P}}^{\mathrm{conf}.}+{\displaystyle \sum _{\mathrm{j},\mathrm{k}=\mathrm{Si},\mathrm{Fe},\mathrm{P}}^{\mathrm{j}\ne \mathrm{k}}\left(n_{\mathrm{jk}}/2\right)\mathsf{\Delta }{g}_{\mathrm{jk}}}$$

(21)

here, ${\mathsf{\Delta }g}_{\mathrm{j}\mathrm{k}}$ ($\mathrm{j},\mathrm{k}=\mathrm{S}\mathrm{i},\mathrm{F}\mathrm{e},\mathrm{P}$) is the pair formation Gibbs energy depending on the thermodynamic symmetry of the ternary system. ${\mathsf{\Delta }g}_{\mathrm{F}\mathrm{e}\mathrm{P}}$ and ${\mathsf{\Delta }g}_{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}$ for the asymmetric Fe–P and Fe–Si systems are respectively calculated using Equations (22) and (23):

$$\mathsf{\Delta }{g}_{\mathrm{FeSi}}=\mathsf{\Delta }{g}_{\mathrm{FeSi}}^{°}+{\displaystyle \sum _{(\mathrm{i}+\mathrm{j})\ge 1}{g}_{\mathrm{FeSi}}^{\mathrm{ij}}}{x}_{\mathrm{FeFe}}^{\mathrm{i}}{\left(x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}\right)}^{\mathrm{j}}+{\displaystyle \sum _{\mathrm{i}\ge 0,\mathrm{j}\ge 0,\mathrm{k}\ge 1}{g}_{\mathrm{FeSi}(\mathrm{P})}^{\mathrm{ijk}}}{x}_{\mathrm{FeFe}}^{\mathrm{i}}{\left(x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}\right)}^{\mathrm{j}}{\left(\frac{Y_{\mathrm{P}}}{Y_{\mathrm{Si}}+Y_{\mathrm{P}}}\right)}^{\mathrm{k}}$$

(22)

$$\mathsf{\Delta }{g}_{\mathrm{FeP}}=\mathsf{\Delta }{g}_{\mathrm{FeP}}^{°}+{\displaystyle \sum _{(\mathrm{i}+\mathrm{j})\ge 1}{g}_{\mathrm{FeP}}^{\mathrm{ij}}}{x}_{\mathrm{FeFe}}^{\mathrm{i}}{\left(x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}\right)}^{\mathrm{j}}+{\displaystyle \sum _{\mathrm{i}\ge 0,\mathrm{j}\ge 0,\mathrm{k}\ge 1}{g}_{\mathrm{FeP}(\mathrm{Si})}^{\mathrm{ijk}}}{x}_{\mathrm{FeFe}}^{\mathrm{i}}{\left(x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}\right)}^{\mathrm{j}}{\left(\frac{Y_{\mathrm{Si}}}{Y_{\mathrm{Si}}+Y_{\mathrm{P}}}\right)}^{\mathrm{k}}$$

(23)

while ${\mathsf{\Delta }g}_{\mathrm{S}\mathrm{i}\mathrm{P}}$ for the symmetric Si–P system was calculated using Equation (24):

$$\mathsf{\Delta }{g}_{\mathrm{SiP}}=\mathsf{\Delta }{g}_{\mathrm{SiP}}^{°}+{\displaystyle \sum _{(\mathrm{i}+\mathrm{j})\ge 1}{g}_{\mathrm{SiP}}^{\mathrm{ij}}}{\left(\frac{x_{\mathrm{SiSi}}}{x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}}\right)}^{\mathrm{i}}{\left(\frac{x_{\mathrm{PP}}}{x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}}\right)}^{\mathrm{j}}+{\displaystyle \sum _{\mathrm{i}\ge 0,\mathrm{j}\ge 0,\mathrm{k}\ge 1}{g}_{\mathrm{SiP}(\mathrm{Fe})}^{\mathrm{ijk}}}{\left(\frac{x_{\mathrm{SiSi}}}{x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}}\right)}^{\mathrm{i}}{\left(\frac{x_{\mathrm{PP}}}{x_{\mathrm{SiSi}}+x_{\mathrm{SiP}}+x_{\mathrm{PP}}}\right)}^{\mathrm{j}}{Y}_{\mathrm{Fe}}^{\mathrm{k}}$$

(24)

where ${\mathrm{g}}_{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}^{\mathrm{i}\mathrm{j}}$, ${\mathrm{g}}_{\mathrm{F}\mathrm{e}\mathrm{P}}^{\mathrm{i}\mathrm{j}}$, and ${\mathrm{g}}_{\mathrm{S}\mathrm{i}\mathrm{P}}^{\mathrm{i}\mathrm{j}}$ are binary liquid model parameters (J/mol) and ${\mathrm{g}}_{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{P}\right)}^{\mathrm{i}\mathrm{j}\mathrm{k}}$, ${\mathrm{g}}_{\mathrm{F}\mathrm{e}\mathrm{P}\left(\mathrm{S}\mathrm{i}\right)}^{\mathrm{i}\mathrm{j}\mathrm{k}}$, and ${\mathrm{g}}_{\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{F}\mathrm{e}\right)}^{\mathrm{i}\mathrm{j}\mathrm{k}}$ are ternary liquid model parameters (J/mol).

3. Critical Evaluation and Thermodynamic Optimization

Thermodynamic optimization of the Si–P and Si–Fe–P systems was performed using the CALPHAD approach based on the critical evaluation of all available phase equilibria and thermodynamic property data. The liquid and solid solutions of all sub-systems were modeled using the MQM [39,40] and CEF [38], respectively. The optimized model parameters of the Si-Fe-P system are summarized in Table 1 and the crystal structure information of all solid phases of this system in

Table 2. Summary of the crystal structure information of all solid phases of the Si–Fe–P system.

Phase	Structure	Prototype	Space Group	Pearson Symbol
FCC_A1	Cubic	Cu	$\stackrel{-}{F}$m$\stackrel{-}{3}$m	cF4
BCC_A2	Cubic	W	Im$\stackrel{-}{3}$m	cI2
BCC_B2	Cubic	CsCl	Pm$\stackrel{-}{3}$m	cP8
Diamond_A4 Si	Cubic	C(dia.)	Fd$\stackrel{-}{3}$m	cF8
Me3P	Tetragonal	Ni3P	I$\stackrel{-}{4}$	tI32
Me2P	Hexagonal	Fe2P	P$\stackrel{-}{6}$2m	hP9
MeP	Orthorhombic	MnP	Pnma	oP8
FeP2	Orthorhombic	FeS2	Pnnm	oP6
SiP	Orthorhombic	SiP	Cmc21	oS24
SiP2	Orthorhombic	GeAs2	Pbam	oP24
Fe2Si	Cubic	CsCl	Pm$\stackrel{-}{3}$m	cP2
Fe5Si3	Hexagonal	Mn5Si3	P63/mcm	hP16
FeSi	Cubic	FeSi	P213	cP8
FeSi2	Orthorhombic	FeSi2	Cmca	oC48
Fe3Si7	Tetragonal	Fe3Si7	P4/mmm	tP3
FeSi4P4	Triclinic	FeSi4P4	P1	---
White P	Cubic	P4	I$\stackrel{-}{4}$3m	C*8
Red P	---	P	---	C*66

.

3.1. The Si–P System

The Si–P system has been well-reviewed by Mostafa [42], Jung and Zhang [30], and Liang and Schmid-Fetzer [31]. According to the literature, liquid solution, solid Si (diamond_A4), red P, SiP, and SiP₂ are stable condensed phases in the Si–P system and are also adopted by the present study.

3.1.1. Phase Diagram

The calculated Si–P phase diagrams from previous assessments [30,31] and the present study are plotted in Figure 1, along with the experimental data [43,44,45,46,47,48,49,50,51,52,53,54,55]. The phase diagrams with suppression of the gas phase are plotted in Figure 1a, and those with gas phase in a total pressure of 0.01, 0.1, 0.5, and 1 atm are in Figure 1b. In the thermodynamic assessment by Jung and Zhang [30], a constant MQM parameter was used to describe the Si–P liquid solution; however, their calculation of liquidus boundaries shows some deviations from the experimental results of the particularly high-P region [43,55]. The assessed eutectic composition ($\mathrm{w}\mathrm{t}.\%\mathrm{P}=42.9$) for the reaction $\mathrm{l}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}+\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{s}\right)$ in their study was found to be much overestimated, which is why the liquidus data nearing the SiP compound could not be reproduced. As a result, SiP was calculated to melt congruently at 1410 K to reproduce the data (1412 ± 2 K) suggested by Safarian and Tangstad [54]. However, the reliability of this data is doubtful for given reasons [31]. In the subsequent assessment by Liang and Schmid-Fetzer [31] using the substitutional solution model, the eutectic reaction $\mathrm{l}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}+\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{s}\right)$ was modified to occur at $\mathrm{w}\mathrm{t}.\%\mathrm{P}=36.6$ and 1404 K, and the melting points of SiP and SiP₂ were increased respectively to 1434 K and 1443 K to bridge the gap with the experimental data of Ugai et al. [55]. These data were measured from samples in the actual composition range of the compounds and considered more reliable. Compared to the assessment of Jung and Zhang [30], the liquidus boundaries calculated by Liang and Schmid-Fetzer [31] show certain improvements, as shown in Figure 1a. The eutectic composition and temperature proposed by Liang and Schmid-Fetzer [31] were inherited by the present study. Meanwhile, the liquidus temperatures were increased to fit better with the experimental results [43,55]. In particular, the melting temperatures of SiP and SiP₂ were further optimized to 1440 K and 1451 K, respectively, to match the data of Ugai et al. [55]. However, it is worth noting that the discrepancies between the solidus and liquidus temperatures in the region of SiP to SiP₂ could not be resolved [30,31]. The eutectic reaction $\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{s}\right)+{\mathrm{S}\mathrm{i}\mathrm{P}}_2\left(s\right)$ was reported by Ugai et al. [55] to occur at 1398 K, which is 53 K lower than the melting point of SiP₂. If such a big difference was constrainedly reconciled, then P-rich side parameters with high-temperature dependence had to be introduced to the liquid phase. This would, unfortunately, result in a “two-liquid phase” immiscibility gap at very high temperatures and also conflicts with phase equilibria of the Si-rich region. In the present study, the eutectic temperature was determined to be 1436 K, exhibiting a difference of 15 K from the SiP₂ melting point.

Figure 2 shows calculated Si–P phase diagrams of the Si-rich region compared to experimental data [43,44,45,46,47,48,49,50,51,52,53,54]. Since these phase diagram data have been reviewed previously [30,31,42], it is not necessary to review them again in the present study. The solubility data for P in solid Si are scattered significantly below 1473 K, as shown in Figure 2. It is hard to judge the accuracy of these data just from their experimental techniques. In the modeling of Liang and Schmid-Fetzer [31], more weight was given to the solidus data of Safarian and Tangstad [54] at eutectic temperature and solvus data of Nobili [51] and Borisenko and Yudin [52] in the determination of solid Si boundary. As a result, a maximum solubility of wt.%P = 1.2 in Si at the eutectic temperature (1404 K) was calculated. According to Jung and Zhang [30], a much higher P solubility limit of wt.%P = 4.2 at 1400 K was calculated to reproduce the higher-temperature solidus data of Trumbore [44], Kooi [46], and Safarian and Tangstad [54], which show good consistency within experimental errors. From a thermodynamic modeling point of view, it is not possible to reproduce the high-temperature and low-temperature solidus data of Safarian and Tangstad [54] simultaneously. The inconsistency between these data [54] was not clarified in both previous assessments. After careful examination of their experiments [54], it was found that mass loss happened continuously, even after the dissociation of SiP, due to the evaporation of P from the samples. Therefore, delayed chemical analysis after TG/DSC analysis would lead to an underestimation of the P content of the solid Si phase. In the present optimization, the Si phase boundaries above the eutectic temperature by Jung and Zhang [30] were adopted with slight modification, while the boundaries below the eutectic temperature were pushed to the P-richer side using a much negative interaction parameter $L_{\mathrm{S}\mathrm{i},\mathrm{P}:\mathrm{V}\mathrm{a}}^{\mathrm{D}\mathrm{i}\mathrm{a}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{d}\_\mathrm{A}4}=-50208$ J/mol (Table 1) to reproduce the experimental data of Trumbore [44], Kooi [46], and Soimi et al. [53], which are in good consistency, as shown in Figure 2.

The distribution coefficient of P in Si, represented by $L_{\mathrm{P}}=C_{\mathrm{P},\mathrm{S}}/C_{\mathrm{P},\mathrm{L}}$ ($C_{\mathrm{P},\mathrm{S}}$ and $C_{\mathrm{P},\mathrm{L}}$: concentration of P in solid and liquid Si, respectively), is a key factor in describing the behavior of P during the crystallizing and melting of Si. The equilibrium $L_{\mathrm{P}}$ was reported to be 0.038 by Hall [56] based on conductivity measurements, 0.09 by Safarian and tangstad [54] from the solidus and liquidus boundaries of Si–P alloys, and 0.123 by Borisenko and Yudin [52] from the enthalpy change of P in solid and liquid Si phases. A much higher distribution coefficient ($L_{\mathrm{P}}=0.35$) was originally obtained by Struthers [57] by means of the Czochralski crystal growth and radiochemical analysis. This value was cited as the equilibrium distribution coefficient in various literature [44,58,59,60,61], which was even mistaken as “new” sources of experimental data [5,54,62]. Huff et al. [62] investigated the distribution coefficient of P in Si under different Czochralski crystal pulling rates. As a result, an “effective” distribution coefficient $L_{\mathrm{P}}=0.32$ at the pulling rate of 1.1 mils/sec and $L_{\mathrm{P}}=0.42$ at the pulling rate of 1.8 mils/sec, were derived. According to Huff et al. [62], the effective $L_{\mathrm{P}}$ is dependent on the crystal growth rate, impurity concentration, interface orientation, temperature gradients, stirring conditions, etc. Recently, Li et al. [6] determined the effective $L_{\mathrm{P}}=0.31$ and 0.33 from directional solidification ingots grown at the rate of ${2.08\times 10}^{-6}$ m/s and ${3.08\times 10}^{-6}$ m/s, respectively. It is noticeable that $L_{\mathrm{P}}$ values resulting from the Czochralski crystal growth method [57,62] and the directional solidification method [6] are typically higher than those from conductivity and thermochemical equilibrium measurements [52,54,56]. Because such Czochralski crystal growth and directional solidification experiments [6,57,62] proceeded typically at non-equilibrium conditions depending on the Czochralski rate and temperature gradient, so very slight supercooling of dilute Si–P alloys can lead to distinct concentration of P in the primary Si crystals. Therefore, these higher $L_{\mathrm{P}}$ values, which are “so-called” effective distribution coefficients, should be overestimated compared to equilibrium ones. The equilibrium $L_{\mathrm{P}}$ can be determined from the solidus and liquidus boundaries depending on the temperature. Based on the optimized Si–P phase diagram, as presented above, $L_{\mathrm{P}}$ increases with the decline in temperature and was determined to be 0.06 at the Si melting temperature (1686.95 K) and 0.11 at the eutectic temperature (1404 K).

3.1.2. Thermodynamic Stability of Si Phosphides

The Gibbs energies of intermediate silicon phosphides SiP and SiP₂, which were calculated from their heat capacity ($C_P$), standard enthalpy of formation (${∆H}_{298.15K}^{°}$), and standard entropy ($S_{298.15K}^{°}$) as expressed by Equation (2), were reoptimized in this study to improve the thermodynamic description of the Si–P system. The $C_P$ of SiP and SiP₂ were directly taken from the assessment of Jung and Zhang [30]. Ugai et al. [63] determined $S_{298.15K}^{°}=34.78$ J/(mol·K) for SiP through its low-temperature $C_P$ data. This value was adopted by the present study without modification. The standard enthalpy of formation for SiP was optimized to $-\mathrm{64,000}$ J/mol to reproduce its melting point data from Ugai et al. [55], as shown in Figure 1a, and SiP dissociation pressure data [55,64,65] presented in Figure 3. On the other hand, since no reliable Gibbs energy data of SiP₂ are available in the literature, so ${∆H}_{298.15K}^{°}$ and $S_{298.15K}^{°}$ for SiP₂ were optimized to be $-\mathrm{79,950}$ J/mol and $64$ J/(mol·K), respectively, in the present study to reproduce its melting point data [55] and make SiP₂ stable down to the ice point. The optimized thermodynamic properties of SiP and SiP₂ are given in Table 1.

3.1.3. Thermodynamic Properties of the Si–P Liquid Solution

Thermodynamic properties of the Si–P liquid solution at any designated composition and temperature can be predicted from the developed Si–P database. As solar-grade Si demands super high purity with P content down to ${10}^{-5}$ wt.% level, so much attention has been paid to the thermodynamic properties of the dilute region because it is of primary importance to the Si refining process. Based on the present thermodynamic optimization, the Henrian activity coefficient of P in Si(l), ${\gamma }_{P\mathrm{i}\mathrm{n}\mathrm{S}\mathrm{i}\left(\mathrm{l}\right)}^{°}$, depending on the temperature, was determined:

$$\mathrm{l}\mathrm{n}{\gamma }_{P\mathrm{i}\mathrm{n}\mathrm{S}\mathrm{i}\left(\mathrm{l}\right)}^{°}=-\frac{4025}{T}+0.7548,1687\mathrm{K}<T<2173\mathrm{K}$$

(25)

where T is the temperature in Kevin (K). Based on the optimized ${\gamma }_{P\mathrm{i}\mathrm{n}\mathrm{S}\mathrm{i}\left(\mathrm{l}\right)}^{°}$, the molar Gibbs energy for the dissolution of P₂(g) into liquid Si (1 wt.% standard state) was determined as follows:

$$\begin{gathered} 0.5P_2\left(g\right)=P_{\mathrm{i}\mathrm{n}\mathrm{S}\mathrm{i}\left(\mathrm{l}\right)}(\mathrm{w}\mathrm{t}.\%); \\ {∆G}_T^{°}=-92455+17.875T (\mathrm{J}/\mathrm{mol}); 1687\mathrm{K}<T<2173\mathrm{K} \end{gathered}$$

(26)

The calculated Gibbs energies of P₂(g) dissolution in Si(l) from previous assessments and the present study are plotted in Figure 4, along with the experimental data of Miki et al. [66] measured using the transportation method and those of Zaitsev et al. [67] with the Knudsen effusion mass spectrometry. These two sets of data show a distinct difference of about 15 kJ/mol. As has been pointed out [31], the activity data of Zaitsev et al. [67] are less self-consistent. Thus, their Gibbs energy data from the same experiments were not considered in the present thermodynamic modeling either. Instead, the experimental data of Miki et al. [66] were favored by the present study, which agrees with the calculation result of Liang and Vassilev-Urumov [31]. However, the modeling result of Jung and Zhang [30] shows a large deviation from both sets of experimental data [66,67], as shown in Figure 4.

Based on the optimized model parameters of the Si–P system, the iso-activity contours of P in a pure liquid standard state between 1000 K and 2200 K are predicted in Figure 5. The iso-activities at $a_{\mathrm{P}\left(\mathrm{l}\right)}={1\times 10}^{-7}$, ${1\times 10}^{-6}$, ${1\times 10}^{-5}$, ${1\times 10}^{-4}$, ${1\times 10}^{-3}$, ${1\times 10}^{-2}$, and ${1\times 10}^{-1}$ are plotted for P(l). According to the present optimization, $a_{\mathrm{P}\left(\mathrm{l}\right)}$ has to be no more than ${1.99\times 10}^{-8}$ and ${1.87\times 10}^{-7}$ to control P in liquid Si within ${1\times 10}^{-5}\mathrm{w}\mathrm{t}.\%$ at 1773 K and ${1\times 10}^{-4}\mathrm{w}\mathrm{t}.\%$ at 1723 K, respectively.

3.2. The Fe–P and Fe–Si Systems

The Fe–P and Fe–Si systems were recently optimized by the present authors [28,29]. In the modeling, gas phase, liquid, BCC_A2, FCC_A1, Fe₃P, Fe₂P, FeP, and FeP₂ of the Fe–P system and liquid, BCC_A2, BCC_B2, FCC_A1, solid Si, Fe₂Si, Fe₅Si₃, FeSi, FeSi₂, and Fe₃Si₇ of the Fe–Si system were taken into account. The optimized model parameters of all these phases can be found elsewhere [28,29] and were adopted in this study.

3.3. The Si–Fe–P Systems

Based on thermodynamic descriptions of the binary Fe-P, Fe–Si, and Si–P systems, gas phase, red P, liquid solution, solid solutions including FCC_A1, BCC_A2, BCC_B2, and solid Si, and stoichiometric compounds including SiP, SiP₂, FeP₂, Fe₂Si, Fe₅Si₃, FeSi, FeSi₂, and Fe₃Si₇ were treated as stable phases in the Si–Fe–P system. Moreover, the dissolution of Si in Fe₃P and Fe₂P to form Me₃P and Me₂P solid solutions in the formulas of (Fe)₃(P, Si) and (Fe)₂(P, Si), which have been confirmed in the experiments [68,69,70], were also considered in current thermodynamic modeling. Meanwhile, the MeP solid solution in the formula of (Fe)(P, Si), which was assumed in the assessment of Miettinen and Vassilev-Urumov [33], and the ternary FeSi₄P₄ compound [71,72,73] were taken into account as well. The model parameters of all binary sub-systems were combined to describe the ternary Si–Fe–P system.

3.3.1. Phase Diagram

The phase equilibria for various vertical sections of Si–Fe–P alloys were measured by Vogel and Giessen [71] and Hummitzsch and Sauerwald [74] using thermal analysis, microscopic examination, and chemical analysis. In the former experiments [71], a full composition range of FeSi alloys containing P up to $\mathrm{w}\mathrm{t}.\%\mathrm{P}=32.5$ was used, and a ternary stoichiometric compound FeSi₄P₄ was found to melt at 1483 K. The experimental data for pseudobinary diagrams of FeP–FeSi, Fe₂P–FeSi, Fe₃P–FeSi, FeSi–Fe₄Si₄P₄, FeSi₂–FeSi₄P₄, and isopleths of $\mathrm{w}\mathrm{t}.\%\mathrm{P}=13,8,5$ and $\mathrm{w}\mathrm{t}.\%\mathrm{S}\mathrm{i}=7$ are compared with the previous and present calculations in Figure 6, Figure 7 and Figure 8. Basically, the calculation results of Miettinen and Vassilev-Urumov [33] are in less satisfactory agreement with most of the phase diagram data. The discrepancies in the eutectic reaction $\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{M}\mathrm{e}\mathrm{P}+\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)$ of the FeP–FeSi section and solidus boundaries of the Fe₃P–FeSi section have been resolved by the present study, as shown in Figure 6a,c, with careful optimization in the Gibbs energies of liquid, MeP, Me₂P, and Me₃P phases. Moreover, the ${∆H}_{298.15\mathrm{K}}^{°}$ and $S_{298.15\mathrm{K}}^{°}$ of FeSi₄P₄ compound were determined to be $-\mathrm{33,4600}$ J/mol and $175$ J/(mol·K) to reproduce its melting point [71] and improve the phase equilibria of the isopleth at $\mathrm{w}\mathrm{t}.\%\mathrm{P}=13$, as shown in Figure 7 and Figure 8a. In particular, the present optimization shows significant improvement in the liquidus and solidus boundaries of the low-Si region for $\mathrm{w}\mathrm{t}.\%\mathrm{P}=5$ and $\mathrm{w}\mathrm{t}.\%\mathrm{P}=8$ isopleths, compared to the previous assessment [33], as shown in Figure 8b,c.

The solubility of P in α-Fe (BCC_A2) with the addition of Si up to 4 wt.% at 1273 K was investigated by Kaneko et al. [75] using chemical analysis and X-ray diffraction. Their experimental data are compared to the present calculation in Figure 9. The solubility data are accurately reproduced using one interaction parameter $L_{\mathrm{S}\mathrm{i},\mathrm{P}:\mathrm{V}\mathrm{a}}^{\mathrm{B}\mathrm{C}\mathrm{C}\_\mathrm{A}2}=-\mathrm{52,300}$ J/mol, as shown in Table 1. Based on the present optimization, the solubility of P in Fe–Si alloys equilibrated Me₃P was calculated to decrease from 2.28 wt.% to 1.13 wt.% with added Si increasing up to 5 wt.%, as shown in Figure 9.

The Ni₃P-type iron phosphide containing Si up to 2.1 wt.% was detected in heat-treated Si steel at 1073 K by Kaneko et al. [76] using the XRD and electrolytic separation method. Figure 10 shows the calculated isothermal section of the Si–Fe–P phase diagram at 1073 K in comparison with experimental data [75,76]. Based on the present optimization, a maximum solubility of $\mathrm{w}\mathrm{t}.\%\mathrm{S}\mathrm{i}=4.5$ in Fe₃P through substitution of P at 1073 K was calculated using an additional model parameter $G_{\mathrm{F}\mathrm{e}:\mathrm{S}\mathrm{i}}^{{\mathrm{M}\mathrm{e}}_3\mathrm{P}}=3G_{\mathrm{F}\mathrm{e}\left(\mathrm{B}\mathrm{C}\mathrm{C}\right)}^{°}+G_{\mathrm{S}\mathrm{i}(\mathrm{D}\mathrm{i}\mathrm{a}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{d}\_\mathrm{A}4)}^{°}-94,500$ J/mol, as listed in Table 1. On the other hand, Me₂P and MeP in the formulas of (Fe)₂(P, Si) and (Fe)(P, Si) were optimized to be complete solid solutions to reproduce other phase diagram data in Figure 8, Figure 9 and Figure 10.

The predicted liquidus surface projection of the Si–Fe–P system between 1073 K and 1873 K is presented in Figure 11, along with experimental data of Vogel and Giesson [71]. The invariant reactions are summarized in

Table 3. Invariant reactions of the ternary Si–Fe–P system with experimental data.

Code	Invariant Reactions	Wt.%Fe	Wt.%Si	Wt.%P	T, K
E1	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{M}\mathrm{e}\mathrm{P}+{\mathrm{M}\mathrm{e}}_2\mathrm{P}+\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)$	71.93	9.90	18.17	1438 [*]
		71.00	14.00	15.00	1439 [71]
E2	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)+\mathrm{M}\mathrm{e}\mathrm{P}$	48.58	29.38	22.04	1366 [*]
		47.10	34.00	18.90	1368 [71]
E3	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}}_3{\mathrm{S}\mathrm{i}}_7$(s)	43.37	38.99	17.64	1371 [*]
		46.00	36.40	17.60	1369 [71]
E4	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(s\right)+\mathrm{S}\mathrm{i}$	2.55	61.63	35.82	1400 [*]
		6.00	61.50	32.50	1389 [71]
E5	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}\mathrm{P}\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)+{\mathrm{S}\mathrm{i}\mathrm{P}}_2$(s)	3.92	38.57	57.51	1428 [*]
E6	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}={\mathrm{M}\mathrm{e}}_2\mathrm{P}+{\mathrm{M}\mathrm{e}}_3\mathrm{P}+\mathrm{B}\mathrm{C}\mathrm{C}\_\mathrm{B}2$	84.93	7.84	7.23	1335 [*]
E7	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}={\mathrm{F}\mathrm{e}}_2\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)+{\mathrm{M}\mathrm{e}}_2\mathrm{P}+\mathrm{B}\mathrm{C}\mathrm{C}\_\mathrm{B}2$	81.01	14.19	4.80	1358 [*]
E8	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{R}\mathrm{e}\mathrm{d}\mathrm{P}+{\mathrm{S}\mathrm{i}\mathrm{P}}_2\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{P}}_2\left(\mathrm{s}\right)$	0.075	0.005	99.92	852 [*]
U1	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+\mathrm{S}\mathrm{i}={\mathrm{F}\mathrm{e}}_3{\mathrm{S}\mathrm{i}}_7\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4$(s)	34.87	47.28	17.85	1383 [*]
EX	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}+{\mathrm{F}\mathrm{e}}_3{\mathrm{S}\mathrm{i}}_7\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(s\right)$	33.50	49.20	17.30	1386 [71]
U2	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)={\mathrm{M}\mathrm{e}}_2\mathrm{P}+{\mathrm{F}\mathrm{e}}_2\mathrm{S}\mathrm{i}\left(\mathrm{s}\right)$	80.03	15.21	4.76	1365 [*]
U3	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+\mathrm{B}\mathrm{C}\mathrm{C}\_\mathrm{A}2={\mathrm{M}\mathrm{e}}_3\mathrm{P}+\mathrm{B}\mathrm{C}\mathrm{C}\_\mathrm{B}2$	86.54	5.92	7.54	1342 [*]
$\mathbf{U}4$	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+\mathrm{M}\mathrm{e}\mathrm{P}=\mathrm{F}\mathrm{e}{\mathrm{P}}_2\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(s\right)$	39.22	12.97	47.81	1410 [*]
U5	$\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)={\mathrm{S}\mathrm{i}\mathrm{P}}_2\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{P}}_2\left(s\right)$	19.77	5.32	74.91	1336 [*]

* optimized in the present study.

. The liquidus isothermals from 1773 K to 1073 K were plotted in colorful lines. The calculated invariant points E1, E2, E3, and E4 agree reasonably with the data within experimental errors. However, the eutectic reaction $\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}=\mathrm{S}\mathrm{i}+{\mathrm{F}\mathrm{e}}_3{\mathrm{S}\mathrm{i}}_7\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)$ in point EX proposed by Vogel and Giessen [71] was supposed to be a mistake and has been corrected to the peritectic reaction $\mathrm{L}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}+\mathrm{S}\mathrm{i}={\mathrm{F}\mathrm{e}}_3{\mathrm{S}\mathrm{i}}_7\left(\mathrm{s}\right)+{\mathrm{F}\mathrm{e}\mathrm{S}\mathrm{i}}_4{\mathrm{P}}_4\left(\mathrm{s}\right)$, which is labeled as U1 in Figure 11.

3.3.2. Thermodynamic Properties of the Si–Fe–P Liquid Solution

The optimized Gibbs energies of binary Si–P, Fe-P, and Fe–Si liquid solutions were interpolated using a “Toop-like” approximation (Fe as the asymmetric component) to describe the ternary Si–Fe–P liquid solution, as discussed in Section 2.4. Moreover, three small MQM parameters, as given in Table 1, are still necessary to reproduce the phase diagram and thermodynamic property data simultaneously.

The solubility of P in a wide composition range of molten Si–Fe alloys was measured at 0.163~0.184 Pa of P₂(g) pressure and 1723 K by Ueda et al. [77] using the transportation method and chemical analysis. Their experimental data, together with the P solubility data in Si(l) by Miki et al. [66], are compared to the previous assessment and present optimization results in Figure 12. As shown in the figure, the solubility of P in liquid Si decreases slightly first and then increases sharply with the increase of Fe content, and the minimum P solubility was calculated to be $\mathrm{w}\mathrm{t}.\%\mathrm{P}=0.0212$ at $\mathrm{w}\mathrm{t}.\%\mathrm{F}\mathrm{e}=50.96$ from the present study. This turning is due to a maximized negative interaction between Fe and Si in the liquid solution. The calculation results from the assessment of Miettinen and Vassilev-Urumov [33] agree well with the experimental data at $\mathrm{w}\mathrm{t}.\%\mathrm{F}\mathrm{e}<64$ but deviate largely from the higher-Fe data. This discrepancy has been successfully resolved based on the present thermodynamic optimization using the MQM, as shown in Figure 12.

Figure 13 shows the calculated natural logarithm activity coefficient of P affected by Si ($\mathrm{l}\mathrm{n}{\gamma }_{\mathrm{P}}^{\mathrm{S}\mathrm{i}}$) in Fe-based Fe–Si–P liquid solution at 1673 K and 1873 K compared to experimental data [78,79,80]. Yamada and Kato [78] investigated the activity coefficient of P in various molten Fe–Si–P alloys at 1873 K using the Knudsen effusion method. The Si content ranged from $1\mathrm{w}\mathrm{t}.\% \mathrm{t}\mathrm{o} 7\mathrm{w}\mathrm{t}.\%,$ but the P content was maintained at $1\mathrm{w}\mathrm{t}.\%$ for all the samples. As a result, the activity coefficient interaction parameter was determined to be ${\epsilon }_{\mathrm{P}}^{\mathrm{S}\mathrm{i}}=11.9\pm 0.6$ at 1873 K. Ban-ya et al. [79] carried out transportation experiments to measure the vapor pressure of phosphorus above the Fe–Si–P melts ($\mathrm{w}\mathrm{t}.\%\mathrm{P}=6.2~13.2$) at 1673 K and calculated ${\epsilon }_{\mathrm{P}}^{\mathrm{S}\mathrm{i}}=7.68\pm 0.44$ for this temperature. The experimental results of Ban-ya et al. [79] were not favored because they assumed only P₂(g) in the gas phase. However, the vaporization of Fe and the formation of other gas species, such as P(g) and P₄(g), at such conditions cannot be neglected. According to the present calculations, the partial pressure of Fe(g), P₂(g), P(g), and P₄(g) are ${2.11\times 10}^{-6}$ atm, ${2.84\times 10}^{-6}$ atm, ${3.79\times 10}^{-8}$ atm and ${1.55\times 10}^{-12}$ atm respectively at $\mathrm{w}\mathrm{t}.\%\mathrm{P}=9$, $\mathrm{w}\mathrm{t}.\%\mathrm{S}\mathrm{i}=6(x_{Si}=0.105)$ and 1673 K. Schenck et al. [80] studied thermodynamic behavior of P in high-P Fe–Si–P alloys ($\mathrm{w}\mathrm{t}.\%\mathrm{P}=14.6~29.8$) equilibrated with the P₂(g) gas at 1788 K using the X-ray fluorescence and chemical analysis. They proposed an average of ${\epsilon }_{\mathrm{P}}^{\mathrm{S}\mathrm{i}}=14.2$ for 1788 K based on the vapor pressure data, which are significantly scattered, as shown in Figure 13, and not considered in the present thermodynamic modeling. The experimental data of Yamada and Kato [78] were adopted to determine the thermodynamic properties of the Fe–Si–P liquid solution. ${\epsilon }_{\mathrm{P}}^{\mathrm{S}\mathrm{i}}$ was calculated to be 11.09 at 1873 K and 14.05 at 1673 K from the optimized Fe–Si–P database.

3.3.3. Predicted Phase Diagram of the Si–Fe–P System

Based on the optimized model parameters for the Si–Fe–P system, the iso-activity contours of Si(l), Fe(l), and P(l) in pure liquid standard state at 1873 K are predicted in Figure 14. As can be seen in the figure, the iso-activities at ${1\times 10}^{-4}$, ${1\times 10}^{-3}$, ${1\times 10}^{-2}$, ${1\times 10}^{-1}$, 0.3, 0.5, 0.7, and 0.9 are plotted for Si(l) and Fe(l) and the iso-activities at ${1\times 10}^{-8}$, ${1\times 10}^{-7}$, ${1\times 10}^{-6}$, ${1\times 10}^{-5}$ ${1\times 10}^{-4}$, ${1\times 10}^{-3}$, ${1\times 10}^{-2},$ and ${1\times 10}^{-1}$ for P(l).

4. Summary

The Si–P and Si–Fe–P systems in the entire composition range were thermodynamically modeled based on the critical evaluation of all available experimental data. The liquid phases and solid solutions were described using the Modified Quasichemical Model (MQM) and Compound Energy Formalism (CEF), respectively. The liquid solution, solid Si, SiP, and SiP₂ were reoptimized to resolve the discrepancies left in previous assessments of the Si–P system. Moreover, the Gibbs energies of Me₃P, Me₂P, and MeP solid solutions with substitution of P in Fe₃P, Fe₂P, and FeP with Si, respectively, and FeSi₄P₄ compound were well determined to reproduce the phase diagram data more accurately. According to the present optimization, a consistent and accurate thermodynamic database of the Si–Fe–P system has been constructed and used to predict unexplored thermodynamic properties and phase diagrams. The present database can be applied to process optimization of Si refining and alloy design.