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

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)3(P,Si)1, (Fe)2(P,Si)1, and (Fe)1(P,Si)1 solid solutions and FeSi4P4 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.


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 × 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 where x i is the mole fraction of species i, G • 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 θ is the atmospheric pressure (=1 atm), and f is the gas fugacity and is identical to the gas pressure (in atm) at normal pressure.

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 3 P, Fe 2 P, FeP, and FeP 2 of the Fe-P system, Fe 2 Si, Fe 5 Si 3 , FeSi, FeSi 2 , and Fe 3 Si 7 of the Fe-Si system, SiP, SiP 2 of the Si-P system, and FeSi 4 P 4 of the Si-Fe-P system were determined from their corresponding heat capacity C P (J/(mol·K)), standard enthalpy of formation ∆H In the cases of pure elements and stoichiometric compounds exhibiting magnetic behavior, an additional Gibbs energy of magnetic contribution G mg (J/mol) will be applied. In this study, G mg was applied to Fe (BCC_A2, FCC_A1), Fe 3 P, and Fe 5 Si 3 and determined from an empirical expression proposed by Inden [36] and modified by Hillert and Jarl [37]: where τ is expressed by T/T * , T * is the Curie temperature T C (K) for ferromagnetic ordering or the Néel temperature T N (K) for anti-ferromagnetic ordering, g(τ) 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).

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 3 P, Fe 2 P, and FeP to form (Fe) 3 (P, Si) 1 , (Fe) 2 (P, Si) 1 , and (Fe) 1 (P, Si) 1 solid solutions represented by Me 3 P, Me 2 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].

FCC_A1 and Solid Si Solutions
FCC_A1 and solid Si (diamond_A4) solutions were described with the formula (Fe, Si, P) 1 (Va) 1 , and their Gibbs energies per formula unit were calculated as follows: ..
x Fe x Si x P L q Fe,Si,P + G mg (4) where x 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 m Fe,P , L k Si,P , L p Fe,Si , and L q Fe,Si,P are adjustable interaction parameters of corresponding binary and ternary systems (J/mol), and G mg is the magnetic contribution to the Gibbs energy (J/mol).

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) 1 (Va) 3 and the ordered part with the formula (Fe, Si, P) 0.5 (Fe, Si, P) 0.5 (Va) 3 . 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: G order BCC (y i , y j ) = y Fe y Fe G Fe:Fe + y Si y Si G Si:Si + y P y P G P:P +y Fe y Si G Fe:Si + y Si y Fe G Si:Fe + y Fe y P G Fe:P +y P y Fe G P:Fe + y Si y P G Si:P + y P y Si G P:Si +0.5RT y Fe Iny Fe + y Si Iny Si + y P Iny P +0.5RT y Fe Iny Fe + y Si ln y Si + y P ln y P where i, j, and k represent Fe, Si, P. y i , y j , y k and y i , y j , y 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) 3 , respectively. The Gibbs energy of the BCC solid solution combining both ordered and disordered contributions was determined from Equation (7): when the site fractions of component i in the first sublattice are equal to that in the second sublattice (y i = y i ), then the ordering contribution ∆G order BCC equals nil, and the Gibbs energy of the BCC phase is the same as that of the disordered BCC_A2 (G disorder S ) calculated by Equation (4). In the case of y i = y i , then ordering contribution ∆G order BCC becomes negative, and the Gibbs energy of the BCC_B2 phase can be calculated using Equation (7). The solid solutions Me 3 P, Me 2 P, and MeP were described with the formula (Fe) n (P, Si) 1 , where n = 3, 2, 1 for Me 3 P, Me 2 P, and MeP, respectively. Their Gibbs energies per formula unit were calculated based on the CEF [38] as follows: Me n P = y P G • Fe n P + y Si G Fe n Si + RT(y P ln y P + y Si ln y Si ) + ∑ m=0,1,2... y P y Si L m Fe:P,Si + G mg (8) here, G • Fe 3 P , G • Fe 2 P , and G • FeP are optimized Gibbs energies (J/mol) of stoichiometric Fe 3 P, Fe 2 P, and FeP compounds of the Fe-P system [28]; G Fe 3 Si , G Fe 2 Si , and G FeSi are Gibbs energies (J/mol) of Fe 3 Si, Fe 2 Si, and FeSi combinations respectively, which need to be optimized in this work; y P and y Si are site fractions of P and Si in the second sublattice, respectively; L m Fe:P,Si is the adjustable interaction parameter (J/mol), and G mg is the magnetic contribution to the Gibbs energy (J/mol).

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: 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 ∆g AB 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): where n A and n B are the mole numbers of A atoms and B atoms (mol), G • A and G • B are the molar Gibbs energies of pure A and B in a liquid state (J/mol), and ∆S conf. AB is the configurational entropy of mixing (J/(mol·K)) given by Equation (11): where n AA , n BB , and n AB represent the mole numbers of (A-A), (B-B) and (A-B) pairs (mol), X AA , X AB , and X AB are pair fractions of corresponding atom pairs, and Y A and Y B are coordination equivalent fractions of A atoms and B atoms. The pair fractions X AA , X BB , X AB and coordination equivalent fractions Y A , Y B can be calculated using Equations (12)-(16): ∆g AB 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 AA and X BB .
where ∆g • AB , g i0 AB and g 0j AB 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 A and Z B , can be varied with the composition to reproduce the short-range ordering.
here Z A AA is the value Z A when all nearest neighbors of an A atom are A atoms, and Z A AB is the value of Z A when all nearest neighbors of the A atom are B atoms. Z B BB and Z B BA are defined in an analogous manner. In the present study, Z Fe FeFe = Z Si SiSi = Z P PP = 6 [28,29], [28][29][30], and Z Fe FeP = 3 [28], as listed in Table 1. Table 1. Optimized model parameters for the Si-Fe-P system. Heat capacity C P (J/(mol·K)), standard enthalpy of formation ∆H

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.

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 (dia-mond_A4), red P, SiP, and SiP 2 are stable condensed phases in the Si-P system and are also adopted by the present study. Table 2. Summary of the crystal structure information of all solid phases of the Si-Fe-P system.
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 (wt.%P = 42.9) for the reaction liquid = Si + SiP(s) 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 liquid = Si + SiP(s) was modified to occur at wt.%P = 36.6 and 1404 K, and the melting points of SiP and SiP 2 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 2 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 2 could not be resolved [30,31]. The eutectic reaction Liquid = SiP(s) + SiP 2 (s) was reported by Ugai et al. [55] to occur at 1398 K, which is 53 K lower than the melting point of SiP 2 . 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 2 melting point. 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 SiP2 could not be resolved [30,31]. The eutectic reaction Liquid = SiP(s) + SiP ( ) was reported by Ugai et al. [55] to occur at 1398 K, which is 53 K lower than the melting point of SiP2. 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 SiP2 melting point.
(a) (b) Figure 1. The phase diagrams of the Si-P system (a) with suppression of gas phase and (b) with gas phase at a total pressure of 0.01, 0.1, 0.5, and 1 atm, compared to available experimental data [43][44][45][46][47][48][49][50][51][52][53][54][55]. For abbreviations, see Abbreviations Section. 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], Figure 1. The phase diagrams of the Si-P system (a) with suppression of gas phase and (b) with gas phase at a total pressure of 0.01, 0.1, 0.5, and 1 atm, compared to available experimental data [43][44][45][46][47][48][49][50][51][52][53][54][55]. For abbreviations, see Abbreviations Section. 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 Diamond_A4 Si,P:Va = −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 Si phase boundaries above the eutectic temperature by Jung and Z adopted with slight modification, while the boundaries below the eutect were pushed to the P-richer side using a much negative interact , : _ 50208 J/mol (Table 1) to reproduce the experimental dat [44], Kooi [46], and Soimi et al. [53], which are in good consistency, as sho The distribution coefficient of P in Si, represented by , , ⁄ concentration of P in solid and liquid Si, respectively), is a key factor in behavior of P during the crystallizing and melting of Si. The equilibrium to be 0.038 by Hall [56] based on conductivity measurements, 0.09 b tangstad [54] from the solidus and liquidus boundaries of Si-P alloys, and enko and Yudin [52] from the enthalpy change of P in solid and liquid Si p higher distribution coefficient ( 0.35 ) was originally obtained by St means of the Czochralski crystal growth and radiochemical analysis. This as the equilibrium distribution coefficient in various literature [44,[58][59][60][61], w mistaken as "new" sources of experimental data [5,54,62]. Huff et al. [62] i distribution coefficient of P in Si under different Czochralski crystal pul result, an "effective" distribution coefficient 0.32 at the pulling rate and 0.42 at the pulling rate of 1.8 mils/sec, were derived. Accordin [62], the effective is dependent on the crystal growth rate, impurity co terface orientation, temperature gradients, stirring conditions, etc. Recen determined the effective 0.31 and 0.33 from directional solidificatio at the rate of 2.08 10 m/s and 3.08 10 m/s, respectively. It is not values resulting from the Czochralski crystal growth method [57,62] and solidification method [6] are typically higher than those from conductivi chemical equilibrium measurements [52,54,56]. Because such Czochralski and directional solidification experiments [6,57,62] proceeded typically rium conditions depending on the Czochralski rate and temperature gr slight supercooling of dilute Si-P alloys can lead to distinct concentration The distribution coefficient of P in Si, represented by L P = C P,S /C P,L (C P,S and C P,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 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 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 P = 0.32 at the pulling rate of 1.1 mils/sec and L P = 0.42 at the pulling rate of 1.8 mils/sec, were derived. According to Huff et al. [62], the effective L 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 P = 0.31 and 0.33 from directional solidification ingots grown at the rate of 2.08 × 10 −6 m/s and 3.08 × 10 −6 m/s, respectively. It is noticeable that L 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 P values, which are "so-called" effective distribution coefficients, should be overestimated compared to equilibrium ones. The equilibrium L 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 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).

Thermodynamic Stability of Si Phosphides
The Gibbs energies of intermediate silicon phosphides SiP and SiP 2 , 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 2 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 −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 2 are available in the literature, so ∆H • 298.15K and S • 298.15K for SiP 2 were optimized to be −79, 950 J/mol and 64 J/(mol·K), respectively, in the present study to reproduce its melting point data [55] and make SiP 2 stable down to the ice point. The optimized thermodynamic properties of SiP and SiP 2 are given in Table 1.

Thermodynamic Stability of Si Phosphides
The Gibbs energies of intermediate silicon phosphides SiP and SiP2, culated from their heat capacity ( ), standard enthalpy of formation standard entropy ( . °) as expressed by Equation (2), were reoptimiz to improve the thermodynamic description of the Si-P system. The were directly taken from the assessment of Jung and Zhang [30]. Ugai mined . ° 34.78 J/(mol·K) for SiP through its low-temperatur value was adopted by the present study without modification. The stand formation for SiP was optimized to 64,000 J/mol to reproduce its me from Ugai et al. [55], as shown in Figure 1a, and SiP dissociation pressur presented in Figure 3. On the other hand, since no reliable Gibbs energy available in the literature, so ∆ . ° and . ° for SiP2 were o 79,950 J/mol and 64 J/(mol·K), respectively, in the present study to rep ing point data [55] and make SiP2 stable down to the ice point. The optim namic properties of SiP and SiP2 are given in Table 1.

Thermodynamic Properties of the Si-P Liquid Solution
Thermodynamic properties of the Si-P liquid solution at any designa and temperature can be predicted from the developed Si-P database. A demands super high purity with P content down to 10 wt.% level, so has been paid to the thermodynamic properties of the dilute region bec mary importance to the Si refining process. Based on the present thermod zation, the Henrian activity coefficient of P in Si(l), °, depending ture, was determined: where T is the temperature in Kevin (K). Based on the optimized °, energy for the dissolution of P2(g) into liquid Si (1 wt.% standard state) as follows: Figure 3. The optimized dissociation pressure (P 2 (g) + P 4 (g)) of SiP from the present study in comparison with experimental data [55,64,65]. For abbreviations, see Abbreviations Section.

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), γ • PinSi(l) , depending on the temperature, was determined: where T is the temperature in Kevin (K). Based on the optimized γ • PinSi(l) , the molar Gibbs energy for the dissolution of P 2 (g) into liquid Si (1 wt.% standard state) was determined as follows: 0.5P 2 (g) = P inSi(l) (wt.%); ∆G The calculated Gibbs energies of P 2 (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. The calculated Gibbs energies of P2(g) dissolution in Si(l) from previous the present study are plotted in Figure 4, along with the experimental d [66] measured using the transportation method and those of Zaitsev et Knudsen effusion mass spectrometry. These two sets of data show a disti about 15 kJ/mol. As has been pointed out [31], the activity data of Zaits less self-consistent. Thus, their Gibbs energy data from the same exper considered in the present thermodynamic modeling either. Instead, the ex of Miki et al. [66] were favored by the present study, which agrees with result of Liang and Vassilev-Urumov [31]. However, the modeling res Zhang [30] shows a large deviation from both sets of experimental data [ in Figure 4. Based on the optimized model parameters of the Si-P system, the tours of P in a pure liquid standard state between 1000 K and 2200 K Figure 5. The iso-activities at 1 10 , 1 10 , 1 10 , 1 1 10 , and 1 10 are plotted for P(l). According to the present op has to be no more than 1.99 10 and 1.87 10 to control P in 1 10 wt. % at 1773 K and 1 10 wt. % at 1723 K, respectively. 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 P(l) = 1 × 10 −7 , 1 × 10 −6 , 1 × 10 −5 , 1 × 10 −4 , 1 × 10 −3 , 1 × 10 −2 , and 1 × 10 −1 are plotted for P(l). According to the present optimization, a P(l) has to be no more than 1.99 × 10 −8 and 1.87 × 10 −7 to control P in liquid Si within 1 × 10 −5 wt.% at 1773 K and 1 × 10 −4 wt.% at 1723 K, respectively.

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 3 P, Fe 2 P, FeP, and FeP 2 of the Fe-P system and liquid, BCC_A2, BCC_B2, FCC_A1, solid Si, Fe 2 Si, Fe 5 Si 3 , FeSi, FeSi 2 , and Fe 3 Si 7 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.

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 2 , FeP 2 , Fe 2 Si, Fe 5 Si 3 , FeSi, FeSi 2 , and Fe 3 Si 7 were treated as stable phases in the Si-Fe-P system. Moreover, the dissolution of Si in Fe 3 P and Fe 2 P to form Me 3 P and Me 2 P solid solutions in the formulas of (Fe) 3 (P, Si) and (Fe) 2 (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 4 P 4 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.
Based on the optimized model parameters of the Si-P system, the tours of P in a pure liquid standard state between 1000 K and 2200 K Figure 5. The iso-activities at 1 10 , 1 10 , 1 10 , 1 1 10 , and 1 10 are plotted for P(l). According to the present op has to be no more than 1.99 10 and 1.87 10 to control P in 1 10 wt. % at 1773 K and 1 10 wt. % at 1723 K, respectively. Figure 5. Predicted iso-activity contours of the Si-P system for a P(l) = 1 × 10 −7 ∼ 1 × 10 −1 (pure P(l) as the reference state).

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 wt.%P = 32.5 was used, and a ternary stoichiometric compound FeSi 4 P 4 was found to melt at 1483 K. The experimental data for pseudobinary diagrams of FeP-FeSi, Fe 2 P-FeSi, Fe 3 P-FeSi, FeSi-Fe 4 Si 4 P 4 , FeSi 2 -FeSi 4 P 4 , and isopleths of wt.%P = 13, 8, 5 and wt.%Si = 7 are compared with the previous and present calculations in Figures 6-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 Liquid = MeP + FeSi(s) of the FeP-FeSi section and solidus boundaries of the Fe 3 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 2 P, and Me 3 P phases. Moreover, the ∆H • 298.15K and S • 298.15K of FeSi 4 P 4 compound were determined to be −33, 4600 J/mol and 175 J/(mol·K) to reproduce its melting point [71] and improve the phase equilibria of the isopleth at wt.%P = 13, as shown in Figures 7 and 8a. In particular, the present optimization shows significant improvement in the liquidus and solidus boundaries of the low-Si region for wt.%P = 5 and wt.%P = 8 isopleths, compared to the previous assessment [33], as shown in Figure 8b    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 , : _ 52,300 J/mol, as shown in Table 1. Based on the present optimization, the solubility of P in Fe-Si alloys equilibrated Me3P 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 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 BCC_A2 Si,P:Va = −52, 300 J/mol, as shown in Table 1. Based on the present optimization, the solubility of P in Fe-Si alloys equilibrated Me 3 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 3 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 wt.%Si = 4.5 in Fe 3 P through substitution of P at 1073 K was calculated using an additional model parameter G  Table 1. On the other hand, Me 2 P and MeP in the formulas of (Fe) 2 (P, Si) and (Fe)(P, Si) were optimized to be complete solid solutions to reproduce other phase diagram data in Figures 8-10. aterials 2023, 16, x FOR PEER REVIEW Figure 9. Isothermal diagram of the Si-Fe-P system on the Fe-rich corner at 1273 experimental data [75]. For abbreviations, see Abbreviations Section. The Ni3P-type iron phosphide containing Si up to 2.1 wt.% was de treated Si steel at 1073 K by Kaneko et al. [76] using the XRD and electrol method. Figure 10 shows the calculated isothermal section of the Si-Fe-P at 1073 K in comparison with experimental data [75,76]. Based on the pre tion, a maximum solubility of wt. %Si 4.5 in Fe3P through substitution was calculated using an additional model parameter : _ ° 94,500 J/mol, as listed in Table 1. On the other hand, Me the formulas of (Fe)2(P, Si) and (Fe)(P, Si) were optimized to be complete to reproduce other phase diagram data in Figures 8-10.  . Isothermal diagram of the Si-Fe-P system on the Fe-rich corner at 1273 K, compared to experimental data [75]. For abbreviations, see Abbreviations Section. Figure 9. Isothermal diagram of the Si-Fe-P system on the Fe-rich corner at 1273 experimental data [75]. For abbreviations, see Abbreviations Section. The Ni3P-type iron phosphide containing Si up to 2.1 wt.% was de treated Si steel at 1073 K by Kaneko et al. [76] using the XRD and electro method. Figure 10 shows the calculated isothermal section of the Si-Fe-P at 1073 K in comparison with experimental data [75,76]. Based on the pr tion, a maximum solubility of wt. %Si 4.5 in Fe3P through substitution was calculated using an additional model parameter : _ ° 94,500 J/mol, as listed in Table 1. On the other hand, Me the formulas of (Fe)2(P, Si) and (Fe)(P, Si) were optimized to be complete to reproduce other phase diagram data in Figures 8-10. The predicted liquidus surface projection of the Si-Fe-P system betw 1873 K is presented in Figure 11, along with experimental data of Vogel an The invariant reactions are summarized in Table 3. The liquidus isotherm to 1073 K were plotted in colorful lines. The calculated invariant points E1 agree reasonably with the data within experimental errors. However, the e Liquid Si Fe Si s FeSi P s in point EX proposed by Vogel and G supposed to be a mistake and has been corrected to the peritectic reactio Fe Si s FeSi P s , which is labeled as U1 in Figure 11. 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. 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 Liquid = Si + Fe 3 Si 7 (s) + FeSi 4 P 4 (s) in point EX proposed by Vogel and Giessen [71] was supposed to be a mistake and has been corrected to the peritectic reaction Liquid + Si = Fe 3 Si 7 (s) + FeSi 4 P 4 (s), which is labeled as U1 in Figure 11. Figure 11. Calculated liquidus surface projection of the Fe-Si-P system between 1073 K and 1873 K, compared to experimental data [71]. For abbreviations, see Abbreviations Section.

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 P2(g) pressure and 1723 K by Ueda et al. [77] using the transportation Figure 11. Calculated liquidus surface projection of the Fe-Si-P system between 1073 K and 1873 K, compared to experimental data [71]. For abbreviations, see Abbreviations Section.

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 2 (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 wt.%P = 0.0212 at wt.%Fe = 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 wt.%Fe < 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.
Materials 2023, 16, x FOR PEER REVIEW method and chemical analysis. Their experimental data, together with the P s data in Si(l) by Miki et al. [66], are compared to the previous assessment and pr timization results in Figure 12. As shown in the figure, the solubility of P in liqu creases slightly first and then increases sharply with the increase of Fe content minimum P solubility was calculated to be wt. %P 0.0212 at wt. %Fe 50.96 present study. This turning is due to a maximized negative interaction between in the liquid solution. The calculation results from the assessment of Miett Vassilev-Urumov [33] agree well with the experimental data at wt. %Fe 64 bu largely from the higher-Fe data. This discrepancy has been successfully resolved the present thermodynamic optimization using the MQM, as shown in Figure 12 [66,77]. For abbreviations, see Abbreviations Figure 13 shows the calculated natural logarithm activity coefficient of P aff Si (ln ) in Fe-based Fe-Si-P liquid solution at 1673 K and 1873 K compared t mental data [78][79][80]. Yamada and Kato [78] investigated the activity coefficient of ious molten Fe-Si-P alloys at 1873 K using the Knudsen effusion method. The S ranged from 1 wt. % to 7 wt. %, but the P content was maintained at 1 wt. % f samples. As a result, the activity coefficient interaction parameter was determin 11.9 0.6 at 1873 K. Ban-ya et al. [79] carried out transportation experi measure the vapor pressure of phosphorus above the Fe-Si-P melts (wt. %P 6 at 1673 K and calculated 7.68 0.44 for this temperature. The experiment of Ban-ya et al. [79] were not favored because they assumed only P2(g) in the g However, the vaporization of Fe and the formation of other gas species, such as P4(g), at such conditions cannot be neglected. According to the present calcula partial pressure of Fe(g), P2(g), P(g), and P4(g) are 2.11 10 atm, 2.84 1 3. 79 10 atm and 1.55 10 atm respectively at wt. %P 9 , wt. %Si 0.105 and 1673 K. Schenck et al. [80] studied thermodynamic behavior of P in h Si-P alloys (wt. %P 14.6~29.8) equilibrated with the P2(g) gas at 1788 K using fluorescence and chemical analysis. They proposed an average of 14.2 fo based on the vapor pressure data, which are significantly scattered, as shown in F and not considered in the present thermodynamic modeling. The experimenta Yamada and Kato [78] were adopted to determine the thermodynamic propert Fe-Si-P liquid solution.
was calculated to be 11.09 at 1873 K and 14.05 at 167 Figure 12. Calculated P solubility in various molten Si-Fe alloys at P P 2 (g) = 0.163 ∼ 0.184 Pa and T = 1723 K, compared to experimental data [66,77]. For abbreviations, see Abbreviations Section. Figure 13 shows the calculated natural logarithm activity coefficient of P affected by Si (lnγ Si P ) 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 1wt.% to 7wt.%, but the P content was maintained at 1wt.% for all the samples. As a result, the activity coefficient interaction parameter was determined to be ε Si P = 11.9 ± 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 ( wt.%P = 6.2 ∼ 13.2) at 1673 K and calculated ε Si P = 7.68 ± 0.44 for this temperature. The experimental results of Ban-ya et al. [79] were not favored because they assumed only P 2 (g) in the gas phase. However, the vaporization of Fe and the formation of other gas species, such as P(g) and P 4 (g), at such conditions cannot be neglected. According to the present calculations, the partial pressure of Fe(g), P 2 (g), P(g), and P 4 (g) are 2.11 × 10 −6 atm, 2.84 × 10 −6 atm, 3.79 × 10 −8 atm and 1.55 × 10 −12 atm respectively at wt.%P = 9, wt.%Si = 6(x Si = 0.105) and 1673 K. Schenck et al. [80] studied thermodynamic behavior of P in high-P Fe-Si-P alloys ( wt.%P = 14.6 ∼ 29.8) equilibrated with the P 2 (g) gas at 1788 K using the X-ray fluorescence and chemical analysis. They proposed an average of ε Si P = 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. ε Si P was calculated to be 11.09 at 1873 K and 14.05 at 1673 K from the optimized Fe-Si-P database.
Materials 2023, 16, x FOR PEER REVIEW Figure 13. Effect of Si on the activity coefficient of P ( ) in molten Fe-Si-P alloys from 1873 K, compared to experimental data [78][79][80]. For abbreviations, see Abbreviations Sect

Predicted Phase Diagram of the Si-Fe-P System
Based on the optimized model parameters for the Si-Fe-P system, the iso contours of Si(l), Fe(l), and P(l) in pure liquid standard state at 1873 K are pre Figure 14. As can be seen in the figure, the iso-activities at 1 10 , 1 10 , 1 10 , 0.3, 0.5, 0.7, and 0.9 are plotted for Si(l) and Fe(l) and the iso-activities at 1 10 , 1 10 , 1 10 1 10 , 1 10 , 1 10 , and 1 10 for P( (a) (b) Figure 13. Effect of Si on the activity coefficient of P (γ Si P ) in molten Fe-Si-P alloys from 1673 K to 1873 K, compared to experimental data [78][79][80]. For abbreviations, see Abbreviations Section.

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 10 , 1 10 , 1 10 , 1 10 , 0.3, 0.5, 0.7, and 0.9 are plotted for Si(l) and Fe(l) and the iso-activities at 1 10 , 1 10 , 1 10 , 1 10 1 10 , 1 10 , 1 10 , and 1 10 for P(l).

Summary
The Si-P and Si-Fe-P systems in the entire composition range were thermodyna cally modeled based on the critical evaluation of all available experimental data. The uid phases and solid solutions were described using the Modified Quasichemical Mo (MQM) and Compound Energy Formalism (CEF), respectively. The liquid solution, s Si, SiP, and SiP2 were reoptimized to resolve the discrepancies left in previous assessme of the Si-P system. Moreover, the Gibbs energies of Me3P, Me2P, and MeP solid soluti with substitution of P in Fe3P, Fe2P, and FeP with Si, respectively, and FeSi4P4 compou were well determined to reproduce the phase diagram data more accurately. Accord to the present optimization, a consistent and accurate thermodynamic database of the Fe-P system has been constructed and used to predict unexplored thermodynamic pr erties and phase diagrams. The present database can be applied to process optimizat of Si refining and alloy design.

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 2 were reoptimized to resolve the discrepancies left in previous assessments of the Si-P system. Moreover, the Gibbs energies of Me 3 P, Me 2 P, and MeP solid solutions with substitution of P in Fe 3 P, Fe 2 P, and FeP with Si, respectively, and FeSi 4 P 4 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.