Self-assembly of Carbon Vacancies in Sub-stoichiometric ZrC1−x

Sub-stoichiometric interstitial compounds, including binary transition metal carbides (MC1−x), maintain structural stability even if they accommodate abundant anion vacancies. This unique character endows them with variable-composition, diverse-configuration and controllable-performance through composition and structure design. Herein, the evolution of carbon vacancy (VC) configuration in sub-stoichiometric ZrC1−x is investigated by combining the cluster expansion method and first-principles calculations. We report the interesting self-assembly of VCs and the fingerprint VC configuration (VC triplet constructed by 3rd nearest neighboring vacancies) in all the low energy structures of ZrC1−x. When VC concentration is higher than the critical value of 0.5 (x > 0.5), the 2nd nearest neighboring VC configurations with strongly repulsive interaction inevitably appear, and meanwhile, the system energy (or formation enthalpy) of ZrC1−x increases sharply which suggests the material may lose phase stability. The present results clarify why ZrC1−x bears a huge amount of VCs, tends towards VC ordering, and retains stability up to a stoichiometry of x = 0.5.

In this paper, the energetics of ZrC 1−x at various V C concentrations is studied using the state-of-the-art first principles and cluster expansion method. We report a fingerprint structural unit, namely the V C triplet constructed by 3 rd nearest neighbor carbon vacancies in all the predicted low energy structures of ZrC 1−x . The defective structure would lose stability when V C concentration reaches over 0.5 and simultaneously, highly repulsive V C s with the 2 nd nearest neighbor coordination unavoidably appear. The results explain some longstanding puzzles in non-stoichiometric interstitial compounds, and this study may also shed lights on how to design or tailor the performances of promising transition carbides.
The energetics of sub-stoichiometric ZrC 1−x presents more information for vacancy tolerance and ordering capability. Firstly, ZrC 1−x displays significant tolerance to high concentration of V C s. The mixing enthalpy of ZrC 1−x with random V C distribution is shown by the dashed curve in Fig. 1. The predicted mixing enthalpies retain negative in the composition range of 0 < x < 0.59. This result suggests that sub-stoichiometric ZrC 1−x phases with a huge amount of vacancies are energetically favorable. Otherwise, sub-stoichiometric ZrC 1−x with positive mixing enthalpy would spontaneously decompose into f.c.c. Zr and B1-ZrC competition phases. Furthermore, sub-stoichiometric ZrC 1−x demonstrates obvious tendency of V C ordering because numerous V C configurations have lower mixing enthalpies than the disordered V C distribution in Fig. 1 Fig. 1, are always negative. Meanwhile, with increase of V C concentration, the ordering enthalpy continuously decreases and reaches the minimum value at around x ~ 0.5. This suggests the fact that the higher the V C concentration, the stronger the ordering tendency of V C s, and in addition, Zr 2 C has the most obviously ordering tendency in all studied sub-stoichiometric ZrC 1−x structures. We expect that ordered Zr 2 C should be the most possible ordered phase synthesized in experiments. This result is consistent with the discovery of ordered Zr 2 C in experiments 7,18 . Besides, the calculated ordering enthalpies of Zr 8 C 7 and Zr 6 C 5 are -30 meV/cation and -55 meV/cation, respectively, which are comparable to that of V 8 C 7 and V 6 C 5 ordered phases (around -20 meV/cation) 3 . Because all the GS structures have negative ordering enthalpies, other predicted GS phases would be fabricated by careful controlling of experimental conditions. Vacancy configurations. Besides the predicted ground states, the structural characteristics of low energy configurations between neighboring GSs are also important to depict the evolution features of V C s in sub-stoichiomitric ZrC 1−x . Figure 2 illustrates the radius distribution function of V C pairs in selected low energy structures with various compositions. V C configurations are displayed by the presence or absence of certain neighboring V C pairs, including 1NN (1 st nearest neighboring), 2NN (2 nd nearest neighboring), 3NN (3 rd nearest neighboring), and 4NN (4 th nearest neighboring) V C pairs. It is striking to notice that 3NN V C pair is found in all structures, especially only 3NN V C pair appears in the GS structures of Zr 8 C 7 and Zr 6 C 5 which have low V C concentrations. When V C concentration increases, 1NN V C pair subsequently presents in the GS of Zr 4 C 3 , then 4NN V C pair emerges in the GSs of Zr 3 C 2 and Zr 2 C. Besides, 1NN V C pair is also found in the low energy structures near the GSs, such as Zr 32 C 27 and Zr 32 C 26 . At higher V C concentration, 4NN V C pair is identified in low energy structures near the GSs, such as Zr 32 C 23 , Zr 32 C 18 , and Zr 32 C 17 . When V C concentration is higher than 50%, 2NN V C pair inevitably arises, such as Zr 32 C 15 and Zr 32 C 14 . Accordingly, their formation energies increase sharply with the presence of 2NN V C coordination. Figure 2 clearly displays the interesting characteristics of V C configurations: most frequently observed 3NN configuration, available 1NN and 4NN configurations, as well as the unfavorable or forbidden 2NN configuration. It is important to disclose the global evolution characteristics of V C s in binary transition metal carbides. Indeed, it was both experimentally and theoretically found that the vacancies in nonstoichiometric interstitial carbides preferred the 3NN shell and excluded the 2NN shells 3,22,24 , which was consistent with our calculation results. Hart et al. investigated the vacancy distribution in low energy structures of TiC 1−x using cluster expansion method. Although they did not mention the evolution of V C configurations at various concentrations, carbon vacancies were found arranging themselves in [112] rows 25 , which was along the crystal direction of 3NN V C pairs, and 2NN configurations gradually disappeared during annealing or ordering.
We expect that there would be inherent local V C configurations commonly appearing in the low energy structures because their mixing enthalpies are so close to the GS envelope line. In ZrC 1−x with stoichiometry near ordered GS phases (around compositions of Zr 2 C, Zr 4 C 3 , Zr 3 C 2 etc.), it would be unusual for any abrupt change of V C patterns. After careful analysis of V C configurations, the common configuration unit, i.e. the 3NN V C triplet, is identified. Figure 3 illustrates the overall features of V C distributions in C-sublattice. In Zr 8 C 7 (P4 3 32) and Zr 6 C 5 (C2/m) as shown in Fig. 3(a,b), respectively, neighboring V C s only occupy 3NN sites and form the 3NN V C triplets. The schematic of 3NN V C triplet is displayed in Fig. 3(f), which is an equilateral triangle with its side length restricted by 3NN coordination. Although it was known that the preference of occupying 3NN V C shell was typical in ordered structures of sub-stoichiometric transition metal carbides and nitrides, the significance of V C triplet configuration was not reported before. We find that the 3NN V C triplet prevails in all the low energy structures like a fingerprint V C configuration.
The 3NN V C triplets are corner-shared in Zr 8 C 7 (P4 3 32) (as shown in Fig. 3(a)), and edge-shared in Zr 6 C 5 (C2/m) (as shown in Fig. 3(b)). In Zr 6 C 5 (C2/m), one defective carbon layer (filled with 2/3 carbon atoms) and its neighboring perfect carbon layer stack alternatively along < 211 > B1 crystal direction. When V C concentration is higher than 1/6, V C s could not be accommodated by occupying the sites restricted by pure 3NN coordination. Neighboring 3NN V C triplets adjust their relative positions and bring out the 1NN V C configuration to host more V C s, as shown for Zr 4 C 3 (C2/m) in Fig. 3(c). The 3NN V C triplets are predominant in the {111} B1 carbon layers and stack along < 211 > B1 crystal direction. With more V C s incorporated in Zr 3 C 2 (Fddd), 3NN V C triplets link adjacent defective carbon layers, as shown in Fig. 3(d). Additionally, self-assembling of 3NN V C triplets generates the 1NN and 4NN V C configurations. For Zr 2 C (Fd-3m) shown in Fig. 3(e), neighboring carbon layers have 1/3 and 2/3 C-sublattice sites occupied by V C s, respectively, in which crowded 3NN V C triplets are orthocenter overlapped. It is noted that every carbon atom in Zr 2 C (Fd-3m) is surrounded by overlapped 3NN V C triplets. These 3NN V C triplets are completely coordinated by 1NN and 4NN V C configurations. Changing any V C site or forming one more V C will introduce the 2NN coordination. Besides these GS structures, remarkable 3NN V C triplets are also found as the common configuration unit in near GS configurations, such as Zr 32 C 25 (C2) and Zr 32 C 23 (R-3m).
Short-range interactions. The above results illustrate the self-assembling of V C configuration and feature the fingerprint configuration of 3NN V C triplet in defective ZrC 1−x . It is speculated that the redistribution of electrons around vacancies 22,26,27 would affect the short-range interactions among V C s by altering the M-C bonds, which would be the important driven force pushing forward the local ordering pattern of V C s. Therefore, the interaction (or binding) energies of various V C clusters, and their correlations with the evolution of V C configurations, are studied here. Using the 3× 3× 3 supercell with 108 Zr sites, the interaction energies of V C pair and triplet are calculated by the following equations 28 : i C 108 108 108 106 108 107 108 108  108 105  108 107 where the first and last terms on the right side are the total energies of perfect ZrC and the supercell containing an isolated V C , respectively; and the second items in equations (1) and (2) stand for the total energies of the supercells containing various V C pairs and the 3NN V C triplet, respectively. The derived interaction energy indicates the thermal stability of various vacancy clusters relative to isolated V C s with infinite distance. and determines whether isolated vacancies would aggregate together to form certain vacancy configuration. Figure 4 plots the interaction energies of various V C pairs and the 3NN V C triplet. For V C pairs, interaction energies increase in the order of , which agrees with that reported by Razumovskiy et al. 29,30 . The interaction energy of 3NN V C triplet yields -49 meV/vacancy or − 147 meV/triplet, which is lower than that of 3NN V C pair (-22 meV/vacancy or − 44 meV/pair). Therefore, the formation of 3NN V C triplet benefits to lower the energy of defective structure. The inset in Fig. 4 sketches neighboring Zr 6 C   1NN, 2NN, 3NN, 4NN V C pairs (the circles) and 3NN V C triplet (the star); the inset shows that neighboring Zr atoms of 3NN V C triplet are endowed with the square-pyramid coordination (C 4v ).
Scientific RepoRts | 5:18098 | DOI: 10.1038/srep18098 octahedra with a 3NN V C triplet. The C atoms in octahedra are removed to produce V C s, and the 3NN V C triplet configuration ensures each Zr atom coordinated with only one V C . Therefore, the carbon coordination of Zr atom is minimally disturbed, and thereby the bonding feature is least affected. Meanwhile, the self-assembling of V C s allows short-range ordering by coordinating the most preferred V C configuration, i.e. 3NN V C triplet. Local ordering in the configuration of 3NN V C triplets would be universal in ZrC 1−x at any composition. In elastic diffuse neutron scattering experiments, ZrC 0.80 and ZrC 0.64 were identified with the same peak positions although they have different compositions 3 . This result suggests similar short-range ordering in the two compounds. If the short-range ordering covers the whole lattice and the defective structure satisfies a new symmetry at certain composition, then an ordered ZrC 1−x phase would be identified. All the results show that self-assembling of 3NN V C triplets is the key factor to maintain phase stability of defective ZrC 1−x with high concentration of V C s, and to realize short-range and long-range ordering of V C s.
The interaction energy of 1NN V C pair, 13 meV/vacancy, is slightly higher than zero but much lower than those of the 4NN and 2NN V C pairs, which are 70 meV/vacancy and 321 meV/vacancy, respectively. The small value of 1NN V C pair shows very weakly repulsive interaction. Therefore, 1NN V C pair is the second preferred configuration in the low energy structures. With more V C s accommodated in GS structures of Zr 3 C 2 and Zr 2 C, together with near GS structures Zr 32 C 23 , Zr 32 C 18 and Zr 32 C 17 , both 1NN and 4NN V C pairs present in order to efficiently coordinate neighboring 3NN V C triplets. Therefore, the presences of 1NN and 4NN V C configurations are helpful for balancing V C concentration and thermal fluctuation.
The 2NN V C pair has extremely high interaction energy (321 meV/vacancy), which means significantly repulsive interaction between V C s in 2NN configuration. Therefore, 2NN V C pair is unfavorable or forbidden in defective ZrC 1−x and it does not show up in low energy structures when V C concentration is lower than 50%. It's noteworthy that the strongly repulsive interaction among 2NN configurations may prevent the formation of large-scale V C clusters in sub-stoichiometric ZrC 1−x . Unfortunately, 2NN coordination is inevitable when V C concentration is higher than 50%. The appearance of 2NN configuration will abruptly increase the formation enthalpy (shown in Fig. 2) and therefore, the defective ZrC 1−x may undergo phase separation in such case. This might be the reason why the critical V C concentration is limited around 50% in ZrC 1−x .
The short-range interactions among vacancies could well explain the evolution feature of V C configurations and the fingerprint configuration. The result indicates that these are the driven-force of V C self-assembling, i.e. bringing down the system energy via the maximization of attractive 3NN configuration and the exclusion of strongly repulsive 2NN coordination, and simultaneously balancing 3NN V C triplets through the moderately repulsive 1NN and 4NN interactions. Also, the short-range interactions among V C s would be more and more significant at high V C concentration. It may provide hints on the origin of enhanced ordering tendency with increasing V C s.
Electronic structures of ZrC 1−x . Projected density of states (PDOS) and projected crystal orbital Hamilton population (pCOHP) 31 are illustrated in Fig. 5 to describe the Zr-C bonding characters in highly defective ZrC 1−x . At high V C concentration, the electronic structures of ZrC 1−x would be significantly affected by vacancy-vacancy interaction and/or vacancy ordering. This was verified by Pickett and Klein as they found complex differences between the calculated electronic structure of an isolated carbon vacancy in B1-NbC by the muffin-tin Green's-function method and the experimental electronic structure of NbC 0.85 from X-ray photoemission spectrum 32 . Therefore, we investigate the electronic structures of ordered Zr 2 C and Zr 16 C 15 phases to understand the electronic structures of ZrC 1−x around the critical composition of x = 0.5. This would illustrate the mechanism of electronic redistribution at high V C concentration.
Firstly, we study the bonding characters in perfect ZrC (B1) as shown in Fig. 5(a). In the perfect unit cell, Zr atom is octahedral coordinated by six C atoms with a high site symmetry (O h ), which ensures good capability of band overlapping. The bonding region in PDOS shows a strong hybridization of C p-and Zr d-states, which promotes the formation of strong covalent Zr-C bond. Besides, the orbital-wise pCOHP analysis indicates that Zr-C bond is dominated by the p-d σ and π bonding with significant share of C(p z ) and Zr ( ) − d z r 3 2 2 interactions. These bonding characters are consistent with other transition metal carbides like TiC 26,33 .
With removing of C atoms, i.e. the electron acceptors in ZrC 1−x , the excess d electrons redistribute on electronic states in high energy level, as clearly shown in Fig. 5(b). For Zr 2 C at the critical concentration, the bonding states between -2.4 eV and -4.8 eV are mainly dominated by C(p z ) and Zr ( ) − d z r 3 2 2 , and the PDOSs locating from -1.9 eV to Fermi level (E F ) correspond to the d-d bonding among Zr atoms. When the concentration of V C s is higher than 50%, for example Zr 32 C 15 (R-3m) with inevitable 2NN V C coordination, anti-bonding states originated from C(p z )-Zr(d xz ) interactions obviously show up near E F in Fig. 5(c). As a result, the formation enthalpy is quite high for ZrC 1−x with x > 0.5, which contains unfavorable 2NN V C configuration.

Discussion
We found neighboring V C s have totally different values of short-range interactions, which increase in the order of (2NN). It goes along with the evolution characteristics of V C configurations. The moderately attractive and strongly repulsive interactions between V C s in 3NN and 2NN configurations, respectively, provide the driven-force of self-assembling of V C s. Meanwhile, 3NN V C triplet is more stable than other V C configurations, and it serves as the fingerprint block in low energy ZrC 1−x structures. At high V C concentration, neighboring 3NN V C triplets modify relative positions by coordinating 1NN and 4NN V C configurations. The present results clearly demonstrate that ordering of V C s in ZrC 1−x is not an abrupt structural change. The local ordered configurations, i.e. 3NN V C triplets, already present in low energy structures at any composition; and at special composition, long-range ordering spreads throughout the whole lattice and one could observe ordered phase with new space group. The underlying mechanism falls into the concept of self-assembling of V C s. Short-range interactions among V C s are the driven-force of local or long-range ordering; and the formation enthalpy is reduced via the generation of attractive 3NN interactions and the exclusion of strongly repulsive 2NN interactions.
It is crucially important to realize the local ordering of V C s and to maintain phase stability in highly sub-stoichiometric ZrC 1−x . On one hand, local ordering of V C configuration with 3NN triplets, instead of the disordered distribution, stabilizes the defective structures by minimally affecting bonding features. On the other hand, the strongly repulsive interaction between V C s in 2NN configuration weakens the formation of large scale vacancy clusters or voids as well. More encouragingly, the self-assembling of V C s provides us the opportunity to tailor overall properties through defect engineering. One exciting perspective is the improvement of radiation resistance of ZrC 1−x via the mechanism of vacancy mediated performance optimization. Local ordering of V C s benefits the accommodation, annihilation and recombination of radiation-induced C-related point defects. Therefore, ZrC 1−x with optimal composition may be promising for next generation nuclear fuel coating and cladding material. Another important prospect is to precisely tune mechanical and thermal properties through adjusting the chemical composition of ZrC 1−x . In the near future, it demands extensive investigations on the effects of V C concentrations and configurations on radiation resistance and thermo-physical properties.
We herein proposed the physical insight of the critical V C concentration, x ~ 0.5, in ZrC 1−x . At relatively low V C concentration, excess d electrons redistribute on d-d bonding states after the removal of C atoms; whilst neighboring Zr-C bonding is not affected significantly. This mechanism is helpful to maintain the structural stability of defective ZrC 1−x . When the V C concentration is higher than the critical composition, 2NN V C configuration inevitably appears which brings out anti-bonding states between Zr and C atomic interactions. The anti-bonding states lead to extremely high formation enthalpy of defective ZrC 1−x with 2NN V C coordination and thereby, the ZrC 1−x compound loses stability and undergoes phase separation.
The current investigation suggests new clue for better understanding of the phase stability and defect engineering in non-stoichiometric interstitial compounds. Firstly, the phase stability and ordering mechanism in ZrC 1−x , especially the discovery of self-assembling of 3NN triplets driven by short-range interactions among V C s, may be common in binary transition metal carbides. There are solid evidences on the similar V C configurations in non-stoichiometric carbides, for instances, the isotopic structures of V 8 C 7 (P4 3 32) 9 and Ti 6 C 5 (C2/m) 22 with Zr 8 C 7 (P4 3 32) and Zr 6 C 5 (C2/m), respectively; the widespread occupation of 1NN and 3NN shells, but not 2NN shells in . This work presents hints to understand following puzzles, like why common characteristics of vacancy configurations exist and what are mechanisms of phase stability and V C ordering in sub-stoichiometric binary transition metal carbides. Secondly, it needs thorough inquiry whether the present results could be extended to other sub-stoichiometric binary transition metal nitrides and oxides because short-range interactions among anion vacancies are possibly different. Various self-assembling mechanisms of anion vacancies may dominate the diversity of the ordering phenomena in sub-stiochiometric nitrides and oxides.
In summary, vacancy configurations and their evolution in sub-stoichiometric ZrC 1−x were investigated by combining first-principles calculations, cluster expansion and supercell methods. Firstly, the negative mixing enthalpy and the negative ordering enthalpy are direct energetic proofs of the tolerance and the ordering tendency of V C s in sub-stoichiometric ZrC 1−x . Besides, the energetically preferred 3NN V C triplet is the fingerprint structural unit in sub-stoichiometric ZrC 1−x . To balance V C concentration and thermal fluctuation, 1NN and 4NN V C configurations with moderately repulsive interactions would show up. The tolerance of vacancies is limited at around 50% because the presence of unfavorable 2NN V C coordination with strongly repulsive interaction energy would lead to structural instability or phase separation.
It's noteworthy that the self-assembling of 3NN V C triplets driven by short-range interactions is the fundamental essence of the phase stability, short-range and long-range ordering in sub-stoichiometric ZrC 1−x . For the first time, we disclose the significant self-assembling mechanism of V C s in sub-stoichiometric ZrC 1−x , which provides us the opportunity to tailor its properties, such as the radiation resistance, mechanical and thermal properties, through defect engineering.

Methods
Sub-stoichiometric ZrC 1−x is taken as a pseudo "binary solid solution" consisting fully occupied Zr-sublattice and mixing of carbon atoms and vacancies at C-sublattice. This system can be handled by supercell 4,35 , order parameter functional method (OPF) 36 , coherent-potential approximation 37 , or cluster expansion (CE) 38 methods. Three ordered phases, Zr 2 C, Zr 3 C 2 and Zr 6 C 5 , have been predicted using OPF method 36 , but the results did not provide structural information. Theoretical calculations using supercell method 4,35 has emphasized the importance of lattice relaxation effect and the vacancy distribution, however it may leave out important configurations because of limited trial structures. The CE method which is the choice here, considers the effects from short-range interactions and local relaxation, and therefore provides both the structural characteristics and corresponding structural energy with DFT (Density Functional Theory) accuracy 39 .
Within the CE approach, the configuration on C-sublattice is described by configurational pseudospin variable σ i for lattice site i, which takes value + 1 or − 1 depending on the occupation of lattice site i by carbon atom or carbon vacancy. A particular arrangement of spins is called a configuration (or structure), which is represented by a vector {σ i } containing the value of configurational pseudospin variable for each site i. The energy of a given composition and configuration is expressed as 39,40 : where V α and ϕ α are called the effective cluster interaction (ECIs) and the correlation function of cluster α (corresponding to a set of sites i: pairs, triplets, tetrahedral, etc.), respectively. For disordered phases, the energy is only composition (x) dependent 41 : where α n stands for one kind of n-body cluster. A principal objective of the CE method is to evaluate unknown ECIs so that the configuration characters can be predicted with the accuracy of DFT calculations. Here, we determine the ECIs in the framework of the structure inversion method [42][43][44] . When energies of ordered structures are prepared from DFT calculations, the ECIs are determined by fitting a truncated form of equation (3) to the DFT energies using the least-squares technique. The accuracy of CE prediction is assessed by the leave-one-out and leave-many-out cross-validation (CV) score [42][43][44] . The set of clusters is optimized using genetic algorithm via minimizing the CV score. The DFT structures include initial selected random structures, low energy structures refined during construction of CE, and reported superstructures in non-stoichiometric interstitial compounds 9,45 . Finally, 25 clusters and 250 structures were selected to construct the CE formula. The average prediction error of the final CE is less than 4.3 meV/cation (per Zr/C-sublattice site). Thereafter, two types of calculations are used to search low energy structures (LESs) via the constructed CE. One is an exhaustive search by calculating the energies of all possible configurations in a finite-sized cell. The other is the simulation annealing (SA) method 46 within a larger structural space. Here, exhaustive search was used for the supercell containing 32 or less Zr sites and the SA method was used for the supercell containing up to 1728 Zr sites. Detailed procedure of constructing CE can be found in ref. 47. In this work, CLUPAN code 42,48 is used for the construction of CE and the searches of LESs.
First-principles calculations were done using VASP 49 code, in which the projector augmented-wave (PAW) method 50,51 within generalized gradient approximation 52 was employed. The plane-wave basis cutoff and the k-mesh separation were set as 500 eV and 0.04 Å −1 , respectively. Full structural relaxations (atomic positions and lattice constants) were performed until the energy difference converges to less than 10 −6 eV.