Introduction
Diamond has long been studied due both to its hardness and natural beauty. In regards to the latter, it is well known that certain imperfections give rise to unique colors, and it is this capacity to host a number of color centers that has brought diamond under increased scrutiny of late. In particular, the fields of quantum sensing and quantum information processing (QIP) have recently prompted considerable investigation, and diamond is a viable candidate for both. For QIP, a bright, photostable, room temperature single-photon source (SPS) is desirable.[Reference Aharonovich, Englund and Toth1] Nitrogen is one of the most common dopants in diamond, being present during various synthesis pathways, and a nitrogen-vacancy (NV) center that is capable of single-photon emission has been thoroughly studied.[Reference Doherty, Manson, Delaney, Jelezko, Wrachtrup and Hollenberg2] NV centers come in two varieties, NV0, with an optical zero phonon line (ZPL) at 2.156 eV,[Reference Davies3] and NV−, which shows a ZPL at 1.945 eV.[Reference Davies and Hamer4] While the discovery of these NV centers marked a significant advancement for the development of quantum technologies, their use as SPSs is compromised by the existence of a large phonon sideband, which creates challenges for technologic implementation.[Reference Orwa, Greentree, Aharonovich, Alves, van Donkelaar, Stacey and Prawer5]
Another promising series of diamond color centers are formed by the introduction of nickel. Often accompanied by nitrogen doping, this transition metal can form a variety of defect sites, characterized by the degree of nitrogen co-doping, vacancy count, and the relative position of the atoms.[Reference Nadolinny, Yelisseyev, Baker, Newton, Twitchen, Lawson, Yuryeva and Feigelson6] Electron paramagnetic resonance measurements show that NE8 is comprised of a nickel atom located interstitially, with four nearest neighbor nitrogen atoms. This NE8 color center serves as a bright, photostable single-photon emitter at 1.56 eV,[Reference Gaebel, Popa, Gruber, Domhan, Jelezko and Wrachtrup7, Reference Aharonovich, Castello, Simpson, Su, Greentree and Prawer8] and it has generated interest for quantum photonics due to its narrowband emission. Prior synthesis approaches include chemical vapor deposition (CVD)[Reference Rabeau, Chin, Prawer, Jelezko and Gaebel9, Reference Wu, Rabeau, Roger, Treussart, Zeng, Grangier, Prawer and Roch10] and ion implantation,[Reference Aharonovich, Zhou, Orwa, Castelletto, Simpson, Greentree, Treussart, Roch and Prawer11] which are often combined with high-temperature annealing processes,[Reference Orwa, Greentree, Aharonovich, Alves, van Donkelaar, Stacey and Prawer5] but precise control of the resultant defect structure remains elusive. The need for a scalable bottom-up growth process points to microwave plasma-enhanced CVD (MPCVD) as a method for the production of color centers embedded in nanodiamond particles.[Reference Zyubin, Mebel, Change and Lin12]
Much of the theoretical work related to diamond color centers has understandably been directed toward developing a detailed understanding of specific color centers in diamond. However, the inability to precisely control the defect structure during synthesis combined with the lack of knowledge regarding reaction kinetics in the MPCVD environment prompts the need for thorough computational study of potential defect compositions, structures, and optical properties. To initiate such a study, we investigate nickel–nitrogen defect sites in nanodiamond, focusing on energies of formation and optical absorption. This work elucidates thermodynamically stable nickel–nitrogen defects in nanodiamonds, irrespective of the synthesis method. These results should prove useful for focusing theoretical work on as-yet unstudied target compositions that are thermodynamically realistic.
Method
Diamond defects may be modeled by supercell and cluster[Reference Zyubin, Mebel, Change and Lin12] approaches. The current work takes a cluster approach for the dual purposes of avoiding concerns regarding the defect–defect interactions common to supercell approaches while also developing a model of the end product of the MPCVD synthesis effort—namely, an isolated color center embedded in nanodiamond. Since computational resources scale non-linearly with cluster size, a size study has been performed on pristine diamond clusters, and C87H76 with a diameter of ~0.9 nm has been found to represent a good balance between computational requirements and quantitative results.[Reference Gothard, Dudis and Bissell13] This also agrees with previous studies that indicate quantum confinement effects begin to vanish about 1 nm sizes for pristine nanodiamond.[Reference Raty and Galli14]
The base C87H76 cluster we have developed has T d symmetry, which matches the site symmetry in bulk diamond. Into the pristine C87H76 carbon clusters, we introduce nickel via two optimization pathways: (1) substitutionally and (2) at an interstitial semi-vacancy, i.e., between two lattice sites. Then, we systematically vary the composition of nickel-nearest neighbors (either nitrogen atoms or vacancies) to generate a thorough sampling of compositions and defect structures. In all cases, the cluster surfaces were passivated by hydrogen atoms.
Substitutionally doped clusters take stoichiometries of the general form NiNxVyC86−x-y, where x indicates the number of nitrogen-nearest neighbors and V indicates nearest neighbor vacancies. Substitutional placement of a nickel atom creates tetrahedral bonding to four nearest neighbors. Due to the high symmetry of the underlying C87H76 cluster, there are then 15 possible unique clusters. Placing a nickel atom at a semi-vacancy site, however, breaks the symmetry and changes the nearest neighbor count for nickel from four to six. The stoichiometry of these clusters is also NiNxVyC86−x-y, where nickel is now at a semi-vacancy site. If we neglect defect orientation, this allows for 28 possible interstitial defect site configurations. With such a large structural space, we only nominally consider variations of charge and spin. Except where specified, we treat the nickel-vacancy stoichiometries as charge neutral. Further, since we do not specifically address spin effects, electrons are coupled such that there is at most one unpaired electron. Figure 1 depicts representative substitutional and interstitial nanodiamond clusters. Hydrogen terminations are not depicted.
(a) Substitutional NiN4C82H76 cluster, (b) interstitial NiN4C81H76 cluster. Carbon atoms are pictured in dark gray, nitrogen atoms blue, and nickel atoms light gray. Terminating hydrogen atoms are not depicted for clarity.

Following cluster development, the geometries were structurally relaxed in the ground state using restricted Hartree–Fock wave functions and Pople's 6–31G* split valence basis set, using the GAMESS software package.[Reference Schmidt, Baldridge, Boatz, Elbert, Gordon, Jensen, Koseki, Matsunaga, Nguyen, Su, Windus, Dupuis and Montgomery15] Since introducing symmetry during structural relaxation represents an external constraint on the system, substitutional clusters were structurally relaxed both with and without symmetry constraints to allow nickel to migrate away from a lattice site. In the case of the interstitial defects, the defect sites break the long-range symmetry in nearly all cases, and no point group constraints were therefore assigned during structural relaxation.
The defect site energy of formation was estimated by comparing the energy of the doped cluster to the energy of the pristine C87H76 cluster, using
where E PROD is the binding energy of the doped cluster, N C is the number of carbon atoms ejected from C87H76 in order to form the doped cluster, N N and N Ni are the number of nitrogen and nickel atoms introduced into the pristine cluster, and E C, E N, and E Ni are the carbon, nitrogen, and nickel self-energies of the isolated atoms.
Time-dependent density functional theory (TDDFT) was chosen as a suitable rapid screening method for characterizing the optical properties of these defect sites. DFT-based approaches are, however, sensitive to the choice of functional, and previous benchmarking of TDDFT results has shown that while pure generalized gradient approximation functionals such as PBE[Reference Perdew, Burke and Ernzerhof16] approximate well the bulk diamond energies for C87H76, a hybrid functional is needed for diamond clusters that incorporate transition metals. Benchmarking has shown good agreement between three such functionals: B3LYP,[Reference Stephens, Devlin, Chabalowski and Frisch17] PBE0,[Reference Adamo and Barone18] and X3LYP,[Reference Xu, Zhang, Muller and Goddard19] and we chose the PBE0 functional for TDDFT excitation.
Discussion
Energy of formation is found to increase with vacancy count, as shown in Fig. 2, while nitrogen counteracts this effect when present at moderate concentrations. Formation energies are lowest for clusters with zero vacancies that have 1–3 nitrogen atoms, and the average formation energy of 13.86 ± 0.18 eV reflects nearly equal likelihood of formation for any structures with these compositions. A similar effect manifests for clusters with single vacancies that have two to four nitrogen atoms. When more than one vacancy is introduced, the formation energy continues to increase past 30 eV, representing increased thermodynamic unlikelihood. We also note that variations in formation energy between clusters with the same composition are due to the imposition of symmetry during structural relaxation. NiC86H76 is a remarkable example of this effect. The dramatic increase in formation energy for this cluster when symmetry is enforced indicates that the actual structure prefers a substantial reduction in symmetry. On the other hand, introducing vacancies or nitrogen atoms significantly decreases the overall cluster symmetry, resulting in a smaller formation energy difference between the optimization pathways. We note also that formation energies do not appear to vary substantially with cluster size, as shown in Table I.
Formation energies as a function of vacancy and nitrogen count, for both substitutional and interstitial clusters. Substitutional clusters were relaxed with and without symmetry constraints. As the formation energy increases with vacancy count, clusters with more than two vacancies are not presented.

Variation of formation and TDDFT energies as a function of cluster size, for nickel at interstitial sites with and without nitrogen co-doping.

Table II presents data summarizing the interplay between defect structures and formation energies. Due to the differences in stoichiometry and defect design, the defect sites exhibit varying degrees of structural rearrangement when allowed to relax. A key factor in gauging the degree of local structural rearrangement is the ability for nickel to diffuse through the cluster upon relaxation. When high-symmetry point groups are employed, or when there are no vacancies in the structure, fewer degrees of freedom exist for nickel diffusion during structural relaxation, while the introduction of vacancies creates additional diffusion pathways for the nickel atom. In the latter case, the nickel atom will drift toward an interstitial site, as constrained by the choice of point group. But when symmetry is removed, nickel tends to migrate into a fully interstitial position, regardless of the starting structure. As the Ni–C bond length is shorter than the Ni–N bond length, the precise structural configuration of the relaxed defect site depends on the number of vacancies as well as the number and location of nitrogen atoms. While the tetrahedral carbon–carbon bond distance in pristine diamond is 1.54 Å, placing a nickel atom at a semi-vacancy interstitial site gives the nickel atom six nearest neighbors at a distance of 1.94 Å. When the nickel-nearest neighbors are all carbon atoms, the nickel–carbon-nearest neighbor separation relaxes to 2.066 ± 0.009 Å. On the other hand, when all six nearest neighbors are nitrogen, the average nearest neighbor separation becomes 2.137 ± 0.017 Å. Thus, when undergoing structural relaxation, nickel will attempt to diffuse preferentially toward carbon rather than nitrogen.
Effect of structural relaxation on nickel atom position for symmetry-enabled substitutional (S), symmetry-disabled substitutional (SD), and interstitial (I) optimization pathways.

Starting nickel position for substitutional clusters is (0, 0, 0) and for interstitial clusters (0.45, 0.45, 0.45). All lengths are given in angstroms.
The preferential migration of nickel toward carbon is also seen for compositions with partial nitrogen filling of the nickel-nearest neighbor sites, such as the NiN2VC830 or NiN4VC810. In these clusters, the average nickel–carbon bond lengths approach 1.95 ± 0.04 Å, while the average nickel–nitrogen bond length approaches 2.11 ± 0.02 Å. The actual equilibrium position of the nickel favors proximity to the carbon-nearest neighbor. When the nearest neighbor composition forces comparatively larger structural rearrangements, as occurs for the NiNVC840 or NiN6VC790 sites, there is a corresponding increase in formation energy, which may reflect longer range structural distortions. For defects following the semi-vacancy interstitial configuration, there is little adjustment to the starting structure during the relaxation process. Combining this fact with the comparatively low energy of formation for this series of clusters indicates these to be likely for spontaneous formation. Nevertheless, there are a number of stoichiometries and structural arrangements in the range of NiN2VC83H76–NiN4VC81H76 with nearly identical structures and formation energies.
The TDDFT absorption energy results for the low formation energy structures are presented in Fig. 3. Variance of TDDFT energies with cluster size has been previously considered,[Reference Gothard, Dudis and Bissell13] and convergence to near the bulk value observed for pristine C87H76 clusters. While both the zero- and single-vacancy clusters show low formation energies and similar structural evolutions for a number of stoichiometries, the absorption energies vary over a range of nearly 2 eV for both zero- and single-vacancy structures. The response of the absorption energy to nitrogen concentration is complicated by the spin multiplicity: introducing an odd number of nitrogen atoms creates a spin doublet, while even nitrogen counts represent singlet or triplet states. Many interdependent parameters influence absorption energies: defect charge, spin state, dangling bonds, and structural distortions. A thorough analysis of these dependencies is beyond the current scope.
TDDFT excitations for nickel–nitrogen complexes with fewer than three vacancies. Energies for zero-vacancy structures correspond to nickel at a lattice site, while others correspond to nickel at a semi-vacancy site.

Table III compares absorption energies estimated here to accepted values. The AB5 (NiNC85−) site is in good agreement with experimental absorption data,[Reference Neves, Pereira, Sobolev, Nazare, Gehlhoff, Naser and Kanda20] which show strong absorption near 1 eV. With regard to the NE8 center, we note that literature values of 1.56 eV correspond to the ZPL in bulk diamond.[Reference Gaebel, Popa, Gruber, Domhan, Jelezko and Wrachtrup7] Previous size effect studies of the NE8 center in nanodiamond have shown absorption to approach 2.2 eV as the size is reduced as far as 1.1 nm.[Reference Thiering, Londero and Gali21] Our data showing optical absorption at 2.45 eV for a 0.9 nm cluster are therefore in reasonable agreement. In addition, while the NE4* site should undergo absorption above 2.5 eV, our study shows the onset of absorption near 2.14 eV. Preliminary studies of other transition metal complexes show similar differences in absorption and emission energies.
Comparison of TDDFT absorption energies to literature results.

There is considerable work still to be done to understand the interplay between defect site compositions, structures, and photoluminescence in nanodiamonds. Nitrogen is a common impurity during synthesis and is shown to lower the formation energy when incorporated in moderate amounts. This approach represents a rapid screening of potential nickel–nitrogen sites in diamond. While it does not predict the emissive properties, it nevertheless represents a rapid-screening approach to the more thorough investigation of other potential color centers in diamond, since absorption energies represent an upper bound to luminescence energies.
Conclusions
A formation energy/TDDFT approach was used to systematically screen nickel–NV complexes in hydrogen-terminated nanodiamond. We find that the introduction of a nickel atom favors the creation of a single-vacancy, although more may be created depending on the reaction kinetics, and that nitrogen is a favorable co-dopant with nickel. When introduced at a lattice site, nickel prefers to migrate into an interstitial position, and a nickel–nitrogen semi-vacancy structure appears both energetically favorable and structurally preferred. TDDFT of ground state structures provides an upper bound to luminescence energies via an estimate of the absorption energies, which are useful from a standpoint of broadly guiding more rigorous cluster development.
More detailed analysis of thermodynamically likely stoichiometries is ongoing and includes factors such as charge state and spin polarization, as these can have significant impacts on the optical properties. The present results provide a useful baseline for understanding low-spin systems. Additionally, multi-configuration approaches are needed to deduce the fine electronic structure and luminescence properties, which are of more interest for quantum information applications. In addition, there is considerable modeling to be done in order to elucidate the reaction kinetics in a plasma environment, as well as to determine the stability of these defect sites during subsequent annealing. Finally, we note that size effects are expected in the absorption and luminescence of color centers in nanodiamond, and these results should be compared with care to color centers in bulk nanodiamond. Future work is directed toward modifying and applying this approach to initialize a detailed study of other color centers in diamond.
Acknowledgments
The authors cordially acknowledge funding by AFOSR, “Single Color Center Engineering,” Program Manager: H. Weinstock; and the DoD Quantum Science and Engineering Program. N. Gothard additionally acknowledges support from the AFOSR Summer Faculty Fellowship Program. Provision of computational resources by the DoD High Performance Modernization Office is gratefully acknowledged.





