Pair Fraction Distribution of the One-Sublattice Modified Quasichemical Model—Application to the Cu-Fe-S Liquid Solution

Within the formalism of the one-sublattice modified quasichemical model an equation for the computation of pair fraction distributions is derived considering distinct model modifications reported in the literature. For generalization an equation system is presented to determine pair fraction distributions of multicomponent solutions. Consistency of the derived equation with the analogous relationship of former versions of the modified quasichemical model is demonstrated. The derived equation is applied to the Cu-Fe-S liquid solution and its binary subsystems. Pair fraction distributions calculated exemplarily at 1200 °C (1473 K) are used to interpret recently reported predictions of the high-temperature Cu-Fe-S system.


Introduction
The classical quasichemical model according to Guggenheim [1] and Fowler and Guggenheim [2] has been modified by Pelton and Blander [3,4] in such a way that (1) the composition of maximum short-range order in a binary system can be chosen asymmetrically, (2) the energy of pair formation is allowed to be composition dependent, and (3) the model is applicable for multicomponent systems.
Considering a binary solution phase, the resulting modified quasichemical model (MQM) follows the concept of bonding of solution species A and B to form (A-B) pairs according to the pair-exchange reaction Using the thermodynamic equilibrium criterion of the Gibbs energy minimum being subject to the balance of amount of substance following equation can be derived in the framework of MQM: Equation 2 contains the equilibrium values of the three pair fractions X AA , X BB and X AB . They correspond to the ratios of moles of the (i-j) pairs n ij and the total amount of pairs according to X ij = n ij n AA þn BB þn AB ð Þ with i, j = A, B. The left side of Eq 2 resembles an equilibrium constant of reaction 1 which is according to Pelton and Blander [3] the reason why the model is called 'quasichemical.' At a certain temperature and pressure, the right side of Eq 2 represents a constant containing the Gibbs energy change Dg AB with respect to reaction (1) for the formation of two moles of (A-B) pairs from one mole of (A-A) and one mole of (B-B) pairs. Equation 2 determines the pair fraction distribution within MQM.
Further modifications and extensions to MQM were undertaken by Pelton et al., [5] Pelton and Chartrand, [6] Chartrand and Pelton, [7] and Pelton et al. [8] The first of these four articles deals with binary solutions whereas the second one refers to multicomponent solutions with one lattice or sublattice. The subsequent two studies focus on extensions to liquids with two 'sublattices'. According to the first two articles [5,6] two distinctive modifications of MQM in the pair approximation for short-range ordering in liquid and solid solutions result in a new model version denoted in this study as M'Q'M: (a) The coordination numbers are allowed to be composition dependent. (b) The Gibbs energy change Dg AB for the formation of (A-B) pairs of reaction (1) is expressed as a function of pair fractions instead of component fractions.
The purpose of this paper is the derivation of an equation analogous to Eq 2 with respect to M'Q'M, that is, considering both alterations (a) and (b). Generalization to multicomponent solutions with one sublattice is provided. Application of the derived equation to the Cu-Fe-S liquid solution is carried out: calculated pair fraction distributions are used to interpret promising predictions of the hightemperature Cu-Fe-S system recently presented by Waldner. [9] Validation of the derived equation is based on the binary subsystems, these are Cu-Fe [10] and the metal-sulfur systems Fe-S [11] and Cu-S [12] where M'Q'M was used for thermodynamic modeling of the liquid phase. The binary system Fe-S [11] was subject of a thermodynamic modeling study where M'Q'M was successfully applied for the first time to a metal-sulfur system.

Model Equations
The derivation of the equation ruling the pair fraction distribution within M'Q'M is developed for the ternary liquid solution of a system A-B-C. Thereby the generalization to multicomponent systems A-B-C-…-K-… on one hand and the reduction to a binary system A-B on the other hand is seen as straightforward.

Basic Equations of M'Q'M
According to Pelton et al. [5] and Pelton and Chartrand [6] the one-sublattice modified quasichemical model provides following total Gibbs energy of a ternary solution phase as a function of temperature, pressure and composition where the molar Gibbs energies g o A , g o B and g o C of the pure components are weighted with the numbers of moles of A, B and C denoted with n A , n B and n C .
The configurational entropy DS config describes random mixing of the (A-A), (B-B), (C-C), (A-B), (A-C), and (B-C) pairs in the one-dimensional Ising approximation: [5,6] The mole fractions X i and the coordination equivalent fractions Y i are defined according to with the coordination numbers Z i where i = A, B and C. The pair fraction X ij is given by the ratio number of moles of the (i-j) pairs n ij divided by the total number of moles of all occurring pairs with i, j = A, B and C as follows: The balance of amount of substance requires three relations between n A , n B and n C weighted by the corresponding coordination number and the mole amount of specific pairs: One distinctive alteration of M'Q'M compared to MQM is related to the coordination numbers which can vary with composition as follows: ðEq 12Þ where Z A AA stands for the coordination number of A when all nearest neighbors of A are exclusively A atoms, Z A AB represents the quantity when all nearest neighbors are B atoms and Z A AC when all nearest neighbors are C atoms. The other coordination numbers in Eq 12 and 13 are defined in the same way when B and C are the central atoms within their corresponding neighborhood.
Since the Gibbs energy function for the A-B-C liquid solution is based on a one-sublattice approach within the framework of M'Q'M first-nearest-neighbor interactions between A, B and C can be considered. The total excess Gibbs energy G ex as third term of Eq 3 contains the three types of first-nearest neighbor interactions in the ternary solution phase: The first nearest neighbor interaction between A and B species, expressed by the pair exchange reaction (A-A) ? (B-B) = 2(A-B) of Eq 1, is described by the term G ex AB of Eq 14. It contains the non-configurational Gibbs energy change Dg AB for the formation of two moles of (A-B) pairs according to reaction (1): The last two contributions of Eq 14, G ex AC and G ex BC stand for the other types of first-nearest neighbor interaction, e.g. between A and C species according to the pair exchange reaction: Their relationship to the corresponding non-configurational Gibbs energy changes Dg AC and Dg BC is analogous to Eq 15. The model quantities Dg AB , Dg AC and Dg BC may contain-if necessary-adjustable temperature and composition dependent parameters to be calculated in the framework of thermodynamic optimizations. Two examples of non-configurational Gibbs energy functions applied for real systems are given in chapter 3 by Eq 27 and 28 together with optimization data in Table 1. can be chosen as independent variables since the balance of amount of substance of Eq 8-10 together with Eq 11-13 allow substitution of n AA , n BB and n CC as follows: For finding the equilibrium distribution of the pair fractions following three partial derivatives of Eq 3 must be determined and set to zero: Introducing following ratios containing the coordination numbers of the Eq 11-13 as the partial derivatives of Eq 20-22 result in Table 1 Excess Gibbs energy parameters (in J mol -1 ) of the Cu-Fe-S liquid phase model according to binary optimization studies by Ansara and Janson, [10] Waldner and Pelton, [11] and Waldner [12] Dg o FeS = -104,888.10 ? 0.338 T Ref. [11] g 01 FeS = ? 35,043.32 -9.880 T Ref. [11] g 02 FeS = ? 23,972.27 Ref. [11] g 03 FeS = ? 30,436.82 Ref. [11] g 10 FeS = ? 8626.26 Ref. [11] g 20 FeS = ? 72,954.29 -26.178 T Ref. [11] g 40 FeS = ? 25,106 Ref. [11] Dg 0 CuS = -62,467.12 -11.071 T Ref. [12] g RT ln X a B AA X b A BB X AB Solving this system of the three Eq 23-25 allows the computation of the independent variables n AB , n AC and n BC . Consequently Eq 17-19 can then be used to determine n AA , n BB and n CC : Finally Eq 7 provides the equilibrium distribution of the six pair fractions at a distinct pressure, temperature and composition. Inserting these pair fractions e.g. into Eq (4) the configurational entropy can be calculated.

Generalization to Multicomponent Solutions
Generalization to a multicomponent solution A-B-C-…-K-…. with n system components and consequently n 2 À Á binary subsystems A-B, A-C,…, A-K,…, B-C,…, is accompanied with n relations due to the balance of amount of substance here noted e.g., for n AA , n BB and n KK and with n 2 À Á partial derivatives of the total Gibbs energy due to its minimization. Here the partial derivatives for the A-B, A-C and A-K subsystem are exemplarily noted of the resulting equation system as with the constants a K ¼ À 1 containing the coordination numbers of the subsystem A-K. Taking into consideration the corresponding Eq 5-7 extended for multicomponent solutions usage of this equation system yields the equilibrium distribution of n 2 À Á þ n Â Ã pair fractions of all possible (i-j) pairs at a distinct pressure, temperature and composition.

Reduction to Binary Solutions
Reduction of the equation system of the preceding section to a binary solution A-B results in Substitution of n AA and n BB according to

RT
Pelton and Blander [3] and Pelton et al. [5] pointed out, that the left side of Eq 2 serves as name-giver for the quasichemical model due to its analogy to the 'equilibrium' constant of the 'quasichemical reaction'

Application to the Cu-Fe-S Liquid Solution
Thermodynamic modeling of the high-temperature Cu-Fe-S liquid solution was recently carried in the framework of M'Q'M. [9] One single Gibbs energy expression at 1 bar total pressure describes the manifold nature of Cu-Fe-S liquid solutions from highly metallic via sulfur-rich to pure liquid sulfur. Predictions of sulfur potentials and hightemperature phase equilibria (demonstrated by four isothermal as well as four isopleth sections) agree with a large stock of experimental data available in the literature, solely on the basis of optimization studies of the binary subsystems Cu-Fe, [10] Fe-S [11] and Cu-S [12] without the use of adjustable ternary parameters. With respect to the extrapolation scheme from the binary subsystems into the ternary system two different options for the excess Gibbs energy offered by M'Q'M are thermodynamically analyzed. [9] Comparison of the model predictive ability favors an 'asymmetric' approach according to for the the non-configurational Gibbs energy change Dg CuS analogously to Dg AB of Eq 15. It can be seen that Eq 27 reduces to the excess Gibbs energy of the Cu-S binary system when no iron is present and consequently X CuFe and X FeFe are zero. The adjustable, and in case of necessity temperature dependent parameters Dg o CuS , g io CuS and g oj CuS were optimized in the thermodynamic modeling study of the Cu-S liquid phase [12] . The superscripts noted as a ring, i and j identify the parameters given in the study of Waldner [12] as follows: the parameter Dg o CuS contributes decidedly to the excess Gibbs energy in the composition region where (Cu-S) pairs are dominating. The parameter g io CuS was selected and optimized in the sulfur-rich region where the pair fraction of (S-S) significantly exceeds the pair fractions of (Cu-Cu) and (Cu-S) whereas the parameter g oj CuS is useful in copper-rich regions where the pair fraction of (Cu-Cu) exceeds significantly the pair fractions of (S-S) and (Cu-S).
According to Eq 14 the ternary Cu-Fe-S liquid solution contains two further excess Gibbs energy contributions with respect to the Fe-S and Cu-Fe binary subsystems. In an analogous manner as described for Dg CuS of Eq 27 the composition dependency of Dg FeS is based on corresponding parameters Dg o FeS , g io FeS and g oj FeS optimized in the thermodynamic modeling study of the Fe-S system. [11] M'Q'M allows the integration of binary subsystems optimized with a random-mixing Bragg-Williams and Redlich-Kister model where the following composition dependency of the third contribution to the excess Gibbs energy, Dg CuFe , is given as follows: In the Cu-Fe binary subsystem Eq 28 reduces to a 'Redlich-Kister' expression for which the adjustable model parameters i L CuFe were taken from Ansara and Jansson [10] . Table 1 lists all excess Gibbs energy parameters of Dg CuS , Dg FeS and Dg CuFe .

Pair Fraction Distribution of Binary Subsystems
Equation 26 is used to calculate the pair fraction distribution in each liquid phase of the binary Cu-S, Fe-S and Cu-Fe subsystems applying the data of Table 1. Figure 1 shows these three pair distributions over the complete composition range. The temperature is chosen with 1200°C (1473 K) according to the temperature at which predictions of sulfur potentials will be demonstrated here exemplarily. Equation 26 contains a partial derivative of the excess Gibbs energy G ex CuS ¼ n CuS =2 ð ÞDg CuS with the non-configurational Gibbs energy change Dg CuS of Eq 27. To give an example for such a partial quantity the derivative of G ex CuS with respect to n CuS is given where the derivation of the first two terms of Eq 27 are explicitly shown: The thermodynamic software package FactSage [13] where M'Q'M is implemented can be used to calculate pair fractions. Consequently, validation of Eq 26 is demonstrated by comparison with points (depicted as full diamonds) which are calculated with FactSage [13] on the basis of the Gibbs energy data from the optimization studies of the binary subsystems Cu-Fe [10] , Fe-S [11] and Cu-S [12] . It can be seen that full agreement between the computations based on FactSage [13] and Eq 26 exists.
The Cu-Fe liquid phase is described with a randommixing Bragg-Williams and Redlich-Kister model [10] as no short range ordering is taken into consideration. This is reflected in Fig. 1(c) with a much flatter maximum of Cu-Fe bonds in the pair fraction distribution of (Cu-Fe) bonds.
Metal-sulfur liquid solutions differ substantially in their internal chemical nature ranging from purely metallic via matte-like ('matte' is the metallurgical notation for a sulfur-rich liquid metal-sulfur phase) to pure liquid sulfur. Nevertheless, their thermodynamic description could be provided by one single Gibbs energy function over the complete composition range [12] . Consequently pair fraction distributions can also be calculated inside composition ranges where the liquid solution is part of heterogeneous phase equilibria. For example, the Cu-S system shows a miscibility gap at 1200°C (1473 K) where a copper-rich melt with 2.35 at. pct S coexists with a matte with 32.77 at. pct S [12] . From Fig. 1(a) it can be seen that X CuCu of this copper-rich melt is around 0.97 and X CuS of the matte phase is close to 0.92. Within the miscibility gap between 2.35 and 32.77 at. pct S the pair fraction X CuCu = 0.97 of the copper-rich melt should show a horizontal line according to Gibbs phase rule. The same is valid for the pair fraction X CuS = 0.92 of the matte phase which should also be constant from 32.77 to 2.35 at. pct S. However, instead of plotting two horizontal lines at X CuCu = 0.97 and X CuS = 0.92 between 2.35 and 32.77 at. pct S (and further horizontal lines for the other pair fractions) the lines outside the range of 2.35 and 32.77 at. pct S are continued smoothly by calculation the pair fraction distribution of a metastable single liquid phase. Inside heterogeneous phase regimes this procedure is conducted in this study for all other presentations of pair fraction distributions of the liquid phase.

Pair Fraction Distribution of the Ternary Cu-Fe-S System
Using the equation system 23-25 together with the data of Table 1 pair fraction distributions are calculated within the ternary Cu-Fe-S system around a selected ternary point P at 1200°C (1473 K). This point P corresponds to a real sample belonging to a large body of copper-iron-sulfur samples which were equilibrated with H 2 S/H 2 gas mixtures by Bale and Toguri [14] in the framework of a thermogravimetric technique for continuous quantitative analysis of sulfur partial pressures at 1200°C (1473 K). Figure 2 shows the published data from Bale and Toguri [14] together with a selected part of the data by Nagamori et al. [15] who also measured the vapor pressure of sulfur over Cu-Fe-S mattes by equilibrating them under H 2 S-H 2 gas mixtures at 1200°C (1473 K). Figure 2 reproduces a part of the predicted lines of the study by Waldner [9] for comparison with above selected but representative experimental data of sulfur potentials available from the literature. Using FactSage [13] as calculation tool experimental sulfur potentials are expressed as equilibrium pressures of S 2 gas. The experimental points [14,15] cover the composition range from the Cu-S binary through increasing ratios n Fe /(n Fe ? n Cu ) up to the Fe-S binary edge. The data according to the ratio n Fe /(n Fe ? n Cu ) = 0.549 from the study of Bale and Toguri [14] contain the sample P marked with a circle. The predictions of the liquid phase model (demonstrated as solid lines) are in satisfactory agreement with experimental data in the composition range between 1/3 (Cu-S edge) and 1/2 (Fe-S edge) mole fraction of sulfur Fig. 2 Predicted sulfur potentials (expressed as equilibrium pressure of S 2 , solid lines) of the Cu-Fe-S liquid solution at various ratios n Fe / (n Fe ? n Cu ) for 1200°C (1473 K) according to the study of Waldner [9] using FactSage [13] as calculation tool together with experimental data (ratio numbers given after symbols). The experimental point marked with P is selected to play an exemplary role for calculating ternary pair fraction distributions in this study and appears also in Fig. 3 where the sulfur potentials show a strong change over several orders of magnitude. To show the predicted stability of the liquid phase over the complete composition range gas phase formation at higher sulfur content is suppressed during the FactSage computation. Bale and Toguri [14] and Nagamori et al. [15] reported also experimental phase equilibria data which they derived from their measurements of sulfur potentials. Together with these data Fig. 3(a) demonstrates a Gibbs triangle for 1200°C (1473 K) with predicted phase boundaries of the study by Waldner [9] using FactSage [13] as calculation tool. Experimentally determined tie-lines are shown as solid lines of which their two end points are indicated by a certain symbol. Predicted tie-lines placed arbitrarily in two-phase regions are shown as dashed lines. The Gibbs triangle of Fig. 3(a) shows the location of the real sample P of the experimental data from Bale and Toguri [14] . This ternary sample P is depicted as circle. Its location in the  [9] using FactSage [13] as calculation tool together with experimental data. Experimental tie-lines are shown as solid lines; calculated tie-lines are shown as dashed lines. Abbreviation: fcc stands for the face-centeredcubic Cu-Fe-alloy phase. Two paths (solid lines) are calculated with Eq 23-25: Path e 1 -P-S along which the ratio Y Fe /Y Cu of P is kept constant; path e 2 -P-e 3 along which the pair fraction X SS of P remains unchanged. (b) Gibbs triangle of the Cu-Fe-S system at 1200°C (1473 K) with equivalent mole fractions as composition axis: Path e 1 -P-S along which the ratio Y Fe /Y Cu of P is kept constant; Path e 2 -P-e 3 along which the pair fraction X SS of P remains unchanged Gibbs triangle-not too close to a binary subsystem or one component corner-is seen as suitable to play an exemplary role for calculating pair fraction distributions in the ternary system of this study. Two paths emanating from P were calculated by solving the equation system 23-25 along which following conditions for composition variables are chosen to be constant: firstly, the equation system 23-25 was solved being subject to the condition of the constant ratio Y Fe /Y Cu of P which results in a path e 1 -P-S with the endpoints e 1 and the sulfur corner S. Secondly the condition of constant pair fraction X SS of P was set inducing the path e 2 -P-e 3 .
Alternatively, Fig. 3(b) demonstrates the corresponding Gibbs triangle when equivalent mole fractions Y i are used as axes instead of mole fractions x i with i = Cu, Fe and S. Due to the mathematical structure of Eq 28 it follows that the line of path e 1 -P-S is straightened. The obtained Gibbs triangle of Fig. 3(b) resembles now the well-known Kohler/Toop [16,17] asymmetric scheme being subject to the conditions of constant mole fraction of sulfur x S (corresponding to X SS in this study) and constant mole fraction ratio x Fe /x Cu (corresponding here to Y Fe /Y Cu ).  Fig. 4. Figure 5 reflects the fact that the interaction of iron produces stronger bonds than the interaction of copper with sulfur. The pair fraction distribution calculated being subject to the constraint of constant X SS of P at 1200°C (1473 K) is split into the two Fig. 6(a) and (b) as a function of mole fraction of iron. Figure 6(a) shows an almost linear course of the metal-sulfur pair fractions. Furthermore, the sum X FeS ? X CuS results in a nearly constant function of the metal-sulfur pair fractions. The fractions of the metalmetal pairs X CuCu , X CuCu and X CuFe are given in the enlargement of Fig. 6(b) since their values are an order of magnitude lower than the pairs with sulfur. The similarity of Fig. 6(b) with Fig. 1(c) originates in the fact that the 'isoplethal' plane along e 2 -P-e 3 in Fig. 3 may be interpreted as projection from the binary Cu-Fe subsystem into the ternary system. Figure 7 reproduces the three Gibbs energies of pair exchange reactions Dg CuS , Dg FeS and Dg CuFe at 1200°C (1473 K) along the path e 2 -P-e 3 of Fig. 3 as a function of mole fraction of iron. Due to the mathematical structure of Eq 27 and the analogous one for Dg FeS the constraint of constant X SS provides almost constant exchange reactions of metal-sulfur bonds over the complete composition range between the Cu-S and the Fe-S subsystems. In view of Fig. 5 it can again be seen that the interaction between the metals and sulfur generates stronger bonds than those between copper and iron.

Relationship between Pair Fraction Distributions and with Prediction of Sulfur Potentials
The predictive ability of M'Q'M using an asymmetric approach for the excess Gibbs energy [9] is interpretable by the pair fraction distributions of the Fig. 4 and 6 together with the Gibbs energies of pair reactions given in the Fig. 5 and 7. The appearance of the profound maxima in Fig. 4 correlates with an effective change of sulfur potentials within a narrow composition range due to consumption of sulfur by copper and iron species to form metal-sulfur bonds. The clearly negative values of Dg CuS and Dg FeS in the Fig. 5 and 7 including minima in Fig. 5 reveal that the thermodynamically preferred metal-sulfur pairs lower the activity of sulfur towards metal-rich regimes-as demonstrated in Fig. 2-sufficiently 'rapid' within a narrow composition range. This mechanism is effectively maintained not only for the ratio n Fe /(n Fe ? n Cu ) = 0.549 of point P but also for other ratios given in Fig. 2 due to the nearly constant pair exchange energies Dg CuS and Dg FeS depicted in Fig. 7.
In view of the concept of Toop [17] to treat certain system components differently due to their fundamental chemical disparity the nearly parallel course of the sum X FeS ? X CuS with constant X SS in Fig. 6(a) causes a second mechanism for the reliable predictive ability of M'Q'M using an asymmetric approach [9] : As long as (S-S) pairs remain constant the two metals do not vary strongly in their common capacity to form metal-sulfur pairs and occur almost mutually interchangeable.

Conclusions
Derivation of an equation ruling the pair fraction distribution within the framework of the one-sublattice modified quasichemical model [5,6] is feasible. The found equation can be generalized to be part of an equation system which allows the calculation of pair fraction distributions of multicomponent solution phases. It reduces consistently to the analogous name-giving ('quasichemical') relationship of former model versions when coordination numbers are described composition independent and excess Gibbs energies are functions of component fraction instead of pair fractions. Application of the presented equation system to the Cu-Fe-S liquid solution gives pair fraction distributions which offer a roughly structural insight into a complex solution and an interpretation of the predictive ability of the one-sublattice modified quasichemical model recently reported for the Cu-Fe-S liquid-solid phase range. Usage of the derived equation (system) is independent of Standard Gibbs energy data and Gibbs energy minimizing software.
Acknowledgments This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Funding Open access funding provided by Montanuniversität Leoben.
Data Availability The data required to reproduce these findings are available within the article or cited in the reference list.

Conflict of interest
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.