Nitrogen Dissolution in Liquid Ga and Fe: Comprehensive Ab Initio Analysis, Relevance for Crystallization of GaN

The dissolution of molecular nitrogen in Ga and Fe was investigated by ab initio calculations and some complementary experiments. It was found that the N bonding inside these solvents is fundamentally different. For Ga, it is between Ga4s and Ga4p and N2p states whereas for Fe this is by N2p to Fe4s, Fe4p and Fe3d states. Accordingly, the energy of dissolution of N2 for arbitrarily chosen starting atomic configurations was 0.535 eV/mol and −0.299 eV/mol for Ga and Fe, respectively. For configurations optimized with molecular dynamics, the difference between the corresponding energy values, 1.107 eV/mol and 0.003 eV/mol, was similarly large. Full thermodynamic analysis of chemical potential was made employing entropy-derived terms in a Debye picture. The entropy-dependent terms were obtained via a normal conditions path to avoid singularity of ideal gas entropy at zero K. Nitrogen solubility as a function of temperature and N2 pressure was evaluated, being much higher for Fe than for Ga. For T=1800 K and p=104 bar, the N concentration in Ga was 3×10−3 at. fr. whereas for Fe, it was 9×10−2 at. fr. in very good agreement with experimental data. It indicates that liquid Fe could be a prospective solvent for GaN crystallization from metallic solutions.


Introduction
Growth of large-size single crystals of GaN has been an important issue for the development of semiconductor technology for several decades. That is because GaN is difficult to grow, both from vapor and liquid phases, and to dope effectively to n-type, p-type and semi-insulating. A number of obstacles contribute to these difficulties, including thermodynamic and kinetic factors [1].
GaN is required for the manufacture of a wide variety of optoelectronic and electronic semiconductor devices, such as optoelectronic light-emitting diodes (LEDs) and laser diodes (LDs) [2], to recently developed THz emitters [2], or in electronic high-power transistors, high electron mobility transistors (HEMT) [3], or 2-D gas (2DEG)-based molecular sensors [4]. The performance of these devices critically depends on the quality and the properties of substrates used for epitaxy. Thus, the goal is to develop reliable and cost-efficient methods of crystallization of large-size excellent-quality single crystals. Despite considerable progress attained using several methods, including gas-phase hydride vapor-phase epitaxy (HVPE) [5,6], ammonothermal [7], or solution growth [8], the goal was only partially achieved.
Generally, liquid-phase methods are more suitable for providing high-quality crystals than those based on the vapor source, as the temperature gradients are easier to control in the liquid, the control of point defects is better, and the chemical purity is easier to protect [9,10]. The gas-phase growth of GaN does not include direct synthesis from the elements, as the kinetics blocks the process thus additional transport agents contribute to uncontrolled doping. The liquid-phase growth includes ammonia-based and metal-based

Calculation Method
Density functional theory (DFT)-based code SIESTA [19,20] was employed in all ab initio calculations reported in this study. SIESTA is the code developed under Spanish Initiative for Electronic Simulations with Thousands of Atoms which was designed to simulate systems including thousands of atoms. In the present calculations the more precise version, capable of dealing with smaller number of atoms, was used. The solution procedure used in Generalized Gradient Approximation (GGA) employs norm-conserving pseudopotentials for determination of the system wave-function created using linear combinations of local basis molecular orbitals. All electron calculations are not possible, therefore the norm-conserving Troullier-Martins pseudopotentials, in the Kleinmann-Bylander factorized form [21,22] for Ga, Fe and N atoms were generated using ATOM program written by the authors of the SIESTA code. The exchange-correlation functional was adopted in revised modification of Perdew, Burke and Ernzerhof (PBE) [23,24] functional for solids and surfaces (PBEsol) [25,26].
The cutoff value in the equivalent plane wave for the real space properties and integration grid maximum values was 275 Ry, roughly equal to the grid spacing in the representation of the distribution of any quantity in position space of about 0.1 Å. The solution of Kohn-Sham equations uses finite size atomic orbitals, therefore it is not affected by the grid size. This grid selection affects the real space representation of the properties derived from the grid. The self-consistent field (SCF) loop was terminated when the maximum difference between the output and the input of each element of the density matrix was below 10 −4 . In the relaxation procedure, the positions of atoms were modified when the forces acting on these atoms were higher than 0.005 eV/Å.
The data obtained in this paper also include an assessment of the number of configurations that are used in the entropy determination. As it is shown, not every configuration is considered, thus the results are also burdened by some systematic error. Thus, the entire chemical potential dependence includes the systematic errors that are enhanced in the evaluation of the equilibrium pressure.
The simulations of interaction between Ga, Fe and N needed critical tests. One of the possible tests was a comparison of DFT results and the experimental data for gallium nitride. The lattice constants following from the ab initio total energy minimization of the wurtzite crystalline lattice of GaN were a = 3.194 Å, c = 5.186 Å. These values are in a reasonable agreement with the experimental data for GaN a = 3.1890 ± 0.0003 Å, c = 5.1864 ± 0.0002 Å [27].
In addition, the DFT ab initio simulations of the nitrogen molecule (N 2 ) characteristics could serve as the test of compatibility of the DFT parameterization for molecular species. This is important as the typical choice for simulation of solids and semiconductors in particular, works poorly for small, covalently bonded molecules. The compromise solution is naturally inferior from the ones optimized for molecules only, and naturally less precise than advanced ab initio calculations, such as W1 or CCSD(T) [28][29][30]. The test values for the N 2 molecule, obtained in our approach were dissociation energy ∆E diss and the bond length d. The N 2 dissociation energy obtained in the presently used parameterization is ∆E diss DFT (N 2 ) = 9.801 eV which is in a good agreement with the experimental value ∆E diss exp (N 2 ) = 9.790 eV [30]. The DFT bond length d N-N DFT = 1.092 Å was in an excellent agreement with the experimental value d N-N exp = 1.097 Å [28]. An additional insight may be obtained from the diagram of bonding states in N 2 molecule, presented in Figure 1. The Crystal Orbital Hamiltonian Population (COHP) data is used to determine the bonding/antibonding interactions [31]. The data presented in Figure 1 indicate that the nitrogen N2s bonding and antibonding states have the lowest energy and accordingly, they are occupied. In addition the bonding N2p states (both 2π and 2σ bonds) are occupied as well so that they contribute to high dissociation energy of the molecule. The integration of COHP peaks gives the following data: (i) E N2s−N2s = 4.002 eV; (ii) E N2s−N2p z = 2.366 eV; (iii) E N2p z −N2p z = 1.623 eV (iv) E N2p xy −N2p xy = 1.616 eV. These values of the sum of the covalent interaction overlaps is 9.607 eV which is in good agreement with the N 2 dissociation energy ∆E diss exp (N 2 ) = 9.790 eV [30]. As a next testing step, basic physical properties of pure gallium and iron were determined. These properties include the cohesive energies, the band structure and the work functions. The atomization energy of gallium, obtained from ab initio calculations was ∆E atom DFT (Ga) = 2.898 eV/atom which is in a reasonable agreement with the atomization energy of rhombic Ga structure, equal to ∆E atom exp (Ga) = 2.81 eV/atom as determined in reference [32]. The density of Ga was evaluated as equal to ρ DFT (Ga) = 6.304 g/cm 3 , which is favorably compared to ρ exp (Ga) = 6.095 g/cm 3 at the melting point (close to 300 K). In summary, the properties of metallic Ga are reasonably recovered by the parameterization of Ga atom orbitals used in SIESTA.
Similar analysis could be made for iron bonding properties. The atomization energy of Fe, obtained from spin-polarized ab initio calculations, was ∆E atom DFT (Fe) = 4.978 eV/atom which is in reasonable agreement with the experimental value of ∆E atom exp (Fe) = 4.28 eV/atom [32]. The DFT density of Fe was equal to ρ DFT (Fe) = 7.825 g/cm 3 , which is not far from ρ exp (Fe) = 7.874 g/cm 3 of iron at normal conditions. For evaluation of Ga and Fe metal properties, the supercells have been chosen. They are shown in Figure 2. Gallium was simulated using 144 Ga atoms supercell ( Figure 2a) representing orthorhombic solid Ga lattice , the most stable Ga structure at low temperatures. Iron was simulated using 108 atom supercell representing face centered cubic crystalline lattice. A cluster of 108 atoms of pure Fe is similar to the cluster shown in Figure 2c with the difference that the Ga atom is substituted by the Fe atom (substitution of the green ball to the gray ball).

Experimental
For experimental evaluation of N solubility in Fe at high nitrogen pressure (HP) and at high temperature (HT), a series of annealing experiments in HP-HT gas (N 2 ) pressure reactors was performed. In this study, the pressure dependence of N solubility in the Fe metal liquid, with the N 2 gas as a source of the solute, has been determined.
For this purpose, a series of iron samples (ARMCO Pure Iron, Grade 2) in the form of cylinder (11 mm in diameter) inserted into BN crucible as shown schematically in Figure 3a was used. Initial mass for each sample (m 0 ) was measured. Then the crucible with the sample was placed into the two-zone furnace and loaded into the high-pressure chamber. The system was evacuated, filled with nitrogen (6 N) and the gas compressed until a required pressure. This preparatory procedure was used in all experiments. The annealing experiments were performed at quasi-isothermal conditions. Annealing temperature for all experiments was 1708 ± 1 K, duration of the annealing at constant temperature was 0.5 h. The system was heated with a rate 800 K/h and cooled by quenching (rate of cooling ca. 3000 K/h) to freeze the high temperature state of the Fe-N solution. The N 2 pressure sequence of 5000, 6000, 7000, 8000, 9800 bar was used to determine the pressure dependence of N solubility in the Fe metal. After the experiments, all samples changed their shape and mass (see scheme in Figure 3b). It was assumed that the increase of the sample mass was due to N dissolved in the metal. A background for such an assumption was our previous analysis reported in references [15,33]. For the reported analysis, the composition of the Fe-N samples treated in a similar way as in this work was studied by the Energy Dispersive X-ray Spectroscopy (EDX-ZEISS SEM with BRUKER detector, Oberkochen, Germany) and the Inert Gas Fusion spectroscopy (ON836 LECO, St. Joseph, MI, USA). We have shown that the observed increase of the sample mass corresponded exactly to the mass of dissolved nitrogen determined by the indicated spectroscopy methods.
The results obtained in this work were used in Section 3.5 for comparison to theoretically calculated solubility values. For iron, the electronic bonding may be analyzed by plots of DOS of a single Fe atom and the Fe supercell consisting of 108 atoms. The results are compared in Figure 4c,d. As it is shown, bonding between the iron atoms arises mainly from the extensive overlap of the Fe3d orbitals. This is drastically different from Ga where 3d orbitals do not participate in the bonding. In summary, the properties of metallic Fe show important differences with respect to the Ga bonding.

Results and Discussion
Finally, a single Ga atom immersed in the Fe matrix was modeled using total 108 atom cluster, composed of 107 Fe atoms and of a single Ga atom. The cluster with Ga atom incorporated is presented in Figure 2c whereas the cluster of pure Fe is inferred from Figure 2c   The results of the calculations are shown in Figure 4e. As it is shown, the 3d states of Ga atoms do not participate in the bonding with Fe matrix. Also Fe3d states are not involved. The bonding is due to interaction of Fe and Ga 4s and 4p states. The difference is due to different behavior of Ga3d states which remain unaffected by surrounding Ga and Fe neighbors. This is demonstrated by the PDOS of Ga/107 Fe SC presented in Figure 4e where Ga3d preserved its molecular character (sharp line), not affected by the overlap with the neighboring Fe atoms.
The energetic effect of dissolution of single Ga atom in Fe was calculated taking into account that an Fe atom is replaced by a Ga atom, according to the formula: where E DFT (107 Fe-1 Ga) denotes total DFT energy of the cluster of atoms (numbers and symbols denote the cluster size and the type of constituting atoms). From the performed DFT calculations, the energy of Ga dissolution in the Fe matrix was determined as ∆E dis DFT (Fe-Ga) = 3.698 eV/atom. The value is relatively large indicating strong interaction between Ga and Fe atoms. The experimental data on liquid Fe:Ga solutions confirm qualitatively, these results by a drastic decrease of Fe melting temperature induced by small concentrations of gallium in iron [15].

Interaction of N Atom Immersed in the Metal (Ga, Fe) Periodic Clusters
Similar analysis could be made for dissolution of nitrogen in both Fe and Ga solvents. As shown recently by Ponomareva et al., the calculated N dissolution energies in a solid Fe cluster depend on the cluster configuration [34]. The pure Fe cluster could be inferred from Figure 2c by replacement of Ga by Fe atom while Fe cluster with a single immersed N atom is presented in Figure 2d. These clusters are used in our ab initio calculations.
The ab initio calculation results for nitrogen atom immersed in the Fe matrix are presented in Figure 5. As it is shown the Fe atoms in the matrix are bonded by Fe3d orbitals while the N atom could be bonded by its own N2s and N2p orbitals. The bonding has complex character, the N2s overlap with Fe4s and Fe4p states is bonding while N2s states are bonded to Fe3d states very weakly. Thus, the bonding of N2s states is essentially absent. In contrast to that, the bonding of N2p states is strong, to both Fe4s and Fe4p and to Fe3d orbitals. In summary, the interstitial nitrogen atom bonding to the Fe matrix is essentially only via N2p states. Please note that the approach is different from the one applied to the Ga-Fe cluster. The number of Fe atoms is preserved, and a single N atom is added into an interstitial position. Therefore, the energy of dissolution of single nitrogen atom in the Fe matrix is calculated using equation different than Equation (1): Using the calculated DFT values for the right-hand side of the equation, the energy of atomic N dissolution in the Fe matrix is: ∆E dis DFT (Fe-N) = −5.050 eV/atom. Noting that the energy of molecular nitrogen is lower by the molecule dissociation contribution ∆E diss DFT (N 2 ) = 9.801 eV, the energy of dissolution of single N 2 molecule could be obtained as: Hence the resulting energy of dissolution of molecular nitrogen in Fe is ∆E dis DFT (Fe-N 2 ) = −0.299 eV/mol. The negative value, indicates a reduction in energy at dissolution of molecular nitrogen, thus promoting high solubility of nitrogen in Fe.
The cluster used for ab initio calculations describing N interaction with the Ga matrix is presented in Figure 2b.
The results of the calculations for nitrogen atom in the Ga matrix are presented in Figure 6. As it was already shown, gallium matrix atoms are bonded by Ga4s and Ga4p orbitals while N atom could be bonded by its own N2s and N2p orbitals. The N-Ga bonding is drastically different from the N-Fe case. The N2s bond with Ga states is molecular in character with small energy dispersion and small magnitude while N2p states are bonded by large dispersion state, i.e., they are extended in real space- Figure 6d. There is a small contribution of N2p-Ga3d bonding, due to extended nature of N2p state- Figure 6f. The bonding and antibonding overlaps between N2s and Ga3d states are compensated. Therefore as it is shown above, bonding of interstitial N atom to Fe and Ga matrices has different character. Please note that the approach is similar to the one applied above to the Fe-N cluster. The number of Ga atoms is preserved, and a single N atom is added. The energy of dissolution of nitrogen single atom in the Fe matrix is calculated using Equation (2). From the obtained DFT values, the energy of atomic N dissolution in the Ga matrix is: ∆E dis DFT (Ga-N) = 4.633 eV/atom. Then the energy of dissolution of molecular nitrogen is obtained via Equation (3). Thus, the energy of dissolution of molecular nitrogen in Ga is: ∆E dis DFT (Ga-N 2 ) = 0.535 eV/mol. The energy is positive and significantly higher than for Fe, therefore the energy increase should be observed at dissolution of molecular nitrogen, and accordingly, the solubility of nitrogen in liquid Ga should be drastically lower than in Fe.
It was found however that the obtained total energy of the Ga-N system changes considerably, depending on the configuration used [34]. Therefore, in order to average over configurational degrees of freedom and obtain representation relevant for the liquid, the molecular dynamics (MD) ab initio simulations of the 54 Ga atom cluster, with and without a single N atom immersed in, have been performed. In such a model, the corresponding concentration of nitrogen was 0.01851 atomic fraction (at. fr.) The average temperature of the simulation was set to 500 K. The configurations of the system, the initial one and the ones after prescribed number of MD steps, is presented in Figure 7. In fact, the temperature 500 K is above the melting point of Ga (T M Ga = 303 K at normal pressure). Therefore, after a short initial period, the configuration of atoms is changed completely. Thus, the initial configuration has no influence on the time averaged data. The time evolution of the total energy of the system for the 54 Ga atom cluster, with and without a single N atom immersed in, is presented in Figure 8. In the simulation, the averaging was made using a sequence of time steps until a steady state has been achieved, as shown in Figure 8. The averaging for the final result was undertaken when the error in the simulation was smaller than 0.001 of the total energy value. Due to periodic boundary conditions, the system is not closed. Therefore, the energy of the system can change in time. The initial time evolution of the Ga cluster both with and without N atom is similarly fast. The system attains the energy, close to the average one after several hundred MD steps. As it is visible, the energy of the system fluctuates around the average one. The temporal evolution of the fluctuations is different in both cases. The system without N atom exhibits two different fluctuation types: first short lived one, typical for thermal noise of independent thermal motions of single atoms, and the second one which is long lived fluctuations typical for hydrodynamics, i.e., common motion of higher number of atoms. Due to sheer size of the simulation cell, this number cannot be too large. Anyway, such a phenomenon is clearly observed.
The system with the immersed N atom behaves differently. The thermal fluctuations are present with the noise similar to the previous one. Nevertheless, the hydrodynamic fluctuations are absent which may be attributed to the small size of the simulated system. Nitrogen atom strongly binds several atoms in the first coordination zone, hampering their motion. The remaining atoms are not able to move in coordinated manner, thus the system fluctuations are reduced to the thermal noise only.
The time averaged total energy was E MD DFT (54 Ga-1 N) = −125,169.879 eV and E MD DFT (54 Ga) = −124,896.443 eV for these two clusters, respectively. The DFT energy of single N atom is E DFT (1 N) = −269.422 eV. In order to derive the DFT energy of dissolution, it should be noted that the total energy derived from the MD simulations contains two additional terms, not present in the standard DFT calculations. The first one is the average energy of thermal motion of atoms, which for a single nitrogen atom, is ∆E T (1 N) = 3k B T, the value that arises from both kinetic and potential energy average. The second contribution is the zero-point energy E ZPE Therefore, the resulting energy of dissolution of molecular nitrogen in Ga(l): ∆E dis DFT (Ga-N 2 ) = 1.107 eV/mol. The energy change associated with single atom is half of the one for molecule, i.e., ∆E dis DFT (Ga-N 2 ) = 0.553 eV/atom. The energy is positive, therefore the energy increase occurs at dissolution of molecular nitrogen that leads generally, to relatively low N solubility in Ga as was already observed [11,35,36].
Similar MD simulations were made for N dissolution in Fe using 54 Fe atoms cluster with single N atom immersed in the metal as presented in Figure 9. As above, the concentration of nitrogen was 0.01851 at. fr. The technical details of the simulations were the same as these used for Ga and the temperature was again set at 500 K.
The evolution of the system looks differently than for the Ga case. At 500 K, Fe is in the solid phase, therefore the evolution of the system is typical for the crystal. The configurations shown in both Fe and Fe:N cases are typical for the solid phase where the long-range order is preserved, and the atoms are attached to their lattice sites. The presence of nitrogen does not induce melting, the order is still preserved. Nevertheless, the fluctuations change the lattice in a visible manner.
The time evolution of the system total energy for the Fe and Fe-N clusters, are presented in Figure 10. As before, the final result was undertaken when the error in the simulation was smaller than 0.001 of the total energy value.  As above, due to periodic boundary conditions, the system is not closed so that the energy of the system may change in time. The time evolution of the Fe cluster both with and without N atom is the same as for the Ga case, similarly fast. The system attains the energy, close to average one after several hundred MD steps. The energy of the system fluctuates around the average. The thermal fluctuations are present with the noise similar but larger for N atom present. The hydrodynamic fluctuations are absent as the system is essentially solid.
The time averaged total energy was E MD DFT (54 given by Equation (5) is ∆E dis DFT (Fe-N 2 ) = 0.003 eV/mol, and accordingly the part associated with single atom is: ∆E dis DFT (Fe-N 2 ) = 0.001 eV/atom. The energy is almost zero, i.e., much lower than the value obtained for liquid Ga, indicating that the solubility of nitrogen in Fe should be much higher at the otherwise, similar conditions.

Equilibrium between N 2 Gas and Nitrogen Dissolved in Ga and Fe Metals
The equilibrium between gas and liquid phases entails equality of chemical potentials. In case of dissolution of nitrogen from the diatomic form in its gas phase in a metallic liquid it is recognized that the N 2 molecule dissociates in contact with liquid metals. For Ga solvent, it was demonstrated by Romanowski et al. [12,13] by DFT modeling of both N 2 dissociative chemisorption and dissolution processes. As iron has higher bonding energy to nitrogen than gallium, it is expected that dissolution of nitrogen in liquid Fe is also dissociative that leads to universal relation for equilibrium: Direct comparison of the chemical potential of both species is not possible as there is no direct prescription to obtain chemical potential at given temperature and pressure. The chemical potential equality can be analyzed by calculating the enthalpy and entropy differences for nitrogen in the liquid and gas phases separately and combine them in the Gibbs free energy. That causes considerable difficulties as in the standard approach, the enthalpy difference at zero temperature, equal to the ab initio energy difference between the vapor and the liquid, is considered. The other contributions are treated as entropy term, i.e., their contribution to chemical potential is proportional to the temperature. Typically, the value of the proportionality constant is adjusted by fit to the available experimental data.
The entropy difference could not be obtained at zero K as the vapor-phase entropy has singularity there (for discussion of the singularity please refer to Appendix A). Therefore, the separate calculation path was proposed for the total difference in chemical potential at the condensed/vapor, i.e., liquid (or solid)/vapor (l-v) phase transition for nitrogen. The enthalpy difference at the solid-vapor transition is calculated at zero K. The entropy difference at this transition is obtained at normal conditions. The difference in total chemical potential must be supplemented by changes occurring during the transition to these points. The entire paths are fairly complex, nevertheless, the same results were obtained independently by Jackson and Walsh [37] for bulk and by us for the properties of gas adsorbed at the surfaces of solids [38]. The equations here are written in formulation given in reference [38], with the sign reversed, i.e., using dissolution energy as defined in Section 2.1, for the vapor-solid chemical potential difference: where the normal temperature and pressure are defined as: T 0 = 20°C = 298.15 K and p 0 = 1 bar. The terms in the above equation are defined as: (i) The first term is the enthalpy change at dissolution for single N atom, calculated as ab initio energy difference between the vapor and the solid (liquid): It was shown recently [39] that the energy difference obtained from DFT calculations (Equation (1)) does not correspond to the enthalpy change during vaporization as the thermodynamic state energy is increased by the energy of vibrations at the ground state, the effect called zero-point energy (ZPE) [39]. As these values are not identical, they contribute to the enthalpy difference giving where the zero-point energy difference is: As pointed out by Ponomareva et al. [34] and discussed in Section 2.1, the energy depends on the atomic configuration used for the calculation. Therefore, in our approach, the enthalpy difference was evaluated for configurations resulting from the MD simulations at 500 K. As discussed in Section 2.1, the dissolution energy was obtained using its MD value, the average thermal energy of single N atom at 500 K ∆E T (1 N) and zero-point energy E ZPE l (N-Ga) (see Equation (4)). (ii) Difference in chemical potential related to the entropy change at dissolution, calculated at normal conditions, (iii) Difference in enthalpy change between both phases at the transition from 0 K to normal conditions, i.e., T 0 , (iv) The chemical potential change caused by the temperature change from T 0 to , (v) The pressure induced change of the chemical potential as Term related to the concentration of N in the metal liquid. The latter contribution may be calculated using standard chemical approximations for ideal solutions: where x is the concentration of N atoms in the solution.

Dissolution of Molecular Nitrogen in Liquid Gallium
The enthalpy of dissolution ((i) in the previous Section) for a single N atom consists of two contributions: energy of dissolution ∆E DFT dis (Ga − 1 2 N 2 ) = = 0.0712 eV. Therefore, zero-point energy difference is ∆E ZPE lv (Ga − 1 2 N 2 ) = 0.143 eV. In summary, the enthalpy of dissolution for single N atom in the gallium metal, from N 2 molecule is: ∆H dis (Ga − 1 2 N 2 ) = 0.696 eV/atom.
Please note that the contribution of dissolution entropy to the free energy (ii) ∆G S−vap (Ga − 1 2 N 2 ) could be easily obtained. As both the vapor and the liquid Ncontaining phases are not ordered, the entropy change could be expressed by the appropriate volume ratio: From the ideal gas law, the volume associated with single N atom at normal temperature and pressure is: v v (N) = v v (N 2 ) 2 = 1.861 × 10 4 Å 3 . The volume associated with incorporation of the N atom into the Ga cluster is v l (N) = 2.990 Å 3 . Thus, the dissolution entropy related free energy change described by Equation (16) is: ∆G S−vap (Ga − 1 2 N 2 ) = 0.224 eV. For the vapor phase, the entropy terms could be evaluated quite easily. In a good approximation, the specific heat of nitrogen in the gas phase is determined by the equipartition principle, i.e., it is constant. The thermal contributions (iii) and (iv) could be obtained be separate calculation for both phases. The vapor-phase contribution is obtained directly, as the heat capacity is constant.
The enthalpy difference is ∆H therm (N 2 ) = H v (T 0 ) − H v (0) = 8.67 kJ/mole = 8.986 × 10 −2 eV/molecule. The latter value for single N 2 molecule could be recalculated for single N atom as: .493 × 10 −2 eV/atom. Generally, the contribution ∆H therm (N 2 ) is small in relation to the enthalpy change at dissolution or to the difference in chemical potential related to dissolution entropy change thus not affecting the entire result considerably.
The second term is the nitrogen free energy difference between standard T 0 and a selected temperature T. This could be directly obtained from already published data concerning the chemical potential of gaseous nitrogen at normal pressure [40]: where y ≡ T 1000 is scaled temperature. From this data the free energy difference is readily obtained as: The entropy contribution in the liquid may be obtained from evaluation of phonon related effects, directly. The phonon part of the vibrational energy E vib , the specific heat C vib , entropy S vib , and free energy F vib associated with the nitrogen atom immersed in gallium may be obtained as [39,41]: where x j ≡ ω j k B T and ω j is phonon frequency of a j-th phonon mode. A simplified treatment may be used, based on Debye theory in which the acoustic phonon frequencies are approximated by linear dependence. Accordingly, the maximum phonon energy, known as the Debye energy, and its equivalents: frequency and the temperature, are related as follows: Using this simplified representation, the spectra sums are replaced by integrals, giving: The N-Ga Debye temperature was obtained from the fit of ab initio derived specific heat to the temperature dependence given by Equations (23)- (26). The obtained value equal θ D = 552 K was used to evaluate other thermal contributions given by Equations (23)- (26).
Evaluating the pressure related term ((v) in the previous Section) one must account that in Equation (17), the pressure of molecular nitrogen is used. In our model, the expression scaled to atomic nitrogen has to be used, i.e., ∆G pres (N) = 1 2 ∆G pres (N 2 ) = − k B T 2 ln p p 0 , where p 0 = 1 bar.
At equilibrium, the chemical potential difference vanishes, i.e., ∆µ vl (p, T) = 0 which allows us to express the equilibrium pressure above specified liquid phase via chemical potential contribution in the ideal gas approximation: In the real gas case (at high pressure), the pressure p in Equation (27) must be interpreted as fugacity of N 2 gas (or activity, if divided by p 0 ) accounting for intermolecular interactions due to high density of the gas.
Therefore, the dependence of the concentration of nitrogen in liquid Ga on activity of the N 2 gas over the solution could be obtained. Such a dependence for several temperatures is presented in Figure 11.
As in expressions for thermodynamic functions, the considered activities of N 2 gas replace the pressure where ideal gas equation of state is not fulfilled, the corresponding pressure had to be evaluated using the modified equation of state. Such equation of state of nitrogen derived from experimental data by Jacobsen et al. [42] and confirmed by molecular dynamics and ab initio simulations by Strak et al. [43,44] was used for the activity to pressure transition. In such a way, the nitrogen solubility as a function of pressure of N 2 pressure was evaluated and is presented in Figure 12.
From the presented calculations, it follows that the solubility of nitrogen in liquid Ga is low, and accordingly, high pressures are required to attain technically viable rates of GaN crystallization [11,35,36]. The technically used N 2 pressure was recovered using nitrogen equation of state for high pressures and high temperatures. It is interesting to mention that thermodynamic activity of a compressed N 2 gas at technically relevant pressures of the order of 10 kbar, is nominally much higher than the corresponding pressure itself. For example at 10 kbar of the molecular nitrogen pressure, the N 2 activity is about a N 2 ∼ = 141, 900 at T = 1300 K, i.e., about 14 times higher than the ideal gas value. The solubility data [34] experimentally evaluated for conditions corresponding to p-T coordinates of Ga-GaN-N 2 triple points are in reasonably good agreement with the presented ab initio results. The result is impressive as the obtained concentrations are exponentially dependent on the calculated quantities, i.e., thermodynamic potentials. The agreement is much better for lower concentrations where the ideal solution approximation is more precise. For higher pressures, the calculated data are consistently smaller than the measured ones, which is attributed to deviation from ideal solution approximation.
The calculated concentrations are relatively small, in the range of several atomic promilles, as already determined experimentally. They limit the possible crystallization rates as they are proportional to the equilibrium concentration and the supersaturation.
The latter could be established using temperature difference, i.e., the difference in solubility, normalized to the equilibrium concentration. Thus, the temperature dependence is of importance here. As it is shown in Figure 12, the concentration changes rapidly with the temperature, allowing the obtaining of high supersaturation via temperature difference between the dissolution and the growth zone in the metal. Unfortunately, high supersaturation accelerated growth often leads to pronounced Mullins-Sekerka instability [45,46] and consequently deterioration of the quality of the resulting crystals [47].
From this point it is interesting to analyze the temperature dependence of nitrogen solubility at constant activity, shown in Figure 13. These data could be used in the plot of the Van't Hoff type following the Van't Hoff relation: to determine thermodynamic heat of dissolution ∆H dis therm . From the linear fit, the value: ∆H dis therm = 0.848 ± 0.001 eV/atom was obtained what is far from the ab initio obtained value ∆H dis (Ga − 1 2 N 2 ) = 0.696 eV/atom. As it is shown, the plot deviates from linear regime significantly. Such dependence is not typical for enthalpy dominated phase transition as the entropy related change of the enthalpy of dissolution is about 21% for the temperature increase to 1800 K. That confirms importance of the entropy related terms in the dissolution of nitrogen in liquid Ga, indicating that the above combined energy-entropy approach is necessary to obtain good agreement with the experimental data. = 0.0712 eV, i.e., equal to ∆E ZPE lv (Fe − 1 2 N 2 ) = 0.139 eV. Altogether, the enthalpy of dissolution is ∆H dis (Fe − 1 2 N 2 ) = 0.172 eV. The entropy of dissolution may be obtained from ab initio data. The volume associated with N atom in the Fe cluster is v l (N) = 2.200 Å 3 . Using the ideal gas volume for normal pressure and temperature v v (N) = v v (N 2 ) 2 = 1.861 × 10 4 Å 3 the dissolution free energy change could be obtained from Equation (5) to get ∆G S−vap (Fe − 1 2 N 2 ) = 0.232 eV. Thus, this value is only slightly different from that obtained for liquid Ga.

Dissolution of Molecular Nitrogen in Metallic Iron
The entropy terms could be obtained following these for the previous case. The enthalpy difference for vapor is identical and equal to: .251 × 10 −2 eV/atom. The free energy difference of gaseous nitrogen is given by Equation (7). The remaining solid contributions may be obtained from Equation (16), using the Debye temperature which for N dissolved in metallic Fe was equal θ D = 512 K.
The pressure dependence of the N concentration in the Fe solvent could be obtained. Such data are presented in Figure 14. Please note that in the presented diagrams the pressure related coordinate represents the chemical potential of compressed N 2 gas, i.e., the activity. The explicit pressure dependence should be recovered expressing the pressure as a function of activity via nitrogen equation of state. From the obtained data it follows that the concentration values are more than one order of magnitude higher than in case of gallium. Thus, the Fe-based solution seems to be a much better option for efficient crystallization of large GaN single crystals. More physical insight may be obtained using technically relevant variable i.e., the pressure, which is recalculated from activity using the equation of state [43][44][45]. The corresponding results are presented in Figure 15. The experimental solubility data obtained according to the procedure described in Section 2.2, for a constant temperature of 1708 K, are also included into the diagram. The experiment clearly confirmed high concentration of N in the Fe metal. Also, the character of the isothermal pressure dependence is well reproduced. The obtained theoretical solubility is approximately half of the measured values. This is most likely related to the two factors: first mentioned in Ga case, is related to nonideality of the solution. The effective interaction is much stronger in the case of Fe than in the case of Ga. Therefore the interaction related deviation from ideality is larger. Second factor is related to the enthalpy calculations. In the above ab initio modeling only single N atom was inserted in the Fe cluster. In fact in the concentration range of 10 at.%, additional nitrogen atoms may be present in the neighborhood. The N-N attraction is strong, which leads to the decrease of the system energy and consequently lowering its chemical potential thus increasing the nitrogen concentration as observed experimentally.
The theoretical evaluations indicate that the solubility of nitrogen at technically possible limit conditions for large volume high temperature gas reactors, i.e., 10 kbar attains considerable values of about 5.8 at.% at T = 1800 K. These data are only slightly reduced by the lower temperatures, i.e., about 3.0 at.% at T = 1300 K. Thus, the temperature increase has drastically smaller influence on the concentration of N in iron than in gallium. This can complicate establishing supersaturation in the growth solution necessary for crystallization by application of the temperature gradient.
The temperature dependence may be studied using concentration vs inverse temperature plot for selected values of nitrogen activity as shown in Figure 16. The Van't Hoff relation (Equation (28)) may be again used to determine thermodynamic heat of dissolution ∆H dis therm . From the linear fit the following value was obtained ∆H dis therm = 0.2889 ± 0.001 eV. This is much smaller than the ab initio obtained value ∆H dis (Fe − 1 2 N 2 ) = 0.172 eV that indicates that the entropy related, temperature-dependent term drastically lowers the enthalpy of dissolution. Additionally, the solubility does not follow the linear plot of the Van't Hoff type, similarly to gallium. Even though the interaction term is approximately four times smaller than for gallium, the deviation from linearity is similar that confirms the role of entropy related terms in the dissolution of nitrogen in liquid metals.
The observed temperature dependence poses more stringent requirement for design of crystal growth apparatus in case of Fe solvent. In order to obtain a comparable supersaturation, much larger temperature difference than for gallium, should be applied. This could be beneficial as the best conditions for the growth is small supersaturation. In perspective, large high-quality crystals of GaN could be obtained from solutions in liquid Fe. tropy difference is calculated via normal conditions while enthalpy difference is calculated via zero K. The temperature-dependent contribution was obtained from Debye theory, for which the Debye temperature for Ga was θ D = 552 K and for Fe θ D = 512 K, i.e., similar values. The Debye contributions are large and temperature-dependent (Equation (13)). For N in liquid Ga, they are ∆G S-therm = −0.488 eV and ∆G S-therm = −1.039 eV for T = 1000 K and T = 2000 K, respectively. For N-in-Fe, they are ∆G S-therm = −0.475 eV and ∆G S-therm = −1.006 eV for T = 1000 K and T = 2000 K , respectively. Thus, they are significant and similar for both metals, seriously affecting the equilibrium. The other contributions, such as thermal enthalpy difference (Equation (12)) or concentration dependent are much smaller. The obtained equilibria for the N 2 -Ga and N 2 -Fe systems are much different. The N solubility was obtained for technically amenable part, i.e., for temperatures between 1300 K and 1800 K. The N 2 activity was limited to a(N 2 ) ≤ 2 × 10 5 (corresponding to the N 2 pressure of 10 4 bar at 1300 K). For these conditions, the solubility of nitrogen in liquid Ga attains relatively low values limited by C N < 3 × 10 −3 at. fr. For Fe these concentrations are much higher C N ≤ 9 × 10 −2 at. fr. The theoretical evaluations were clearly supported by available experimental results for both systems.
The solubility is temperature-dependent, as demonstrated by the relatively high enthalpy of dissolution, obtained from the concentration dependence which was ∆H dis therm = 0.848 ± 0.001 eV and ∆H dis therm = 0.2889 ± 0.001 eV for Ga and Fe, respectively. These values are different from the ab initio values which were ∆H dis (Ga − 1 2 N 2 ) = 0.696 eV and ∆H dis (Fe − 1 2 N 2 ) = 0.172 eV, respectively. The values following from the concentration dependence and the ab initio ones are considerably different, which confirms the important role of entropy contributions. It could be noted that much higher values of the enthalpy are observed for liquid Ga, which indicates the steeper temperature induced change of the N concentration for Ga. For Fe, the temperature dependence of the N concentration is much weaker what can be relevant for a proper design of the crystallization experiment where low supersaturation in this otherwise promising, growth solution has to be created. Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.
Sample Availability: Samples of the compound Fe with dissolved N are available from the authors.

Appendix A. Singularity of Entropy Contribution to Solid-Vapor Transition at Zero K
Equilibrium between solid and vapor and entail chemical potential equality of both phases: µ v (p, T) = µ s (p, T) In general, the equality holds for any species exchanged between phases. Here for simplicity, we assume that single species is exchanged, thus for simplicity any reference is neglected. The chemical potential can be divided into enthalpy and entropy contributions, for both phases µ v,s (p, T) = h v,s (p, T) − Ts v,s (p, T) (A2) In ideal gas approximation, the pressure change of the chemical potential may be approximated as: For solid phase, the specific volume is 4 orders of magnitude smaller, thus the solid phase contribution is neglected. Then the equilibrium pressure over the solid may be obtained as Van't Hoff relation: where ∆h vs = h v − h s > 0 and ∆s vs = s v − s s > 0, are positively defined enthalpy and entropy of vaporization at normal conditions, respectively. Generally, these thermodynamic functions are temperature-dependent. In ideal gas approximation the specific heat C p is constant vs temperature, independent of the density (therefore a pressure could be used).
The temperature dependence of the enthalpy of the vapor as: For the solid phase, the specific heat is given by Debye which, for low temperatures, may be approximated by cubic dependence C s (T) ∼ αT 3 which gives the following dependence of the enthalpy of solids at low temperatures In summary, the enthalpies of both phases are regular at T = 0 K and application of zero K reference state for enthalpy is correct. Therefore, zero temperature is selected as the reference state, i.e., T 0 = 0 K.
Quite different conclusions are obtained for the entropy term. For the solid, the entropy contribution is regular at T = 0 K: The vapor entropy term is singular, with logarithmic divergence at zero K: Thus, the entropy at zero K cannot be used in determination of vapor-solid equilibria. Therefore, for entropy, normal pressure and temperature is selected: p 0 = 1 bar, T 0 = 273.15 K.
In general, the latter is a manifestation of well-known paradox, the entropy at zero K should be zero. For the ideal gas, the entropy is divergent, thus the vapor phase could not exist at zero K. It vanishes as the pressure of the vapor vanishes. In order to find the proper dependence of the vapor pressure at zero K the Van't Hoff relation could be used (Equation (A4)). The precise relations is given by Equation (2), but for analysis of low temperature regime, the relation given by Equation (A4) is sufficient. Adopting above selection of zero temperature reference point for enthalpy, the enthalpy change is: ∆h vs = ∆H vap (0) + C p T − αT 4 /4 (A9)