Calculation of Phase Diagrams and First-Principles Study of Germanium Impacts on Phosphorus Distribution in Czochralski Silicon

Phosphorus (P) is one of the most widely used donor dopants for fabricating a low-resistivity silicon (Si) substrate. However, its volatile nature and the relatively small equilibrium segregation coefficient in Si at the melting temperature of Si impede the efficient and effective growth of low-resistivity Czochralski (CZ) Si single crystal. The primary objective of this work is to theoretically perceive the influence of germanium co-doping on the heavily P-doped Si crystal by means of CALculation of PHase Diagrams (CALPHAD) approaches and density functional theory (DFT) calculations. Phase equilibria at the Si-rich corner of the Si-Ge-P system has been thermodynamically extrapolated based on robust thermodynamic descriptions of involved binary systems, where Si-P and Ge-P have been re-assessed in this work. Phase diagram calculation results indicate that at a given P concentration (e.g. 0.33 at.% P) Ge co-doping lowers the solidification temperature of the Si(Ge, P) alloys, as well as the relevant equilibrium segregation coefficients of P in the doped Si. DFT calculations simulated the formation of (i) monovacancy in Si as well as (ii) solutions of Si(P) and Si(Ge) with one dopant substitutionally inserted in 64- and 216-atom Si cubic supercells. Binding energies were calculated and compared for Ge-Ge, Ge-P and P-P bonds positioning at the first nearest-neighbors (1NN) to the third nearest-neighbors (3NN). P-P bonds have the largest bonding energy from 1NN to 3NN configurations. The climbing image nudged elastic band method (CL-NEB) was utilized to calculate the energy barriers of P 1NN jump in the 64-atom Si cubic supercell with/without a neighboring Ge atom. With Ge present, a higher energy barrier for P 1NN jump was obtained than that without involving Ge. This indicates that Ge can impede the P diffusion in Si matrix.


Introduction
With the ongoing trend towards multifunction and miniaturization in electronic devices, contact resistivity associated power consumption is increasingly pronounced. This, together with the demand for low on-resistance of power devices, drives the need for lowering the resistivity of substrate that can then be utilized to fabricate the epitaxial structure. Phosphorus is a favorable dopant species. Despite the mature technology in Czochralski (CZ) Si growth, the growth of heavily P-doped CZ Si crystal with high yield is not without challenges. 1 Basically, issues related to crystal growth are provoked by the small equilibrium distribution (segregation) coefficient (K eq ), the highly volatile nature of P, and the retrograde character of the Si(P) solidus. [1][2][3][4] The small K eq results in the uneven distribution of P along the Si ingot (i.e. concentration of P increasing along the crystal body), thus affecting the yield of the demanding resistivity level of the CZ Si ingots. This phenomenon is marginal in lightly P-doped CZ Si since one can optimize the fabrication parameters to achieve a good balance between the enrichment of P at crystal/melt interface and the loss of P resulted from evaporation. However, the 1 3 evaporation rate of P increases with increasing doping level, which leads to severe loss of P during melting, homogenization, and crystal pulling stages of CZ Si growth process. In addition, the retrograde P solubility in Si results in the non-constant K eq values of P at light and high doping levels. To precisely control the crystal growth process, the relevant K eq values at high P doping levels is indispensable since the widely accepted K eq value of 0.35 is no longer applicable to the growth of low resistivity P-doped Si crystal. Moreover, P atoms enter substitutionally in Si. [5][6][7] Misfit between the covalent radii of P and Si introduces the solute strain in the crystal thus generating an internal stress. [8][9][10] According to Vegard's law, substitutional P atoms in Si generates a lattice contraction with respect to Si crystal. 9,11,12 These issues give rise to potential challenges in controlling dislocation-free silicon crystal growth and fabricating good quality epitaxial layer (epi-layer). For instance, when a lightly doped (n −) epi-layer is grown on a heavily doped (n +) Si substrate, the misfit dislocation can form at the substrate/epi-layer interface. 10 Furthermore, the composition gradient drives the diffusion of P from heavily doped substrate towards the lightly doped epi-layer during subsequent processing, thus likely forming another auto-doped layer between the n + substrate and the n − epi-layer.
Co-doping (i.e. using a third element in addition to the selected dopant) has been proposed as a solution to the aforementioned issues, such as to lower the evaporation of the volatile species, to enlarge the K eq value of selected dopant, to compensate for the misfit due to atomic size differences and to retard the diffusion of dopant. 10,13 Germanium (Ge) is a promising dopant in Si since it has the potential to control the defect formation. 10,14 In this respect, Ge can be an electrically neutral co-dopant in P-doped Si. To explore the Ge influence on the P-doped Si crystal, investigations were performed on the carrier mobility (µ), resistivity (ρ), K eq and process parameters of Ge-P co-doped CZ Si growth. 10,15,16 The Hall-measurements on the P-doped Si x Ge 1−x (0.93 < x < 0.96) single crystal showed that when carrier concentration was higher than 10 18 atoms/cm 3 , the obtained µ and resistivity values were comparable to those in the P-doped Si. 15 The K eq of P in Si x Ge 1−x (0.93 < x < 0.96) single crystal significantly depended on the doping level. 16 Kawazoe et al. 10 proposed the empirical composition relationship between P and Ge in CZ Si single crystal (i) to avoid the abnormal crystal growth and (ii) to prevent the lattice distortion resulted from high doping level of P. 10 These studies suggest that Ge is a plausible co-doping candidate to grow the low-resistivity P-doped CZ Si crystal; however, their results were case dependent. There is no satisfactory understanding of the Ge impact on the P behavior in liquid and solid Si yet. Consequently, the phase equilibria at the Si-rich corner of the Ge-P-Si ternary system is of significant importance.
Here, we implemented calculation of phase diagrams (CALPHAD) and density functional theory (DFT) calculations at the Si-rich corner of the Si-Ge-P ternary system to reveal the influence of Ge co-doping on the heavily P-doped Si. CALPHAD methods in this work serve the purpose of evaluating the solidification path of the Si(P, Ge) ternary alloy and calculating the K eq of P in relevant Si(P, Ge). DFT calculations attempt to discover the effect of Ge on preventing P diffusion in Si.

A Brief Review of Previous Work on Phase Equilibria
To date, neither phase equilibria nor thermodynamic properties of the Ge-P-Si ternary system is available. Nevertheless, phase diagrams of the involved three binary end-member systems, namely, Ge-Si, Ge-P and Si-P, have been experimentally and thermodynamically studied. 2,3,[17][18][19][20][21][22][23][24][25][26] The Ge-Si system has a continuous solid solution phase diagram. The first review and thermodynamic evaluation of its phase equilibria was reported by Olesinski et al. 17 Then in 1992, Bergman et al. 18 determined the partial enthalpy of Si in the Si-Ge liquid state at 1327 K in the concentration range of 0 < x Si < 0.03. Based on the derived excess Gibbs energy of liquid over the entire composition range, they performed the thermodynamic investigation on the Si-Ge phase diagram that agreed well with the assessment work published in Ref. 17. The thermodynamic description of Si-Ge is included in the work of Scientific Group Thermodata Europe (SGTE) where the modification has been built upon the combination of Refs. 17 and 18. In this work, the thermodynamic parameters for the Ge-Si system are directly taken from the SGTE Solution database version 4.9 (SSOL4). 19 Among these three binary systems, the Si-P system has been subjected to the most extensive theoretical investigations. 2,3,[20][21][22][23] The first thermodynamic evaluation of this system was achieved by Olesinski et al. 20 They proposed the equilibrium phase diagram of the Si-P binary system covering the entire composition range. In 2009, Tang et al. developed a thermodynamic database for fabricating solar cell grade silicon materials (SOG-Si). That database covered the reassessment of the Si-rich corner phase equilibria of the Si-P system. 21 Arutyunyan et al. 23 optimized the thermodynamic description of phases at the Si-rich side as well by using their in-house software. Jung and Zhang 2 conducted the first complete thermodynamic assessment of the Si-P system. Later, Liang and Schmid-Fetzer 3 thermodynamically reoptimized the description of the Si-P system through critically reviewing all available original experimental data. In addition to SiP, they introduced a second intermediate compound SiP 2 in the Si-P phase diagram. The challenge of reproducing all the measured solidus and solvus of Si has 1 3 been highlighted. Because of the largely scattered solubility of P in Si solid, they developed two alternative sets of selfconsistent thermodynamic parameters. 3 There is no doubt that these independent sets of thermodynamic descriptions for Si-P can serve their individual targets. However, no clear definition is available for the lattice stability of metastable diamond structural P (P_diamond). It is critical for the description of the Gibbs free energy of the relevant diamond solution phases, such as Si(P) and Ge(P). In this work, DFT calculations were conducted to derive the lattice stability of metastable P_diamond. Accordingly, the modification of Si-P is performed.
The Ge-P binary is not less complex than the Si-P system. Regarding the temperature-composition (T − x) phase diagram of the Ge-P binary system, it was initially constructed based on differential thermal analysis (DTA) results of 12 bulk alloys (i.e. 10 at.%, 20 at.%, 30 at.%, 40 at.%, 45 at.%, 50 at.%, 55 at.%, 60 at.%, 65 at.%, 75 at.%, 80 at.% and 90 at.% P). 24 These alloys were prepared by melting at 960°C for 10-12 h under elevated argon pressure, and subsequently annealed at unspecified temperature for 300 h. DTA were then carried out in a structured chamber into which inert gas (nitrogen and argon) was introduced under a certain pressure (unspecified in Ref. 24). The acquired heating curves were adopted to construct the T − x diagram in the temperature range of 400-1000°C. The T − x diagram includes the liquidus boundaries and two invariant reactions, i.e. (i) peritectic reaction of liquid + Ge → GeP at 725°C (liquid composition ~ 63 at.% P), and (ii) an undefined invariant reaction involving liquid, P and GeP. Although the sample preparation process has been introduced in detail by Goncharov et al., 24 they did not clarify whether the constructed T − x phase diagram of Ge-P was referring to a constant pressure. In their later publication, 25 in addition to the DTA results, they have reported the data measured via a static manometric method. Based on the measured data and thermodynamic analysis of the interaction of components, Ugai et al. 25 constructed the pressure-temperature-composition (P − T − x) diagram of the Ge-P binary system that was projected on the T − x, P − T and P − x planes. The solid solubility of P in Ge_diamond has been determined by means of microhardness measurements 26 and Hall effect measurements. 27 Abrikosov et al. 26 investigated the solubility of P in Ge_diamond at 500°C, 600°C, 700°C, 800°C, 850°C and 900°C. Samples that have been made of single-crystal Ge (impurity level under 10 −5 %) and P_red (99.95%P), were annealed at a given temperature for 320-1000 h to achieve the equilibrium. By combining microstructural and chemical analyses with microhardness measurements, they proposed the solubility boundary of P in Ge_diamond. It was found that the maximum P solubility in Ge_diamond was 0.45 at.% P at 600°C. However, the equilibrium solid solubility of P in Ge at 670°C was 0.12 at.% according to the supersaturated Ge(P) solid solution decomposition investigation. 28 Hall effect measurements were conducted on the P-doped CZ Ge single-crystal to study the solubility of P in Ge_diamond in the temperature range of 300-880°C. 27 Before characterizing the Hall carrier concentration, Fistul et al. have pretreated samples by annealing at 800°C for 2 h and quenching in ethylene glycol (with the cooling rate of 200°C/s). Besides, the measurement error has been controlled to be under 5% by adopting proper shape of sample and high magnetic fields (16 kA/m). Based on the collected kinetic annealing curves, Fistul et al. have found that P in Ge_diamond exhibited retrograde solubility. The maximum concentration of P in Ge_diamond was 7 × 10 19 atom/cm 3 (or, ~ 0.16 at.% P) at 800°C, 27 rather than the previously reported 0.45 at.% at 600°C. 26 Nevertheless, the P solubility in Ge at 650°C was ~ 4.7 × 10 19 atom/cm 3 (or, ~ 0.11 at.% P), which is in good agreement with the experimental data (0.12 at.%) in Ref. 27.
Thermodynamic properties of GeP have been relatively well studied. [29][30][31][32][33][34][35][36][37][38][39] Ugai et al. 31 measured the low temperature heat capacity of GeP by means of a calorimeter in temperature range of 6.99-304.83 K. The high temperature heat capacity of GeP has never been experimentally measured. The C p expression of GeP above room temperature was estimated in literature. [32][33][34] For instance, in the first compilation of thermodynamic properties of inorganic compounds, Barin et al. 32  ). [32][33][34] Zumbusch et al. 29 and Süss et al. 30 have determined the dissociation pressure of GeP in the temperature ranges of 773-832 K and 773-973 K, respectively. However, we noticed an error for the unit of the y-axis in the decomposition pressure graph, i.e. Fig. 3 in Ref. 30. It is in Pa rather than in kPa. One reason is that the experimental data from [2930] was in the level of 10 4 Pa. Those data were superimposed in Fig. 3 of Ref. 30; however, they were presented in the range of 10 4 -10 5 kPa rather than 10-10 2 kPa. The error was also evidenced by proposed expression of the total dissociation pressure of GeP as a function of temperature. 30 Zumbusch et al. 29 and Süss et al. 30 also calculated ΔH 0 298 and S 0 298 of GeP according to the reactions: 4GeP = 4Ge(s) + P 4 (g) at 540°C and Ge(s) + P(s, white) = GeP(s) and, Ge(s) + 1/4P 4 (g) = GeP(s) and P 4 (g) = 2P 2 (g), respectively.
Olesinski et al. 40 summarized and reviewed the phase equilibria [24][25][26][27][28] and thermodynamic properties 29-32 of the Ge-P system that were published before 1985. They assessed the steady-state temperature-composition (T − X) phase diagram of Ge-P at 4600 kPa 40 with respect to the reported data 1 3 in literature. 25,29,30,32 For instance, the liquidus boundary of Ge-P phase diagram at 4600 kPa has been determined based on the determined liquidus data from. 25 However, Olesinski et al. 40 has stated that the evaluated phase diagram was valid only on the condition of that the pressure should guarantee the melting of P rather than sublimation in the given temperature ranges. For instance, the pressure should be no less than 4300 kPa, as proposed in Ref. 24. In addition, the solid solubility of P in Ge_diamond was ignored in Ref. 40. In view of this, the re-assessment of the Ge-P phase diagram is essential to reproduce the Ge(P) solidus.

CALPHAD
In the present work, the phase equilibria of the Ge-P-Si ternary system are thermodynamically extrapolated by combining the Ge-Si database from SGTE, 19 the currently optimized Ge-P and Si-P descriptions. The lattice stability of the element i (i = Ge, P and Si) is attributed to the enthalpy of its stable state at 298.15 K and 100 kPa, H SER i , as recommended by SGTE. 41 The Ge-P-Si ternary system includes three solution phases (i.e. liquid, Si_diamond and Ge_diamond), as well as three intermetallic compounds (namely, GeP, SiP and SiP 2 ). An ordinary substitutional solution model is applied to describe these solution phases. Thermodynamic parameters were evaluated using the PARROT optimization module in the 2018b version of the Thermo-Calc software package. 42 The molar Gibbs energy of a solution phase can be represented as a sum of the Gibbs energy for the involved components i, the ideal entropy term describing a random mixing of components, and the excess Gibbs energy describing the deviation from the ideal behavior, i.e., where R is the gas constant, T is the absolute temperature, x i is the molar fraction of component i (i = Ge, P, Si), o G denotes the molar Gibbs energy of i, with reference state of . The term ex G describes the excess Gibbs energy that has not been represented by the first two terms. For the involved binary systems (i.e. Ge-P, Ge-Si, and P-Si), the molar excess Gibbs energies, Bin.ex G , were expressed as: For solution phases in binary systems, the Redlich-Kister (RK) power series has been adopted to describe the composition dependence for the binary interaction parameters L v ij : The values of L v Ge,Si were directly cited from, 19 L v Si,P and L v Ge,P of both liquid and diamond solution phases have been optimized in this work. For the Ge-P-Si ternary system, no ternary interaction parameters (L i,j,k ) are evaluated in the present work since no experimental data are available. The ternary molar Gibbs energy Tern.ex G is extrapolated from contributions of the three binary subsystems. To add these contributions together, several methods have been developed by applying certain geometric principles, for instance, Muggianu method, Kohler method and Toop method. [43][44][45][46] In this work, Ge and Si are chemically similar, however, P has very different chemical behavior concerning chemical behaviors of components Ge and Si. Therefore, an asymmetric method (i.e. Kohler-Toop module) has been adopted as recommended by Refs. 43-46. Here, P is defined as the Toop constituent. Ge and Si are specified as the Kohler constituents. Accordingly, the excess Gibbs energy of ternary solution phase is written as: Regarding the three intermetallic compounds, descriptions of SiP and SiP 2 are directly cited from, 3 and that of GeP is thermodynamically assessed in this work, as described by following: Evaluated expression of C p above room temperature in Refs. 32-34 generate slightly smaller values (40.79 J/ mol K 32 and 42.89 J/mol K 33,34 ) than the experimentally determined C p of GeP at 298.15 K (44.35 J/mol K) by Ref. 31. In this work, the evaluated expression from 33,34 has been directly adopted, as the discrepancy between evaluated values and measured data was acceptable. Formation enthalpy of GeP at 0 K and normal atmosphere has been widely studied via DFT calculations in several open database platforms, such as Open Quantum Materials Database (OQMD), 37 Material Project 38 and Aflow, 39 as summarized in Table I. However, various lattice parameters and band gaps have been adopted to calculate the total energy of GeP by those platforms. This, together with the various pseudopentials of ground state of Ge and P, can partially interpret the large discrepancy among the theoretically calculated formation enthalpy (ΔH). As such, the comparison of these calculated data was impossible. On the other hand, the DFT calculated

Table I
GeP crystal structure (stable phase with space group of C2/m), and thermodynamic properties at 0 K calculated by DFT calculation Formation energy is described by   Table II. Accordingly, the measured ΔH 0 298 and S 0 298 at 298.15 K in Refs. 29 and 31 have been utilized to derive the Gibbs energy of GeP. In order to reproduce the dissociation of GeP, the formation enthalpy of GeP at 298.15 K has been slightly modified.

Density Functional Theory Calculation
Density functional theory calculations have been applied to find the formation energy of the GeP compound, monovacancy formation energy of Si, substitutional solute energies of dopants (P and Ge) in Si, and the 1st, 2nd, and 3rd neighbored di-doping formation energies. The calculation was performed by using the projector augmented wave (PAW) method, 49,50 as implemented in the highly efficient Vienna ab initio simulation package (VASP). 51,52 GGA was used for the exchange-correction functional. In this work, total energy of GeP (with a space group of Cs/m) has been calculated by using crystal constants from Ref. 48 It is reported as high quality in ICSD. Regarding energies related to Si_diamond structure, two sets of cubic supercells, 64 atoms (2 × 2 × 2) and 216 atoms (3 × 3 × 3), have been utilized in relevant simulations. The total energies of each supercell have been converged to less than 0.1 meV with respect to electronic, ionic and unit cell degrees of freedom by using the constant planewave energy cutoff of 600 eV. The 7 × 7 × 7 and 3 × 3 × 3 k-point meshes were chosen to integrate the supercell structure of 64 and 216 atoms over the Brillouin zone, respectively. These k-points were determined to be sufficiently large to perform the relevant calculations.
In order to gain insight into the Ge impact on the diffusional difficulty of P in a Si matrix, the climbing image nudged elastic band method (CI-NEB) 53,54 has been employed to find the so-called minimum-energy pathways (MEPs) for P migration as a sequence of intermediate configurations. The supercell with 64 atoms was used for the 1NN P atomic jumping geometry calculations. The ground state energies were determined by GGA level of theory. Five intermediate images between nearest local-energy minimum configurations were used in NEB calculation to find the energy barrier for the nearest neighbor P diffusion with/ without co-existing Ge.

Thermodynamic Evaluation of the Ge-P-Si System
In this work, the thermodynamic optimizations were conducted on the Ge-P system and the Si-rich side of the Si-P system. The experimental data of phase equilibria in Ge-P from Refs. 24 and 25 were measured under different pressures. Thus, the current optimization of the Ge-P phase diagram was started by reproducing the more reliable decomposition pressure of GeP. 29,30 Subsequently, the parameters of the liquid phase were optimized based on the experimentally measured invariant reactions, i.e. three-phase equilibrium of liquid-Ge-gas and liquid-GeP-gas. 31 The solubility of P in solid Ge was the last term brought into the optimization. Figure 1 illustrates the decomposition pressure of GeP obtained in this work, along with the experimental data from Refs. 29 and 30. The calculated decomposition temperatures under various pressure conditions fit well with the experimentally determined values. The measured three-phase equilibria (i.e. liquid-Ge-gas, and liquid-GeP-gas) were depicted on the pressure-temperature and pressure-composition projections. As shown in Fig. 2, a satisfactory agreement has been obtained between those experimental results and calculated data of this work.  1 Decomposition pressure of GeP: GeP(s) = Ge(s) + P (gas). X-axis is the total phosphorus pressure, with majority gas species of P 2 and P 4 . 1 3 Figure 3a and b show the calculated phase diagrams of Ge-P along with the published data of phase equilibria. [24][25][26][27][28]40 Less weight was given for data determined by microhardness in Ref. 26 that are much higher than all the others. 27,28 The phase diagram at 4600 kPa (Fig. 3a) has been calculated in which the liquidus evaluated by Olesinski et al. 40 was superimposed. Results from Ref. 24 were also plotted on Fig. 3a for providing available experimental points in the literature, as that has been done in Ref. 40. Comparing the evaluated liquidus by Olesinski et al., 40 the liquidus calculated in this work is slightly closer to the P edge. Nevertheless, such deviation should be ascribed to the different lattice stabilities of elements adopted in this work 41 and that reported in Ref. 40. This is evidenced by the different melting temperature of P_red at 4600 kPa. For instance, P_red had a melting point of 593°C in Ref. 40, while it is 578.9°C according to SGTE 41 adopted in this work. Figure 3b shows the Ge-rich side of the Ge-P binary phase diagram under normal atmosphere. The calculated single-phase boundary (SPB) of the solid solution phase, i.e. Ge_diamond (hereinafter, denoted as Ge(P)), agreed well with the measured solid solubilities of P in Ge_diamond. 27,28 It is evident that a set of self-consistent thermodynamic parameters of Ge-P was achieved here. Therefore, the phase diagram of Ge-P at normal atmosphere has been extrapolated based on the Ge-P thermodynamic description acquired in this work, as shown in Fig. 3c.
Regarding the solid solubility of P in Si_diamond (hereinafter denoted as Si(P)), there are extensive thermodynamic calculations and experimental measurements. 2,3,21-23,55-60 The determined P solubility boundaries in Si solid scattered largely. The experimental phase equilibria have been critically reviewed by Jung et al. 2 and Liang et al. 3 Hereby, no space is taken for repeating the review in this paper.
As described in "A Brief Review of Previous Work on Phase Equilibria" section, the Si-rich side of the Si-P system requires a re-optimization owing to the reliable lattice stability of metastable P_diamond. In this work, the thermodynamic parameters of Si(P) and liquid solution have been re-optimized with reference to the work of Liang and Schmid-Fetzer 3 owing to their insight evaluation into the original experimental works. Here, the author also initially attempted to divide the measured P solubility limits in Si_ diamond into two groups as the strategy applied by Liang and Schemid-Fetzer. 3 Accordingly, two sets of parameters (Model I and II) were achieved to reproduce the scattered single-phase boundary of the Si(P) solution phase. As illustrated in Fig. 4, these calculated P solubility boundaries are in satisfactory agreement with each group of the reported data, except for data from Ref. 26. The P solubility in Si_diamond determined by microhardness measurements 26 represented the lowest solubility boundary among those reported data. During this optimization, a small weight has been given to data from Ref. 26

3
On the other hand, the experimental data, representing the highest solubility limits, 56,57,60 was rooted from the Si-P-O ternary system rather than the pure Si-P binary system. The oxygen originated either from the starting material 56,57 or the subsequent annealing process. 60 For instance, the slightly oxidizing atmosphere (90% N 2 + 10% O 2 ) was applied during the high-temperature heat treatment processes in Ref. 60. The Hall measurement is less accurate for hightemperature species since the concentration of electrically inactive phosphorus increases with temperature due to the phosphorus diffusion in silicon. This phenomenon becomes significant when the temperature is above 750°C. 60 Consequently, we adopted a third model to reproduce the P solubility boundary in Si solid, i.e. the solvus determined by Hall measurement (at relatively low temperature) 58,59 and the solidus determined by DTA. 22 The calculated boundary was superimposed in the inserted graph of Fig. 4. The published experimental data were relatively well reproduced. Model III thus represented the most reliable set of description in this work. It has been utilized in the thermodynamic extrapolation of the Ge-P-Si equilibrium phase diagram.
In the CZ Si growth process, the typical gas pressure is in the range of 1.5-5 kPa. 1 The decomposition temperatures of 1:1 ratio phosphide in both Ge-P and Si-P systems decrease with lowering total pressure. It is, therefore, of great importance to extrapolate the phase diagrams of Ge-P and Si-P to relevant process pressure, for instance 5 kPa. Figure 5 shows the calculated phase diagrams of Ge-P and Si-P at 1.5 and 5 kPa. The diamond solid and gas two-phase equilibrium region presented in both systems. The three-phase

3
equilibria of phosphide-liquid-diamond were replaced by phosphide-gas-diamond in both systems as pressure lowered to 5 kPa. The temperature window for liquid-diamond two-phase co-existing regions shrank under such low-pressure condition, so did the solid solubility limit of P in both Si(P) and Ge(P) solids. These phenomena get more pronounced with further decreasing pressure, as shown in the relevant phase diagrams at 1.5 kPa in Fig. 5b and d. These issues complicate the control of the growth process of heavily P-doped CZ Si single crystal. For instance, the evaporation of P starts at lower temperature when a lower pressure condition is applied.
For the Ge-P-Si ternary system, it was thermodynamically extrapolated based on the thermodynamic descriptions of the binary end-member systems, i.e. Ge-P was according to SGTE, 19 Ge-P and Si-P were according to the current results. As no experimental investigation available, analytical methods were not attempted to assess the ternary interaction parameters. The geometric model has been utilized to extrapolate the thermodynamic description of the ternary solutions (i.e. Si(Ge, P) and liquid phases). In this work, the Kohler-Toop model was the optimal model, in order to treat the very different chemical behavior of P from those of Ge and Si. The obtained thermodynamic parameters of this work are listed in Table III.

Impacts on Heavily P-doped Si during Growth
The solidus lines of Si in Fig. 5c and d indicate that the maximum equilibrium P solubility in Si is around 0.478-0.68 at.% under the typical gas pressure range (1.5-5 kPa) of the CZ Si growth process. These mole fraction values correspond to a P atom density of around 2.39 × 10 20 -4.85 × 10 20 atoms/cm 3 . The conversion between mole fraction (x p ) and atom density (X p ) is as follows: where N A is the Avogadro constant, M Si 1−x P x is the atomic mass of the Si solution in g/mol that is derived from atomic mass of P and Si in this work, i.e. M Si 1−x P x = (1 − x)M Si + xM P . The density of Si, Si , in g/cm 3 is adopted for the conversion, for the following reasons: (i) negligible effect of P dopant level (i.e. ≤ 10 18 atoms/cm 3 ) on the Si density 12 and (ii) no experimentally measured density of Si at high P doping level.
Despite the small solid solubility of P in Si, the P-doped CZ Si crystal with resistivity under 1 mΩ cm has been achieved. For instance, Chiou 9 incorporated about 1.11 × 10 20 atom/cm 3 P in CZ Si crystal and obtained a resistivity value of 0.71 mΩ cm. The equilibrium P concentration in Si solid presented in Fig. 5c and d suggests the potential to further lower the resistivity by incorporating higher P level in Si crystal. For example, to achieve a resistivity value of 0.5 mΩ cm, it requires dissolving around 1.65 × 10 20 atom/ cm 3 (or 0.33 at.%) P in Si, according to the free-access resistivity-composition converter. 62 This P doping level creates a lattice contraction of around 2.97 × 10 −4 , as the lattice strain of substitutional solution is 9,12

3
Here Δa is the modified lattice constant with increasing dopant concentration, a Si is the lattice constant of undoped Si (0.5431 nm). For dopants occupying the substitutional sites, the solute lattice contraction (+) or expansion (−) coefficient β is determined by the covalent radii of dopants (r dopant ) and that of Si (r Si ), as well as the number of Si atoms per cm 3 (X Si is 5 × 10 22 atom/cm 3 ) 13 The covalent radii of Ge, P and Si are: 0.122 nm, 0.11 nm, 0.1173 nm, respectively. 12,13,16 The coefficients β of Ge and P in Si are hereby obtained: β P is around 1.17 × 10 −24 cm 3 /atom, and β Ge is around − 8.91 × 10 −25 cm 3 /atom. Using these coefficients of β, the concentration of Ge that is required to compensate the lattice strain introduced by 1.65 × 10 20 atoms/ cm 3 P in Si can be derived via X Ge β Ge = X P β P . Here, the required concentration of Ge is around 2.17 × 10 20 atoms/ cm 3 (or around 0.45 at.%, referring to Eq. 6).
Using the extrapolated thermodynamic database of Ge-P-Si ternary system, isopleth along 0.33 at.% P at Sirich side of Si 99.67 P 0.33 -Ge 99.67 P 0.33 (at.%) has been calculated. As presented in Fig. 6a, at constant P doping level (0.33 at.%) the solidification temperatures of the Si(Ge, P) ternary alloys decreased with increasing Ge contents. For instance, a slightly lower solidification temperature (1394.25°C) was obtained for SiGe 0.45 P 0.33 comparing with that of the binary alloy SiP 0.33 (1399.95°C). In addition, the single phase boundaries of Si(Ge, P) in the Fig. 5 Calculated Ge-rich side of the Ge-P phase diagram at 5 kPa (a) and 1.5 kPa (b), as well as calculated Si-rich side of the Si-P phase diagram at 5 kPa (c) and 1.5 kPa (d). Here Ge(P) and Si(P) represent the solid solution phase of Ge_diamond and Si_diamond in Ge-P and Si-P binary system, respectively.

3
isopleth also reveal that for 0.33 at.% P the incorporated Ge concentration limits at 3.369 at.%. For a binary system, it is straight forward to get the K eq value from the relevant phase equilibrium, i.e. K eq = C L /C S , C L and C S are the dopant concentrations on the liquidus and solidus lines at a given temperature. Concerning the K eq of P in Si(P, Ge), phase diagram calculation works as well. Firstly, a vertical section that crosses the doping concentration of interest is calculated. The solidification temperature of that alloy is then obtained. Subsequently, a corresponding isothermal section will be calculated. The liquid composition that equilibrates with the solid can be directly taken by selecting the relevant tie-line in the liquid-solid two-phase equilibrium field (as depicted in Fig. 6b), and then the K eq of the dopant at his doping level (e.g. 0.33 at% P in this work) and temperature is acquired. Table IV lists the K values from literature and the current calculation by using the extrapolated Ge-Si-P ternary thermodynamic description.
The K eq values of P in Si at around melting point of Si obtained in this work are comparable with other theoretical data. 3,22,58 However, all the calculated values are way too small than the widely accepted effective segregation coefficient (K eff of 0.35). 56 The reason of such difference can be partially interpreted by the well-known Barton-Prim-Slichter equation 63 and the modifications. 4,64 The K eff is determined not only by the K eq but also by the concentration dependent properties of the dopants (e.g. diffusivity in melt, evaporation rate) and the process parameters (e.g. pulling rate, angular velocity of crystal rotation). In addition, the impurities of oxygen and carbon are unavoidable during CZ Si growth process. The dissolved oxygen impacted the solubility level of P in Si, 65 as also discussed in the last section. Furthermore, the thermal donors, associated with interstitial oxygen in CZ Si and formed during cooling process in temperature range of 400-500°C, contribute to the free carrier concentration, 66 thus influencing the accuracy derivation of dopants K values via resistivity methods. The influence of oxygen on the heavily P-doped CZ Si will be exploited in another publication.
There is limited data available regarding the K eff of P in Si x Ge 1−x 16,67-69 Refs. 67 and 68 reported the K eff values of lightly P-doped Ge-rich SiGe crystal (Si 0.05 Ge 0.95 ); however, the P concentration was not specified. Yonenaga et al. 69 evaluated the K eff of P in Si-rich CZ SiGe crystal based on (i) the carrier concentrations in the melts and the grown crystal via relation K eff = C S /C L , or (ii) the carrier concentrations at the position of growth initiation of the crystal (C S0 ) and carrier concentration at grown SiGe with a solidification ratio (f) of 0.5 (C S ) via normal freezing equation 16 They reported a dramatic increase of the K eff (listed in Table IV). 16 Those values were not without doubt since the oxygen concentration was unignorable, for instance, it was 9 × 10 17 atoms/cm 3 in their un-doped crystal. 16 In this work, the comparison with K eff values reported in Ref. 16 were not attempted. Instead, the K eq of P varying with the Ge fraction in Si 1−x Ge x was calculated (Fig. 7) by referring to the alloy concentrations utilized in Ref.  16. The results demonstrated that: (i) at constant P doping level, the K eq of P decreased with increasing Ge fraction in the Si 1−x Ge x crystal; (ii) at a constant Ge fraction, the higher the P doping level, the smaller K eq of P one gets. As shown in Fig. 7, the K eq values of P are greatly dependent on the Ge fraction. It is observed that the K eq value of P in Si(P, Ge) can increase with increasing P concentration as long as Ge fraction changes accordingly, for instance, under the condition that x(Ge) is less than 0.043 at the P doping level of 6 × 10 19 atoms/cm 3 , and x(Ge) is higher than 0.046 at the P doping level of 8 × 10 17 atoms/cm 3 .

3
The K eq values of P (0.33 at.%) with respect to different Ge co-doping level were theoretically calculated using the currently extrapolated Ge-P-Si thermodynamic database. The relevant solidification temperatures were directly read from the isopleth of 0.33 at.% P (Fig. 6a). At this P doping level (0.33 at.%), the K eq values of P in Si(P, Ge) decreased with increasing Ge co-doping concentrations.

Impacts on P Distribution in Si Matrix
To understand the role of Ge in P distribution in the Si matrix, the current DFT calculations simulated the monovacancy formation energy of un-doped Si, the formation energy of substitutional Si(P) and Si(Ge) solutions, as well as the binding energies of two like (P-P, Ge-Ge) and dislike (Ge-P) dopants positioning as 1NN, 2NN and 3NN in the Si matrix. Generally, the formation energy can be derived from the equation: Accordingly, the vacancy formation energy of Si without solute atoms ΔE V f is determined by the energy difference of the perfect structure and that of monovacancy: where N is the total atoms in the supercells, Va 1 indicates one vacancy on the Si lattice site. E Si N−1 Va 1 and E Si N correspond to the total energy of the supercell including one vacancy and that of perfect cell.
The formation reaction for the Si-based solid solution with substitutional dopants is: One can obtain the substitutional formation energy of a single dopant in Si matrix: Here, E Si N−1 X 1 is the total energy of the cell including one dopant, and E X is the energy of an isolated dopant atom. Table V summarizes the formation energies derived from total energies obtained by using GGA exchange-correlation potential. The monovacancy formation in Si obtained in this work is comparable with data from previous DFT studies, 75,76 i.e. at 0 K for the neutral vacancy,ΔE  Table V show that of the calculated formation energy of the Si(Ge) substitutional solution is less than that of Si(P). It implies a higher tendency to form the Si(Ge) solution phase than to form the Si(P) solution phase, which is supported by the significantly larger solubility of Ge in Si solid phase comparing with that of P. The binding energy for two substitutional solute dopant atoms in Si is obtained from: In Eq. 13, iNN denotes the 1st to 3rd nearest neighbors. X 1 and X 2 are the substitutional dopant atom Ge and P. Table VI lists the binding energies of P-P, Ge-Ge and Ge-P derived by following Eq. 13. The binding energy between Ge and Ge atoms at the 1NN position is smaller than the binding energies of Ge-Ge at 2NN and 3NN positions. It is apparent that Ge atoms tend to gather in the Si matrix. Examining the results of the binding energies of Ge-P bonds, one can notice that Ge-P bonds distributing at 1NN and 2NN positions are more energetically favorable than at 3NN; however, for P-P bonds the opposite tendency is observed. In addition, binding energies of P-P bonds are higher than those of Ge-P and Ge-Ge bonds for each level of neighboring bonds position. This suggests that the P atoms are energetically favorable to interact with Ge atoms rather than to cluster in the Si matrix. The Ge co-doping in P-doped Si can then provide the drag-and-drop effect, and thus to homogenize the P distribution in Si solid. P diffusion in Si poses a technological difficulty in forming well-defined doped regions for devices. In order to check the impact of Ge on the migration of P atoms in the Si matrix, the CL-NEB calculations were performed to simulate the migration barriers for a 1NN P jump in the 64-atom Si cubic supercell. The calculations were executed with/without a Ge atom occupying a site that also forms the nearest neighbor bond with the same corner atom as P neighboring. The barrier energies of the P diffusion between 1NN position are depicted in Fig. 8a. The migration of P atoms in the Si matrix is illustrated in Fig. 8b. Curves in Fig. 8a demonstrates that the barrier energy of P diffusion with Ge presented is larger than that of the case without involving Ge, which is in agreement with previous experimental research on stability of the P-Va pair in n-doped Si 1−x Ge x . 77 This justifies our conclusion that Ge co-doping can slow P diffusion in Si solid.

Conclusion
This work investigated the Ge influence on the P distribution in Si based on thermodynamic considerations by analyzing results from CALPHAD approaches and DFT calculations. Phase diagrams of Ge-P and Si-P binary systems were thermodynamically re-assessed. Thermodynamic description of the Si-P-Ge ternary system has thus been theoretically extrapolated. It was then utilized to calculate the isopleth diagram of SiP x -GeP x (x = 0.33 at.%) and isothermal sections of Si-Ge-P at 5 kPa. Accordingly, the Ge impacts on the solidification temperature and K eq of P in doped Si have been studied based on the calculated diagrams. It was found that at constant P concentration (0.33 at.%) the solidification temperature of Si(Ge, P) decreased with increasing Ge concentration, so did the K eq of P. DFT calculations performed on the formation energy for monovacancy of Si crystal and substitutional Si-dopant solutions, i.e. Si(Ge) and Si(P), and the binding energy of P-P, Ge-Ge and Ge-P bonds loacted at 1NN, 2NN and 3NN positions. Comparing with P-P bonds, both Ge-Ge and Ge-P bonds have larger binding energy at all the three configurations. The energy barriers of P migration between the 1NN bond were calculated by DFT-NEB. The higher energy barriers were obtained when neighbored by a Ge atom than that without Ge occupying a site of another1NN bond. These results led to the conclusion that in heavily P-doped Si, Ge does not have the capability to improve the even distribution of P along the Si ingot during crystal growth process, nevertheless, Ge co-doping in the solid state of the P-doped Si crystal can reduce the P atoms gathering and its diffusion in Si and thus to stabilize the P-doped Si structure.
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

3
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.