Phosphorus Substitution Preference in Ye’elimite: Experiments and Density Functional Theory Simulations

Density functional theory (DFT) simulation has been recently introduced to understand the doping behavior of impurities in clinker phases. P-doped ye’elimite, a typical doping clinker phase, tends to form when phosphogypsum is used to manufacture calcium sulfoaluminate cement (CSA) clinkers. However, the substitution mechanism of P has not been uncovered yet. In this study, the influence of different doping amounts of P on the crystalline and electronic structure of ye’elimite was investigated using backscattered scanning electron microscopy–energy X-ray dispersive spectroscopy, X-ray diffraction tests, Rietveld quantitative phase analysis, and DFT simulations. Furthermore, the substitution preference of P in ye’elimite was revealed. Our results showed that increasing the doping amount of P increased the impurity contents in CSA clinkers, transforming the ye’elimite crystal system from the orthorhombic to the cubic system and decreasing the interplanar spacing of ye’elimite. Based on the calculation results of the defect formation energies, additional energies were required for P atoms to substitute Ca/Al atoms compared with those required for P atoms to substitute S atoms in both orthorhombic and cubic systems of ye’elimite. Combined calculation results of the bond length–bond order and partial density of states showed that the doped P atoms preferably substituted S atoms; the second possible substituted atoms were Al atoms, while there was only a slight possibility for substitution of Ca atoms. The substitution of P atoms for S atoms can be verified based on the elemental distribution in P-doped ye’elimite and the increasing residual CaSO4 contents. The transition of the crystal system and a decrease in the interplanar spacing for ye’elimite can also prove that the substitution of P atoms for Al atoms occurred substantially.


Introduction
Ye'elimite, also known as tetracalcium trialuminate sulfate, is the dominant mineral in calcium sulfoaluminate (CSA) cement and exhibits high tailored dimensional expansion and early strength [1][2][3]. Since the 1970s, CSA cement has been widely used as a rapid repairing material and has recently been promoted for use in numerous domains, including soft soil stabilization and three-dimensional printing [4][5][6]. Moreover, CSA cement has been considered a promising alternative binder owing to its lower CO 2 emissions and energy consumption compared with Portland cement, in which ye'elimite plays an essential role [7][8][9]. The forming conditions of ye'elimite determine the sintering temperatures and prolong time of CSA clinkers, while the main hydration products of CSA cement are produced by reactions between ye'elimite and gypsum [10,11].
The chemical formula of stoichiometric ye'elimite is Ca 4 Al 6 SO 16 which has been identified as an orthorhombic system with a Pcc2 space group at room temperature using diffraction tests and Rietveld method [12]. However, since ye'elimite is a sodalite with a general composition of M 4 (T 6 O 12 )X, certain types of dopants can enter its lattice during the practical manufacturing process of CSA clinkers, forming solid solutions of P-doped ye'elimite clinkers were polished and sprayed with gold for BSEM-EDS analysis. A scanning electron microscope (FEI Quanta 450, FEI, Hillsboro, OR, USA) was used with an acceleration voltage of 20 kV. All synthesized clinkers were ground to pass via 80-µm sieves for XRD tests using an X-ray diffractometer (Bruker D8 advance Davinci design) with Cu Kα 1,2 radiation (λ 1 = 1.5406 Å and λ 2 = 1.5444 Å). Diffraction data were collected at 5 • -80 • (2θ) in the step-scan mode with an operating power of 1600 W (operating voltage was 40 kV, and current was 40 mA). The step size was 0.02 • , and the step duration was set to 0.5 s. To determine the phases in tested clinkers using the crystallography open database (COD) 2013, the Bruker evaluation software was used to analyze the obtained XRD patterns. Thereafter, the Rietveld method was conducted using TOPAS4.2, where the parameters were set according to our previous studies [18][19][20]. Table 2 presents information on all the phases of Rietveld method.

DFT Simulations
The Cambridge sequential total energy package (CASTEP) [41] was used to optimize the geometry of the model and calculate the total energy (E t ), BL-BO distributions and PDOSs of P-doped ye'elimite based on DFT simulations. The generalized functional gradient approximation was combined with the Perdew-Burke-Ernzerhof (PBE) method to calculate the electron exchange correlation [42]. The projection-augmented wave was used for all simulations, and the kinetic energy cutoff of the plane wave base was set to 720 eV [43] (the results of the kinetic energy cutoff test and K-points can be seen in Figure A1). Considering the calculation accuracy, the convergence energy, force, stress and displacement were determined using 10 −6 eV/atom, 0.01 eV/Å, 0.02 GPa and 10 −4 Å, respectively. The lattice relaxation and PDOS calculations were performed on k-point grids with densities of 0.04 and 0.02 (1/Å), respectively. The detailed grids of the two ye'elimite phases are presented in Table A1.
Since dopants can cause transformation of the crystal system of ye'elimite, both orthorhombic and cubic systems were considered and investigated for comparison in this study. The orthorhombic system of ye'elimite exhibits the Pcc2 space group and contains 64 O, 4 S, 24 Al (according to the symmetry, 12 Al in special positions and labeled Al1, while the remaining 12 Al in general positions and labeled Al2 [44]), and 16 Ca atoms in the unit lattice [12]. In the case of the cubic system of ye'elimite, the space group is I43m and the unit lattice comprises 32 O, 2 S, 12 Al, and 8 Ca atoms [16]. Figure 1 shows the crystalline structures of orthorhombic and cubic systems of ye'elimite, and Table A2 presents the lattice parameters.
constructed for the orthorhombic and cubic systems, respectively. It can be noticed that, as all cell/supercell parameters of a, b and c were larger than 9Å (see Table A3), the established cell/supercell totally fulfilled the periodic boundary condition. The crystal structures of pure orthorhombic and cubic ye'elimite used in this paper were obtained through optimizing the crystal structures in Ref. [12] and Ref. [16], respectively. To establish a Pdoped model with a cubic system, one P atom was introduced in the crystalline structure of cubic ye'elimite (labeled as P@SS) at the Ca, Al, and S positions (labeled as P@Ca, P@S, and P@Al, respectively). Furthermore, two P atoms were introduced in the crystalline structure of orthorhombic ye'elimite (labeled as P@ST) at the Ca, S, Al1, and Al2 positions (labeled as P@Ca, P@S, P@Al1, and P@Al2, respectively). All established models for Pdoped ye'elimite can be seen in Sections 3.3 and 4. In this study, the possibility of P doping behavior in ye'elimite was evaluated using Ef, which can be calculated using Equation (1) [48][49][50][51].
where E0 is the bonding energy of pure ye'elimite, and EP is the bonding energy of Pdoped ye'elimite (bonding energy refers to the energy required to separate each pseudopotential atom to infinity, and details of the relationship between bonding energy and total energy can be seen in Table A4). n is the number of P atoms introduced in ye'elimite, and q denotes the net electron number of P-doped ye'elimite (as the raw material used for impurity was Ca3(PO4)2 and the conditions of synthesis were high temperatures, P 5+ preferred to be the most stable charge state of P in this paper). EF represents the Fermi energy relative to EVBM, where EVBM is the energy of the maximum valence band of pure ye'elimite. Note that the doping mechanism can be classified into two types: interstitial and substitutional. Previous studies have shown that the substitutional type can occur more readily than the interstitial type [45][46][47]. Thus, only P doping substitution has been considered in this study. To ensure an equal number of atoms and the same substitutional ratios are involved in calculations, the cells/supercells with sizes of 1 × 1 × 1 and 2 × 1 × 1 were constructed for the orthorhombic and cubic systems, respectively. It can be noticed that, as all cell/supercell parameters of a, b and c were larger than 9Å (see Table A3), the established cell/supercell totally fulfilled the periodic boundary condition. The crystal structures of pure orthorhombic and cubic ye'elimite used in this paper were obtained through optimizing the crystal structures in Ref. [12] and Ref. [16], respectively. To establish a P-doped model with a cubic system, one P atom was introduced in the crystalline structure of cubic ye'elimite (labeled as P@SS) at the Ca, Al, and S positions (labeled as P@Ca, P@S, and P@Al, respectively). Furthermore, two P atoms were introduced in the crystalline structure of orthorhombic ye'elimite (labeled as P@ST) at the Ca, S, Al1, and Al2 positions (labeled as P@Ca, P@S, P@Al1, and P@Al2, respectively). All established models for P-doped ye'elimite can be seen in Sections 3.3 and 4.
In this study, the possibility of P doping behavior in ye'elimite was evaluated using E f , which can be calculated using Equation (1) [48][49][50][51].
where E 0 is the bonding energy of pure ye'elimite, and E P is the bonding energy of P-doped ye'elimite (bonding energy refers to the energy required to separate each pseudopotential atom to infinity, and details of the relationship between bonding energy and total energy can be seen in Table A4). n is the number of P atoms introduced in ye'elimite, and q denotes the net electron number of P-doped ye'elimite (as the raw material used for impurity was Ca 3 (PO 4 ) 2 and the conditions of synthesis were high temperatures, P 5+ preferred to be the most stable charge state of P in this paper). E F represents the Fermi energy relative to E VBM , where E VBM is the energy of the maximum valence band of pure ye'elimite. Furthermore, since the doping rates of P in the orthorhombic and cubic systems of ye'elimite were kept constant in all calculations, normalization was not needed in this study. µ P represents the chemical potentials of impurities (P atoms), and µ 0 represents the chemical potentials of substituted atoms (Ca, Al or S atoms). The terms nµ 0 -nµ P in Equation (1) for different configurations were calculated according to Equations (2)-(4) as below, and details on Gibbs free energies involved in calculations for µ P and µ 0 are presented in Table A2.
As for Equations (2)-(4), µ O was set as the Gibbs free energy of oxygen in the stable pure bulk single-crystal phase, which can be described by Equation (5).

Elemental Distribution in P-Doped Ye'elimite
BSEM-EDS was performed to determine the elemental distribution in P-doped ye'elimite, and the results of P0.10 are shown in Figure 2. Based on the contrasts observed in Figure 2a and the elemental distribution in Figure 2b-e, P-doped ye'elimite accounted for the vast majority of the observed region. In particular, Ca, Al, S, and P were almost evenly distributed in P-doped ye'elimite, indicating that P entered the lattices of ye'elimite. Further, the densities of S in certain areas were slightly attenuated (such as in the area indicated by the dotted line in Figure 2d), while the densities of P slightly increased in the corresponding areas (such as the area indicated by the dotted line in Figure 2e). This phenomenon suggests that P and S exhibit a potential substitutional relationship in P-doped ye'elimite.
Furthermore, since the doping rates of P in the orthorhombic and cubic systems of ye'elimite were kept constant in all calculations, normalization was not needed in this study. μP represents the chemical potentials of impurities (P atoms), and μ0 represents the chemical potentials of substituted atoms (Ca, Al or S atoms). The terms nμ0-nμP in Equation (1) for different configurations were calculated according to Equation (2)-Equation (4) as below, and details on Gibbs free energies involved in calculations for μP and μ0 are presented in Table A2.
As for Equation (2)-Equation (4), μO was set as the Gibbs free energy of oxygen in the stable pure bulk single-crystal phase, which can be described by Equation (5).

Elemental Distribution in P-Doped Ye'elimite
BSEM-EDS was performed to determine the elemental distribution in P-doped ye'elimite, and the results of P0.10 are shown in Figure 2. Based on the contrasts observed in Figure 2a and the elemental distribution in Figure 2b-e, P-doped ye'elimite accounted for the vast majority of the observed region. In particular, Ca, Al, S, and P were almost evenly distributed in P-doped ye'elimite, indicating that P entered the lattices of ye'elimite. Further, the densities of S in certain areas were slightly attenuated (such as in the area indicated by the dotted line in Figure 2d), while the densities of P slightly increased in the corresponding areas (such as the area indicated by the dotted line in Figure 2e). This phenomenon suggests that P and S exhibit a potential substitutional relationship in P-doped ye'elimite.     Figure 3 shows a comparison of the XRD patterns of orthorhombic and cubic systems of ye'elimite to determine differences in their crystalline structures. Except for some low-intensity peaks of tricalcium aluminate (Ca 3 Al 2 O 4 ) and calcium oxide (CaO), all peaks belonged to the orthorhombic and cubic systems of ye'elimite. Furthermore, the orthorhombic and cubic system patterns were generally similar as the major peaks of both the patterns overlapped at a high intensity (such as the peaks at 23.7 • , 33.7 • , and 41.6 • ). However, compared with the cubic system pattern, the orthorhombic system pattern exhibited some additional characteristic peaks (such as those appearing at approximately 18.1 • , 20.6 • , 35.8 • , and 37.3 • ), which can be observed in the magnified image of Figure 1. Thus, the additional characteristic peaks can be used to distinguish the orthorhombic and cubic systems.

Rietveld Quantitative Phases Analysis of Synthesized Clinkers
intensity peaks of tricalcium aluminate (Ca3Al2O4) and calcium oxide (CaO), all longed to the orthorhombic and cubic systems of ye'elimite. Furthermore, the ort bic and cubic system patterns were generally similar as the major peaks of both terns overlapped at a high intensity (such as the peaks at 23.7°, 33.7°, and 41.6°). H compared with the cubic system pattern, the orthorhombic system pattern exhibi additional characteristic peaks (such as those appearing at approximately 18. 35.8°, and 37.3°), which can be observed in the magnified image of Figure 1. T additional characteristic peaks can be used to distinguish the orthorhombic and c tems.    Figure 4 shows the XRD patterns of P-doped ye'elimite in comparison with the orthorhombic system. All diffraction peaks were assigned according to the standard diffraction patterns of ye'elimite (COD#4001772 and COD#4511960), anhydrite (COD#5000040), mayenite (COD#2102955), calcium phosphate (COD#1517238), aluminum phosphate (COD#9006404), monocalcium aluminate (COD#4308075), CaO (COD#7200686), and tricalcium aluminate (COD#9015966).
As shown in Figure 4a, more peaks for impurities, such as Ca 3 (PO 4 ) 2 , AlPO 4 , CaAl 2 O 4 , and CaSO 4 , can be clearly detected in samples of P0.15 and P0.20, compared with ST and P0.10 samples. The peak intensities for Ca 3 (PO 4 ) 2 , AlPO 4 , and CaSO 4 in particular gradually increased as the P doping amount increased. Moreover, according to Figure 4b, higher phosphorus doping amounts yielded lower intensities of characteristic peaks for orthorhombic ye'elimite, indicating that ye'elimite gradually transformed from the orthorhombic to the cubic system. Furthermore, based on the comparison of three major peaks of ye'elimite in Figure 4c, the 2θ values for the corresponding peaks increased as the P doping amounts increased. When combined with the equation of Bragg diffraction (Equation (6)), it can be deduced that higher P doping amounts decrease the ye'elimite interplanar spacing. 2dsinθ = nλ (6) Rietveld quantitative phase analysis (RQPA) was performed to further determine the compositions of clinkers containing P-doped ye'elimite. Figure 5 shows a selected range (15-45 • /2θ) of the Rietveld plots for P0.15, and the results are presented in Table 3. According to the RQPA results, as the P doping amounts increased from 0 mol to 0.20 mol, the total ye'elimite content decreased from 96.4 to 84.6 wt.%, while the impurity contents of calcium aluminate phases and calcium/aluminum phosphate increased from 2.8 and 0 wt.% to 11.6 and 2.7 wt.%, respectively. Note that the CaSO 4 residues accounted for 0.6 and 1.2 wt.% in P0.15 and P0.20 clinkers, respectively. Additionally, based on the RQPA results, the transformation of the crystal system caused by large P doping amounts can be reflected by the ratio of cubic ye'elimite to orthorhombic ye'elimite. Figure 4 shows the XRD patterns of P-doped ye'elimite in comparison with the orthorhombic system. All diffraction peaks were assigned according to the standard diffraction patterns of ye'elimite (COD#4001772 and COD#4511960), anhydrite (COD#5000040), mayenite (COD#2102955), calcium phosphate (COD#1517238), aluminum phosphate (COD#9006404), monocalcium aluminate (COD#4308075), CaO (COD#7200686), and tricalcium aluminate (COD#9015966).   As shown in Figure 4a, more peaks for impurities, such as Ca3(PO4)2 and CaSO4, can be clearly detected in samples of P0.15 and P0.20, comp P0.10 samples. The peak intensities for Ca3(PO4)2, AlPO4, and CaSO4 in ally increased as the P doping amount increased. Moreover, according to phosphorus doping amounts yielded lower intensities of characteristic rhombic ye'elimite, indicating that ye'elimite gradually transformed from bic to the cubic system. Furthermore, based on the comparison of thre ye'elimite in Figure 4c, the 2θ values for the corresponding peaks increa ing amounts increased. When combined with the equation of Bragg diff (6)), it can be deduced that higher P doping amounts decrease the ye'el spacing.

2dsinθ = nλ
Rietveld quantitative phase analysis (RQPA) was performed to furt compositions of clinkers containing P-doped ye'elimite. Figure 5 shows (15-45°/2θ) of the Rietveld plots for P0.15, and the results are presente cording to the RQPA results, as the P doping amounts increased from 0 the total ye'elimite content decreased from 96.4 to 84.6 wt.%, while the i of calcium aluminate phases and calcium/aluminum phosphate increas wt.% to 11.6 and 2.7 wt.%, respectively. Note that the CaSO4 residues and 1.2 wt.% in P0.15 and P0.20 clinkers, respectively. Additionally, ba results, the transformation of the crystal system caused by large P dopin reflected by the ratio of cubic ye'elimite to orthorhombic ye'elimite.  15. Observed data resulted from XRD test, calculated line resulted from Rietveld refinement, and the differential line resulted from the difference between observed data and calculated line.

Defect Formation Energies of P-Doped Ye'elimite
The E f of P-doped ye'elimite with orthorhombic and cubic systems was calculated by considering different configurations to determine the substitution preference of P doping in ye'elimite. The doped models used for calculations are shown in Figures 6 and 7, and the lattice parameters of the doped models are presented in Table A3. The calculations E f results are shown in Figure 8 (refer to Table A4 for details). The relative stabilities of defects and possibilities of reaction can be reflected by the values of formation energies (E f ), in which lower E f represents more stable configuration and the easier reaction [45,47]. Thus, the E f values can be used to estimate the possibility and stability of doping behavior. According to the results presented in Figure 8, P atoms most preferentially tended to substitute S atoms since the E f values of P@S configurations remained to be minimum (1.27 and 1.20 eV for orthorhombic and cubic systems, respectively). The elemental distribution of P and S obtained using BSEM-EDS can be used to identify the substitution of P atoms for S atoms. Additionally, the increasing CaSO 4 residue in the clinkers with P-doped ye'elimite can confirm the substitutional relationship between P and S since CaSO 4 generally cannot be detected in clinkers of pure ye'elimite owing to the inevitable decomposition of CaSO 4 during the temperature-rise period [52,53].
XRD test, calculated line resulted from Rietveld refinement, and the differential line resulted fro the difference between observed data and calculated line.

Defect Formation Energies of P-Doped Ye'elimite
The Ef of P-doped ye'elimite with orthorhombic and cubic systems was calculated b considering different configurations to determine the substitution preference of P dopin in ye'elimite. The doped models used for calculations are shown in Figures 6 and 7, an the lattice parameters of the doped models are presented in Table A3. The calculations results are shown in Figure 8 (refer to Table A4 for details). The relative stabilities of d fects and possibilities of reaction can be reflected by the values of formation energies (E in which lower Ef represents more stable configuration and the easier reaction [45,47 Thus, the Ef values can be used to estimate the possibility and stability of doping behavio According to the results presented in Figure 8, P atoms most preferentially tended to su stitute S atoms since the Ef values of P@S configurations remained to be minimum (1.2 and 1.20 eV for orthorhombic and cubic systems, respectively). The elemental distributio of P and S obtained using BSEM-EDS can be used to identify the substitution of P atom for S atoms. Additionally, the increasing CaSO4 residue in the clinkers with P-dope ye'elimite can confirm the substitutional relationship between P and S since CaSO4 gene ally cannot be detected in clinkers of pure ye'elimite owing to the inevitable decompos tion of CaSO4 during the temperature-rise period [52,53].  Alternatively, when P atoms substituted Ca/Al atoms, the Ef values remained positive and followed the alignment of P@Al (P@SS) < P@Al1 (P@ST) ≈ P@Al2 (P@ST) < P@Ca (P@ST) < P@Ca (P@SS). The alignment indicates that further additional energies were required for P atoms to substitute Ca atoms than for P atoms to substitute Al atoms. However, since P-doped ye'elimite was synthesized at a high temperature (1300 °C in this study), the required additional energies could be easily obtained to overcome the barrier for P atoms to substitute Ca/Al atoms.

Electronic Structural Matching
To further estimate the possibility of different P-doped configurations, the BO-BL distributions in P-doped ye'elimite were calculated according to the Mulliken population rule [54]. Generally, covalent bonds exhibit BO values close to 1, whereas ionic bonds show BO values close to 0. Further, greater overlaps in BO-BL distributions between substituted atoms and dopants indicate a greater possibility for doping behavior to form solid  Alternatively, when P atoms substituted Ca/Al atoms, the Ef values r tive and followed the alignment of P@Al (P@SS) < P@Al1 (P@ST) ≈ P@Al2 (P@ST) < P@Ca (P@SS). The alignment indicates that further additional en quired for P atoms to substitute Ca atoms than for P atoms to substitute A ever, since P-doped ye'elimite was synthesized at a high temperature (1 study), the required additional energies could be easily obtained to overco for P atoms to substitute Ca/Al atoms.

Electronic Structural Matching
To further estimate the possibility of different P-doped configuratio distributions in P-doped ye'elimite were calculated according to the Mullik rule [54]. Generally, covalent bonds exhibit BO values close to 1, where show BO values close to 0. Further, greater overlaps in BO-BL distribution stituted atoms and dopants indicate a greater possibility for doping behavio solutions. Figure 9 shows the BL-BO distributions in pure and P-doped ye'elim ent configurations. For both orthorhombic and cubic ye'elimite, S-O bonds Alternatively, when P atoms substituted Ca/Al atoms, the E f values remained positive and followed the alignment of P@Al (P@SS) < P@Al1 (P@ST) ≈ P@Al2 (P@ST) < P@Ca (P@ST) < P@Ca (P@SS). The alignment indicates that further additional energies were required for P atoms to substitute Ca atoms than for P atoms to substitute Al atoms. However, since P-doped ye'elimite was synthesized at a high temperature (1300 • C in this study), the required additional energies could be easily obtained to overcome the barrier for P atoms to substitute Ca/Al atoms.

Electronic Structural Matching
To further estimate the possibility of different P-doped configurations, the BO-BL distributions in P-doped ye'elimite were calculated according to the Mulliken population rule [54]. Generally, covalent bonds exhibit BO values close to 1, whereas ionic bonds show BO values close to 0. Further, greater overlaps in BO-BL distributions between substituted atoms and dopants indicate a greater possibility for doping behavior to form solid solutions. Figure 9 shows the BL-BO distributions in pure and P-doped ye'elimite with different configurations. For both orthorhombic and cubic ye'elimite, S-O bonds were observed with relatively high BO bonds, Ca-O bonds showed BO values of close to 0, and Al-O bonds were located in the regions between S-O and Ca-O bonds. Furthermore, when P atoms were introduced to substitute S atoms in P-doped ye'elimite, P-O bond regions overlapped well with S-O bond regions. In the configurations where P atoms substituted Al atoms, partial P-O bonds approached closely to the regions of Al-O bonds in both orthorhombic and cubic systems (refer to the marked regions in Figure 9a,b). However, when P atoms substituted Ca atoms, most of the P-O bonds remained away from the regions of Ca-O bonds.
Al atoms, partial P-O bonds approached closely to the regions of Al-O bonds in both orthorhombic and cubic systems (refer to the marked regions in Figure 9a,b). However, when P atoms substituted Ca atoms, most of the P-O bonds remained away from the regions of Ca-O bonds.
The electronic structures of bonds formed between O and substituted/doped atoms can be further investigated using PDOS. Figure 10 shows the PDOS calculation results for various elements in pure and P-doped ye'elimite. According to Figure 10a,c, for pure ye'elimite with both orthorhombic and cubic systems, the O 2p orbitals hybridized with Al 3p and S 3p orbitals to form Al-O and S-O bonds, respectively, while the O 2p orbitals hybridized with Ca 3d orbitals to form Ca-O bonds. To form P-O bonds in P-doped ye'elimite ( Figure 10b,d), the O 2p orbitals hybridized with P 3s and P 3p orbitals. Compared with pure ye'elimite, the PDOS of P@S (in both P@ST and P@SS) exhibited the most similar distribution for all P-doped ye'elimite configurations, whereas partial overlapping (between −4 and −2 eV) could be observed for the PDOS of P@Al, including P@Al in P@SS and P@Al1/P@Al2 in P@ST. Additionally, similar electronic contributions between P-O and Ca-O bonds are seldomly observed in both orthorhombic and cubic structures of ye'elimite.
Combining the BL-BO and PDOS calculation results, for P-doped ye'elimite with both orthorhombic and cubic systems, the P atoms preferably substituted S atoms; the second possible substituted atoms were Al, while the substitution for Ca atoms exhibited only a slight possibility.  The electronic structures of bonds formed between O and substituted/doped atoms can be further investigated using PDOS. Figure 10 shows the PDOS calculation results for various elements in pure and P-doped ye'elimite. According to Figure 10a Combining the BL-BO and PDOS calculation results, for P-doped ye'elimite with both orthorhombic and cubic systems, the P atoms preferably substituted S atoms; the second possible substituted atoms were Al, while the substitution for Ca atoms exhibited only a slight possibility.

Discussion
According to the Ef, BL-BO, and PDOS calculation results, P atoms were more likely to substitute S atoms, which is a spontaneous process, rather than Al/Ca atoms. In addition, the substitution of P atoms for S atoms can be verified using the elemental distribution in P-doped ye'elimite ( Figure 2) and the increasing residual CaSO4 contents (Table 1).
According to the XRD patterns demonstrated in Figure 4 and the RQPA results presented in Table 1, increasing the P doping amount caused the transformation of the crystal system from the orthorhombic to the cubic system. Previous studies have shown that the freezing of anionic groups caused by dopants prompted the transformation of the crystal system of ye'elimite during cooling [16]. The freezing effect can also be reflected by the E values of ye'elimite, where a lower E value generally represents a more stable structure

Discussion
According to the E f , BL-BO, and PDOS calculation results, P atoms were more likely to substitute S atoms, which is a spontaneous process, rather than Al/Ca atoms. In addition, the substitution of P atoms for S atoms can be verified using the elemental distribution in P-doped ye'elimite ( Figure 2) and the increasing residual CaSO 4 contents (Table 1).
According to the XRD patterns demonstrated in Figure 4 and the RQPA results presented in Table 1, increasing the P doping amount caused the transformation of the crystal system from the orthorhombic to the cubic system. Previous studies have shown that the freezing of anionic groups caused by dopants prompted the transformation of the crystal system of ye'elimite during cooling [16]. The freezing effect can also be reflected by the E values of ye'elimite, where a lower E value generally represents a more stable structure that tends to be retained during the crystalline transformation process [12]. To confirm the stability of the P-doped ye'elimite structure, the E values for P substituting different atoms in orthorhombic and cubic structures of ye'elimite are compared in Figure 11. terials 2021, 14, x FOR PEER REVIEW that tends to be retained during the crystalline transformation process stability of the P-doped ye'elimite structure, the E values for P substitu in orthorhombic and cubic structures of ye'elimite are compared in Fi According to the results depicted in Figure 11, for pure ye'elimit system demonstrated more stability at the room temperature than the agrees with the experimental results presented in Figure 3. Similarly with P atoms substituting Ca/S atoms, the orthorhombic system of P-d hibited lower E values than the cubic system of P-doped ye'elimite atoms were introduced to substitute Al atoms, the E value for the cub than that for the orthorhombic system. When all configurations with P ent atoms in P@ST and P@SS were combined, only the substitution of P would result in the transformation of the crystal system of P-doped y Additionally, based on the XRD patterns shown in Figure 4, incr amount resulted in a gradual decrease in the ye'elimite interplanar sp erally explained by smaller ionic radii of dopants compared with [14,18,20]. According to the radii of the ions involved in P-doped ye' the radii of phosphonium ions with different coordination numbers of sulphion, implying that the decrease in the ye'elimite interplana attributed to the substitution of P atoms for S atoms. Considering the substitution for Ca atoms revealed by BL-BO and PDOS results, the oms for Al atoms was the most likely driving force for the decreas planar spacing. Table 4. Radii of ions involved in P-doped ye'elimite [55]. According to the results depicted in Figure 11, for pure ye'elimite, the orthorhombic system demonstrated more stability at the room temperature than the cubic system, which agrees with the experimental results presented in Figure 3. Similarly, for configurations with P atoms substituting Ca/S atoms, the orthorhombic system of P-doped ye'elimite exhibited lower E values than the cubic system of P-doped ye'elimite. However, when P atoms were introduced to substitute Al atoms, the E value for the cubic system was lower than that for the orthorhombic system. When all configurations with P substituting different atoms in P@ST and P@SS were combined, only the substitution of P atoms for Al atoms would result in the transformation of the crystal system of P-doped ye'elimite.
Additionally, based on the XRD patterns shown in Figure 4, increasing the P doping amount resulted in a gradual decrease in the ye'elimite interplanar spacing, which is generally explained by smaller ionic radii of dopants compared with substituted atoms [14,18,20]. According to the radii of the ions involved in P-doped ye'elimite (Table 4), all the radii of phosphonium ions with different coordination numbers are larger than that of sulphion, implying that the decrease in the ye'elimite interplanar spacing cannot be attributed to the substitution of P atoms for S atoms. Considering the slight possibility of substitution for Ca atoms revealed by BL-BO and PDOS results, the substitution of P atoms for Al atoms was the most likely driving force for the decreasing ye'elimite interplanar spacing. Combining the above analysis with the consideration of charge equilibrium in the lattices, the co-substituting P-doped ye'elimite were proposed to be the most probably models (see Figure 12), in which P atoms substituted two S atoms and one Al atom in cell/supercell of orthorhombic and cubic ye'elimite, respectively (P@S&P@Al(ST) refers to the co-substituting P-doped ye'elimite with orthorhombic system, P@S&P@Al(SS) refers to the co-substituting P-doped ye'elimite with cubic system). It should be noted that the E f of the co-substituting P-doped ye'elimite was, respectively, equaled to −0.22 and −0.18 eV for orthorhombic and cubic system (the comparison of E f in this paper and references can be seen in Table A5), which illustrates the co-substituting models were more rational than the above single substituting models. aterials 2021, 14, x FOR PEER REVIEW Combining the above analysis with the consideration of charg lattices, the co-substituting P-doped ye'elimite were proposed to b models (see Figure 12), in which P atoms substituted two S atoms cell/supercell of orthorhombic and cubic ye'elimite, respectively (P@ to the co-substituting P-doped ye'elimite with orthorhombic system fers to the co-substituting P-doped ye'elimite with cubic system). It s the Ef of the co-substituting P-doped ye'elimite was, respectively, e −0.18 eV for orthorhombic and cubic system (the comparison of Ef in ences can be seen in Table A5), which illustrates the co-substituting rational than the above single substituting models. To further examine the rationality of the co-substituting P-dope mental and simulated XRD patterns were compare d in Figure 13. P@ to the mixture of P@S＆P@Al(ST) and P@S＆P@Al(SS), and the ratio 0.89 which is the same as P0.20. It can be noticed that major peaks experiment and simulation were almost coincident. Additionally, the BO and PDOS were also calculated and shown in Figures A2 and A3, w the rationality of the model of co-substituting P-doped ye'elimite fro tronic structure. To further examine the rationality of the co-substituting P-doped ye'elimite, experimental and simulated XRD patterns were compare d in Figure 13. P@S&P@Al(Mix) refers to the mixture of P@S&P@Al(ST) and P@S&P@Al(SS), and the ratio of SS/ST was set as 0.89 which is the same as P0.20. It can be noticed that major peaks of ye'elmite for the experiment and simulation were almost coincident. Additionally, the distributions of BL-BO and PDOS were also calculated and shown in Figures A2 and A3, which can also verify the rationality of the model of co-substituting P-doped ye'elimite from the aspect of electronic structure. 0.89 which is the same as P0.20. It can be noticed that major peaks experiment and simulation were almost coincident. Additionally, the BO and PDOS were also calculated and shown in Figures A2 and A3, w the rationality of the model of co-substituting P-doped ye'elimite from tronic structure.

Conclusions
The effect of different doping amounts of P on the crystalline structures of ye'elimite was investigated, and the doping mechanism was revealed using BSEM-EDS, XRD tests, the Rietveld method, and DFT simulations. The following conclusions were obtained:

1.
Experiments and Rietveld analysis confirmed that doped P entered ye'elimite to form P-doped solid solutions, resulting in increased impurity contents in clinkers, a crystal system transformation from the orthorhombic to the cubic system, and a decrease in the ye'elimite interplanar spacing.

2.
Based on calculation results of E f , additional energies were required for P atoms to substitute Ca/Al atoms compared with those for substituting S atoms for both orthorhombic and cubic structures of ye'elimite. By combining the BL-BO and PDOS calculation results, the doped P atoms preferably substituted S atoms; the second possible substituted atoms were Al atoms, while there was only a slight possibility for substitution of Ca atoms. 3.
The substitution of P atoms for S atoms can be verified using the elemental distribution in P-doped ye'elimite and the increasing residual CaSO 4 contents. The crystal system transformation and a decrease in the ye'elimite interplanar spacing can also imply that the substitution of P atoms for Al atoms occurred substantially.

4.
Based on analysis of phosphorus substitution preference in ye'elimite, the co-substituting P-doped ye'elimite were proposed to be the most probably models. Through values of E f , comparison of XRD patterns and electronic structural matching, the rationality of the model for co-substituting P-doped ye'elimite can be verified. Appendix A data curation, investigation, visualization; C.Y., data curation, investigation; C.C., proje istration; J.C., data curation, supervision. All authors have read and agreed to the publishe of the manuscript. Appendix A Figure A1. Results of kinetic energy cutoff test and K-points (SS). Structure of SS was selected to conduct the kinetic energy cutoff test and K-points test, as the structure contains all pseudopotentials. Table A1. K-point grids for simulation models.

Supercells Relax PDOS
Orthorhombic ye'elimite  1 Lattice parameters of pure orthorhombic and cubic ye'elimite used in this paper were obtained through optimizing the crystal structures in Ref. [12] and Ref. [16], respectively.  1 The formula for calculating the binding energy is: where E b is the bonding energy, E t is the total energy of unit cell, and E Al , E O , E Ca , E S and E P are the energies of pseudo atomic Al, O, Ca, S and P respectively.        Table A5.
Comparison of E f in this paper and references.

Crystal System Configurations
Value of E f /eV Sources