Ab Initio Molecular Dynamics Insight to Structural Phase Transition and Thermal Decomposition of InN

1. Introduction

Indium Nitride (InN) is an intriguing and very important although relatively less known, semiconductor. Its importance is related to revolutionary development of InGaN based light sources – Light Emitting Diodes and Laser Diodes [1], awarded by the Nobel Prize in physics’2014 [2]. Until 2002 InN has been considered as a wide bandgap (about 2 eV) member of the III-N semiconductor family including GaN, AlN and BN. Due to significant improvements in high purity molecular beam epitaxy (MBE) of InN, it was convincingly demonstrated that the bandgap of high structural quality pure InN is as low 0.64eV [3,4]. Therefore the energy spectrum covered by III-N semiconductors: E_g(InN) = 0.64 eV, E_g(GaN) = 3.47 eV [5,6], E_g(AlN) = 6.08 eV [7] and E_g(wBN) = 6.39 eV [8] became even more interesting extending from far infrared to ultraviolet spectral range. Moreover, the electron effective mass m^*(InN) = 0.07m₀ [4], corresponding to the “new” energy gap is much smaller than that of the other nitrides: GaN - m^*(GaN) = 0.20m₀ or AlN - m^*(AlN)=0.32m₀ [4]. It indicates a possibility of hosting very high mobility electron gas in the InN crystal. However, the growth of InN single crystals is still a challenge due to specific thermodynamic properties of this nitride.

InN crystal is relatively soft despite it contains atomic nitrogen. The InN bonding is less covalent and more ionic as well as much weaker than in the other nitrides what is reflected in larger lattice constants and much smaller bandgap in InN [4]. Therefore, as you would expect, the thermodynamic properties of the In-N system are affected by both strong triple bonding in N₂ molecule lowering the free energy of the In+N₂ system, and relatively weak bonding of the InN crystal [9,10]. That leads to very high pressures of nitrogen at relatively low temperatures along the InN(s)-In(l)-N₂(g) triple phase stability line [11]. In addition, the solubility of nitrogen in liquid indium is extremely low so the growth of InN bulk crystal from indium solution is difficult [11,12]. As a consequence, large bulk crystals of InN have not been grown so far and due to the lack of high quality reference material, experimental investigations of exciting physical properties of InN are difficult.

An interesting route facilitating the efforts towards crystallization of InN and also providing valuable insight into the physical properties of the InN(s)-In(l)-N₂(g) system could be obtained from molecular dynamics (MD) ab initio simulations. In this study, we introduce the ab initio MD simulations as a potential tool to evaluate phase behavior of InN in wide pressure and temperature range. It can be useful for further exploration of phase diagram of InN: its experimental verification and its use as a background for crystal growth. An additional, encouraging motivation was a very reasonable agreement of the proposed approach with the experimental data we achieved recently for GaN [13].

1.1. Phase Diagram of InN

Figure 1 summarizes the most consistent experimental data on the InN p-T phase diagram obtained by Saitoh et al. [14] and Soignard et al. [15] based on the in situ X-ray diffraction (XRD) measurements at high pressure – high temperature conditions.

The diagram shown in Figure 1 indicates the following behaviour of InN:

• At the lower pressure range (up to 7-10 GPa) InN decomposes when heated above approx. 710 °C [8] - green arrow in Figure 1. It was also confirmed by our earlier differential thermal analysis (DTA) experiments up to 2 GPa [12] at high pressure of N₂ gas. The independence of decomposition temperature on pressure suggests that the equilibrium temperatures in the InN-In-N₂ system can be even lower than the measured ones.
• At 12 GPa (low T) to 7-10 GPa (high T), a structural phase transition from hexagonal wurtzite to cubic rocksalt phase induced by increasing pressure is observed - blue arrow in Figure 1. The borderline between the two solid phases is inclined towards lower pressures however different linear or non-linear character of this line follows from XRD experiments reported in [14] and [15], respectively. A possible verification could be checking if at i.e. 8 GPa the InN crystal in its wurtzite phase transforms into rocksalt at heating to 800–1000 K – red arrow in Figure 1.
• The InN crystal in its high pressure rocksalt phase also decomposes at high temperatures but the decomposition temperature strongly increases with increasing pressure – violet arrow in Figure 1.
• The decomposition is suppressed only at pressure as high as > 16 GPa and then the congruent melting of InN (without formation of N₂ thus reversible) is possible – magenta arrow in Figure 1.

In the next sections we explore by the ab initio MD simulations the following phase transitions of InN:

• wurtzite-to-rocksalt structural phase transition induced by high pressure (blue arrow);
• wurtzite-to-rocksalt structural phase transition induced by high temperature at 8 GPa (red arrow);
• Thermal decomposition of InN crystal in both wurtzite and rocksalt phases (green and violet arrow, respectively).

In this paper, special attention is given to the determination of the p-T stability range of InN in its wurtzite modification, which is particularly important for the design of future crystallisation experiments in high-pressure reactors with high volume routinely used for diamond or BN growth. Congruent melting of InN (magenta arrow in Figure 1) was not considered.

2. The Simulation Method

Temporal evolution by the ab initio MD simulations requires a huge amount of calculations which is challenging because at each time step, the solutions of full quantum mechanical problem via density functional theory (DFT) Kohn-Sham equations have to be obtained [16]. The forces acting on each atom in the system, are obtained via Hellmann-Feynman theorem used for time integration of MD time evolution equations [17]. The standard time step, 0.2 fs long, was employed in this study. Total simulation time lasted 30 000 steps, that was carried out for each p,T set. For the time dependent MD calculations the optimal choice is Verlet integration algorithm [17] that is precise to fourth order in the time step length. This algorithmically simple method requires a single evaluation of the force acting on each atom at every time step thus it is the most effective numerically. The Verlet algorithm is implemented in the time dependent simulations in SIESTA (Spanish Initiative for Electronic Simulations of Thousands of Atoms) DFT package selected for our simulations [18]. The solution of nonlinear Kohn-Sham equations is expressed as sum of finite radius pseudo-atomic orbital functions set [19] that is reduced in number by application of the Troullier-Martins norm-conserving nitrogen and indium pseudopotentials in the Kleinmann-Bylander formulation [20,21]. The exchange-correlation energy was approximated by Perdew-Burke-Ernzerhof functional developed for solid and surface calculations (PBEsol) [22]. The solutions of the Kohn-Sham nonlinear equations were obtained by self-consistent field (SCF) loop, terminated when the maximum difference between two consecutive values of all elements of the density matrix fell below 10^-4.

In the reported simulations, the finite size InN 3×3×2 supercell was used containing total 72 indium and nitrogen atoms, arranged in the hexagonal wurtzite lattice. In some cases intentionally, the initial arrangement of In-N atoms in the ideal cubic rocksalt lattice stable at high pressures [14,15], was used. Periodic boundary conditions were imposed for atom motion, pseudopotentials and electric potential. For the considered systems, the Brillouin zone k-integration was replaced by sum over 3×3×3 Monkhorst-Pack k-point mesh. The initial supercells used in our study are presented in Figure 4a and b.

The parameterization of the equations and solution basis was verified by direct comparison to the experimental data for wurtzite InN, gaseous N₂ and pure In. The ab initio lattice parameters of wurtzite InN are $a_{InN}^{DFT}=3.533$ Å and $c_{InN}^{DFT}=5.693$ Å that remain in a very good agreement with the x-ray data: $a_{InN}^{exp.}=3.5374\left(1\right)$ Å and $c_{InN}^{exp.}=5.7027\left(5\right)$ Å [23], respectively.

The calculated dissociation energy and bond length of the N₂ molecule were $\Delta E_{diss.}^{DFT}\left({\mathrm{N}}_2\right)=9.801$ eV and $d_{\mathrm{N}-\mathrm{N}}^{DFT}=1.092$ Å, in a good agreement with the experimental $\Delta E_{diss.}^{exp.}\left({\mathrm{N}}_2\right)=9.790$ eV [24] and $d_{\mathrm{N}-\mathrm{N}}^{exp.}=1.097$ Å [25], respectively.

The properties of metallic indium were evaluated including density and cohesive energy. The ab initio energy of In atomization was $\Delta E_{DFT}^{atom}$(In) = 2.39 eV/atom, close to the experimental value of vaporization enthalpy $\Delta E_{exp.}$ = 2.35 eV/atom [26]. The density of the In liquid at the melting temperature T_m(In) = 156.6 °C, is ${\rho }_{exp.}$(In) = 7.02 g/cm³ which is comparable with the ab initio simulation result: ${\rho }_{DFT}$(In) = 7.05 g/cm³.

For the time-integration procedure in the imposed NVT conditions i.e. corresponding to the Canonical Ensemble, the Nose-Hoover thermostat for temperature stabilization [27,28] was employed. From the initial setup of atoms in the form of the wurtzite/rocksalt InN lattice (Figure 1), the target temperature was achieved by rescaling the velocities of all atoms to the kinetic energy obtained from equipartition principle for a selected temperature T. The molecular pressure was determined using Irving-Kirkwood formula [29]. More details related to the ab initio MD simulations of nitrides may be found in Ref. [13]. The pressure – temperature conditions used in this study are included into Figure 1.

3. Result and Discussion

3.1. Pressure Induced Solid-Solid Phase Transition at Low Temperature

A sequence of the InN clusters annealed at 800 K for a period of 6 ps, at different pressure, is shown in Figure 2. The sequence illustrates the wurtzite – rocksalt structural phase transition starting at pressure of 8 GPa. The structure identification is based on the analysis of coordination numbers in the cluster presented in Figure 2c in which both 4 fold (wurtzite) and 6 fold (rocksalt) coordination spheres were found. The transition pressure of 8 GPa is lower than the one following from the high pressure XRD experiments (10-11GPa) [14,15]. The possible causes for this disagreement are hidden in both theoretical and experimental approaches. The mixed coordination observed at the transition point (8 GPa, 800 K) suggests a “nucleation – growth” mechanism of the structural transformation. Therefore a hysteresis is expected what is indeed observed in the experiments making the relation between theoretical and experimental transition pressure values, consistent. Moreover, intrinsic errors in both the simulations and the measurements have to be taken into account.

According to the simulation, at 12 GPa the crystal is fully converted into its high-pressure rocksalt phase – Figure 2d, as confirmed by the coordination analysis, which is in agreement with the experimental data.

Since in the simulation the NVT Canonical Ensemble has been used where pressure was a function of the NVT set, and remembering that the pressure induced wurtzite – rocksalt phase transition includes a significant change of the InN volume [14], the pressure trajectory during the simulation has to be analyzed. The analysis for the wurtzite – rocksalt transition during 6ps indicated an averaged pressure decrease not bigger than 300MPa at intrinsic pressure value scattering of $\pm$ 500 MPa.

3.2. Temperature Induced Solid-Solid Phase Transition at 8 GPa

In Figure 3 we have shown a sequence of clusters corresponding to 8 GPa and different temperatures ranging from 800 to 2000 K. The initial supercell was wurtzite, Figure 4a.

The sequence starts from the cluster with mixed coordination already presented in Figure 3c as a border case between wurtzite and rocksalt phases at 800 K. With increasing temperature (1100 K and 1400 K), an ordering towards cubic rocksalt phase is clearly visible which confirms that the borderline between the two solid phases is negatively inclined and rather non-linear like it was suggested in [15].

For the highest temperature of 2000 K, the crystalline order is lost and the crystal decomposes with the formation of N₂ molecules (decomposition process is analyzed in the next sub-section).

Given that at 8 GPa and 1100 K the phase of InN clearly evolves towards rocksalt, we decided it would be more consistent for 1400 K and 2000 K, to start the simulations from the initial supercell in the rocksalt structure, Figure 4b. The corresponding comparison for 1400 K is shown in Figure 4.

Figure 4. Effect of 6 ps virtual annealing at 8 GPa and 1400 K on the InN cluster: c – with initial wurtzite supercell (a), d – with initial rocksalt supercell (b), green balls – indium atoms, blue balls – nitrogen atoms, red balls - nitrogen atoms forming N₂ molecules.

Figure 4. Effect of 6 ps virtual annealing at 8 GPa and 1400 K on the InN cluster: c – with initial wurtzite supercell (a), d – with initial rocksalt supercell (b), green balls – indium atoms, blue balls – nitrogen atoms, red balls - nitrogen atoms forming N₂ molecules.

The situations are quite different: for the initial supercell in wurtzite, its conversion to the rocksalt structure with no signs of lattice disintegration is visible whereas for the initial supercell in rocksalt, the lattice is thermally destabilized and moreover, N₂ molecules start to emerge what is in agreement with the experimental observations. It indicates that for the applied simulation algorithm and time scale it is reasonable to start the virtual heating of the crystal avoiding structural phase transition during thermal evolution of the system.

3.3. Decomposition of InN

The experimental diagram of InN phase stability (metastability?) shown in Figure 1, indicates that the crystal should decompose at high temperatures to indium and N₂ at pressure up to approx. 16 GPa. In the XRD experiments [14] the thermal decomposition was reflected by irreversible disappearance of the XRD spectral features of the InN crystalline phases. In contrast, at pressure of 19.2 GPa the XRD spectrum corresponding to the rocksalt InN crystal lattice also disappeared however at cooling, its recovery has been observed. Therefore it was suggested that there was a sign of congruent melting of InN and that the congruent melting of InN is possible only in for the rocksalt high pressure phase due to the wurtzite InN decomposition before reaching its melting point in whole pressure range of the wurtzite InN stability.

The ab initio MD simulations reproduced the thermal decomposition of InN in terms of the crystal lattice destruction and formation of N₂ molecules. Figure 5abcd including the results for 2000 K demonstrates formation of the N₂ dimers for all considered pressure conditions: 0 GPa, 6 GPa, 8 GPa and 12 GPa. At lower temperature, T = 1400 K, the decomposition with N₂ formation was observed for 0 GPa, 6 GPa and 8 GPa but for 12 GPa, the formation of N₂ molecules was suppressed (Figure 5h) and the rocksalt structure was preserved. It suggests that at T = 1400 K, the InN crystal (rocksalt) is stable at p = 12 GPa and unstable at p = 8 GPa what is in a very good agreement with experimental suggestions illustrated in Figure 1.

The formation of N₂ molecules was identified by the analysis of distribution of interatomic distances between nitrogen atoms (Figure 6) and further confirmed by the Crystal Orbital Hamilton Population (COHP) analysis of the selected pairs of nitrogen atoms [30,31].

4. Summary and Outlook

The reported ab initio MD approach simulates both pressure and temperature induced structural phase transitions of InN. The results show that the wurtzite-rocksalt borderline is negatively inclined and nonlinear like it follows from [15] but is not supported by [14]. The experimental verification of this controversy is strongly needed for precise evaluation of the wurtzite InN stability range. Then the optimum temperature and corresponding pressure can be found and explored for crystallization of the nitride.

Interestingly, thermal decomposition of InN with formation of N₂ molecules was identified in the simulation what is a great advantage of the ab initio MD. The approach could be especially useful for evaluation of the still unknown melting curve of InN for which suppression of the N₂ formation is expected. For high temperatures ($\ge$ 1400 K) expected for melting, the relatively short simulation time of 6 ps is sufficient to complete the relevant transitions. For the InN melting experiments, the laser heated diamond anvil cell technique coupled with XRD, EXAFS or Raman spectroscopy i.e. [32] will be used.

J.P.—ab initio calculations; S.K.—thermodynamics modeling, supervising theory, manuscript preparation; B.S.—concept, experiment planning and supervising, manuscript preparation; P.S.—experimental N-in-Fe solubility evaluation; S.P.—data analysis, critical evaluation of the manuscript; I.G.—concept, data analysis, manuscript preparation. All authors have read and agreed to the published version of the manuscript.