Electronic features of vacancy, nitrogen, and phosphorus defects in nanodiamonds

Defective nanostructures with a surface termination are the focus of this work. In order to elucidate the influence of the defect on the properties of nanomaterials, we take hydrogen terminated nanodiamonds. Various vacancy defect centers are separately embedded in a nanodiamond at different positions. These include some of the known defects, such as the charged nitrogen-vacancy (NV−), the silicon-vacancy (SiV0), the germanium-vacancy (GeV0), the phosphorous-nitrogen (PN), and the nickel-vacancy (NiV−). For these defective nanodiamonds, we probe the influence of the defect type, its position, as well as the size of the nanodiamond through their structural and electronic features. A detailed and comparative analysis is provided here, based on quantum mechanical simulations. Our results shed light into the inherent differences of these defects in nanodiamonds, allowing for a better understanding of defective nanostructures. In the end, we discuss the potential of tuning their characteristics in view of novel nanotechnological applications.

carbon atom and a nitrogen atom replacing a nearest neighbor or a next-nearest neighbor carbon atom, respectively [32,35]. Both these defects have a ground state with a total spin of S = 1/2.
Overall, the possible applications of defective nanodiamonds beyond the NV center range from biomedicine over quantum computation and information to single photon emitters [36][37][38]. For these, a detailed understanding of the defect centers and their inherent properties, such as their stability, charge and spin states is still needed. Of high importance is the investigation of surface and size effects of the nD [39][40][41][42], the effect of the defect position within the nD [25,43], the type of defect [26,29,44,45], the influence of functionalization and modification of the nD [6,46,47], etc. Along these lines, this investigation focuses on the analysis of the structural and electronic properties of nanodiamonds including different vacancy-related and phosphorus defects. The distinct characteristics of each defective nD are compared in order to identify trends. Our analysis will be presented with respect to the defect type, its position within the nD, and the size of the nD. Accordingly, this paper is organized as follows: the methodology used together with the defects studied in this work are outlined in section 2, the results are discussed in section 3, while in the end we present a summary of our work and an outlook in section 4.

Methodology
Computer simulations based on density functional theory (DFT) [48,49] as implemented in the code SieSta [50] were carried out. For the exchange-correlation functional, the local (spin) density approximation (LSDA) of Ceperley and Alder [51] as parameterized by Perdew and Zunger [52], as well as the generalized gradient approximation (GGA) of Perdew-Burke-Erzernhof (PBE) [53] and non-relativistic norm-conserving Troullier-Martins pseudopotentials [54] were used. In comparing LDA with GGA we found significant differences especially in the energy levels and energy gaps that will be discussed later on, so that in this article, only the GGA results will be presented. For expanding the Kohn-Sham states, we have considered a polarized double-ζ basis-set (DZP) for all atoms and a real space sampling (mesh cutoff) grid above or equal to 250 Ry. For such systems, we have tested the use of a more accurate triple-ζ with polarization basis-set (TZP), which was found only to increase the computation time without changing quantitatively the results. At this level of our proof of principles study, the use of atomic-like orbitals instead of plane waves is sufficient. The results for both cases are qualitatively very similar around the bandgap region and the frontier orbitals are very much localized, justifying the use of atomic-like orbitals. Further studies on a more accurate resolution of the electronic structure of the nDs will utilize a plane wave basis set. Based on this basis sets comparison and the qualitative picture on the trends in defective nDs we aim to provide, we choose the computationally less demanding DZP basis set. A total charge was set in each calculation. The extra charge is added in the whole computational cell and to both spin channels. Note, that partially occupied orbitals are allowed. The addition of a charge q add is done by adjusting the Fermi level so that EF(q=0) EF(q=q add ) DOS( ) d = q add , where DOS is the electronic density of states, E F (q = 0) is the Fermi level of the neutral system, DOS( ) the density of states at energy , and E F (q = q add ) the adjusted Fermi level due to the presence of the additional charge q add . The additional charge affects in particular the separate spin states and thus the obtained eigenenergies. As an example, the ground state of the NV − defect has a total spin S = 1, while the NV 0 has a total spin of S = 1/2, leading to different eigenenergies [24,55]. The structural relaxations were performed until the net forces of each atomic component become smaller than 0.04 eV Å −1 . The structures were built using the software VirtualNaNolab [56] and geometrically optimized using the conjugate gradient (CG) algorithm. All calculations are spin-polarized, with the total spin initially not set. In these calculations, the Kohn-Sham states are split into two sets, one for the spin-up and one for the spindown configurations equally representing these spin configurations. A separate exchange-correlation operator for each spin results in different Kohn-Sham Hamiltonians for each spin and the Kohn-Sham orbitals for each spin are calculated self-consistently. For speeding up the calculations, we have first optimized the structures using the spin-polarized version of LDA. Once this relaxation is performed, the nDs are further relaxed under the PBE approximation.
We have prepared our structures starting from a perfect bulk diamond lattice. For this, we have relaxed its geometry in order to obtain the optimum lattice constant for both exchange-correlation functionals. In order to assess the nD size effect, five quasi-spherical nDs of various radii were cleaved from the bulk crystal created with the corresponding lattice constants. In order to do so, we have first set the desired diameter. Then by starting from a certain carbon atom, we have removed all other carbon atoms at a distance larger than this diameter. In this way, the resulting structures are close to spherical. Note, that thermodynamically stable hydrogenated nDs with a very small diameter are almost spherical [57]. Detonation nanodiamonds of only a few nm in diameter are also close to spherical [58,59], while quasi-spherical NV-rich nanodiamonds less than 10 nm in diameter have also been fabricated for fluorescent applications [60]. The diameters of the nDs range d ≈ 1.2-2.4 nm, corresponding to a total number of atoms between 222 and 1474. Four of the nDs with a defect placed at the center are depicted in figure 1. The nD surface dangling bonds were passivated with hydrogen atoms. Note, that hydrogen-terminated diamonds [61][62][63][64] show a quenching of the nitrogen-vacancy defect emission and care should be taken in choosing a good functionalization [65]. For the current study though, the photostability of the defect is not of relevance. Our aim is to provide a qualitative comparison of the different defect centers and not discuss at this level dynamical or excitation properties. The goal here is only to unravel electronic characteristics by taking the hydrogenated nanodiamond as a generic host for the defect.

Defect types
Briefly, the defects included in our study are the NV − (negatively charged nitrogen-vacancy), the SiV 0 (neutral silicon-vacancy), the GeV 0 (neutral germanium-vacancy), and the NiV − (negatively charged nickel-vacancy) known as the nickel split-vacancy defect or double semi-vacancy [66]. Representative of non-vacancy related defects, also the PN + (positively charged phosphorous nitrogen) both in the nearest neighbor, and next-nearest neighbor configuration of the P and N defect will be considered in this work. The structures of these defect centers after their geometry relaxation are shown in figure 2. Note, that only one defect is embedded in each nanodiamond. The total spin S for the hydrogenated nDs with various defects after the relaxation was found to be S(NV − ) = 1, S(SiV 0 ) = 1, S(GeV 0 ) = 1, S(NiV − ) = 1/2, S(PN + nn ) = 1/2, S(PN + nnn ) = 1/2, which are indeed the ground states as reported by other studies [26,29,32,34,55,67,68]. Note, that there are a vast number of known intrinsic or extrinsic defects in diamond [69]. It was not possible to study here all these. Instead, we have taken a certain number of representative extrinsic effects for our comparative study.
Specifically, the NV − defect consists of a nitrogen atom replacing a carbon atom in the diamond lattice and a neighboring vacant site. The neutral SiV 0 and GeV 0 defects have a geometry similar to that of the NV − center, where the impurity atom is Si or Ge, respectively and not N. Note, that the N atom in the NV − defect relaxes away from the vacancy, while the Si and Ge atoms of SiV and GeV relax towards it, as observed also in figure 2. The NiV − defect is made of a nickel atom that replaces a carbon site in the diamond lattice and has two neighboring empty carbon sites (two vacancies). The nearest neighbor PN + and next-nearest neighbor PN + are defects that replace a pair of carbon atoms in the lattice by a nitrogen and a phosphorous atom at a nearest neighbor or at next nearest neighbor sites, respectively. In order to assess the influence of the defect position and its distance from the nD surface, the defects were embedded in the nD at four distinct radial positions as depicted in figure 3 for the NV − defect. Other defects were shifted along a similar path. Note, that the orientation of each defect with respect to the normal vector of the surface is kept fixed for all defect positions.

Results and discussion
We proceed with a comparison of the different defect centers, being interested in the structural and the electronic features of the hosting nDs. Focus is given on the trends in the electronic energy levels, using also the charge density for a qualitative understanding of the distribution of electrons around the defect centers. The behavior of the frontier orbitals, the highest occupied (HOMO) and lowest unoccupied molecular orbitals (LUMO) will be investigated in more detail, as their difference leads to the technologically relevant energy band-gap. We analyze the results first with respect to the nD size, the defect type, and its position within the nD.

nD size effect
We begin with an NV − defect placed at the center of the nD. The respective electronic energy levels are depicted in figure 4. In order to promote comparison among the different sizes and later the different defect types, we use  (1474), respectively. As a representative defect, the negatively charged nitrogen-vacancy defect at the center is taken. The defect symmetry axis is normal to these views. The carbon, nitrogen, and hydrogen atoms are shown in gray, blue and white, respectively. This color coding will be used throughout.
Electron. Struct. 1 (2019) 025002 the pure (non-defective) nD as a reference. For this, we have shifted in the respective graphs all energy levels, so that the LUMO (E LUMO pure ) of the pure nD lies at 0 eV. This is not the only choice as relevant studies in the literature show. Other energy levels, such as the vacuum energy are also used. Our aim, though, is to reveal changes in the properties of the pure nanodiamond when this is doped. Since the LUMO level can be detected experimentally and is of high relevance to applications involving optical properties and excitations, we chose this as the reference. As a first observation, the electronic structure of the NV − defect center close to E LUMO pure can be identified. This includes the defect levels a 1 (1), a 1 (2), e x , and e y as in the energy level scheme for a bulk defective diamond [44]. Six electrons are occupying the HOMO to HOMO-3 levels for the spin-up and the HOMO and HOMO-1 level for the spin-down configurations. In the spin-up configuration, all the defect levels lie below E LUMO pure . In the spindown configuration, the Fermi energy of the NV − (not shown) lies between the a 1 (2) and (e x , e y )-levels. In the figure, the defect levels are identified through the degenerate e x and e y levels close to the Fermi level (not shown) for the respective defect, and the occupations of the energy levels for spin up and spin down states. Relative shifts are observed in the energy levels for the different nD sizes, but the defect electronic structure is qualitatively kept the same and is only slightly perturbed. Interestingly, the energy levels of the larger nD lie below those of the smaller nDs. This behavior is reversed above the E LUMO pure level. In that regime, there is a clear trend of increasing energies for larger nDs. The spacing between the energy levels, which are not the defect levels increases as the system size decreases. This is known as the quantum confinement effect related to the 'more discrete' energy levels, as the electrons are confined within a nanometer-size nD. A comparison of the spin up and spin down , and (f) PN +1 (next nearest neighbor). The carbon, nitrogen, silicon, germanium, nickel, and phosphorous atoms are shown in gray, blue, ochre, yellow, green and purple, respectively. The vacancies are indicated by light pink spheres (at the coordinates of the missing carbon). This color coding will be used throughout. Only a section of the nD around the defects is shown for visibility. Based on these features, we can extract the electronic band-gap for all nD sizes and both spin configurations. Note, that our methodology does not lead to quantitatively exact values for the band-gaps, but can reproduce trends. In figure 5, the energy gaps for both spin configurations are shown as a function of the nD diameter. While the spin down energy gap seems not to change considerably, the spin up energy gap strongly decreases from around 2.3 eV to 1.8 eV with increasing diameter d. This can be understood by looking again into figure 4. In this, the three apparently size-invariant energy levels corresponding to a 1 (2), e x and e y lie below the NV − Fermi energy (not shown) in the spin up, but in the spin down case, the same Fermi energy lies between the a 1 (2) and (e x , e y )-levels. This leads to a change in the respective energy difference in the latter case. The charge density analysis (not shown) reveals that the electrons localize mainly around the N atom due to its higher electronegativity compared to C. The main difference between the spin up and spin down charge density is found at the carbon atom that neighbors the vacancy. There the charge density is slightly higher for the spin up than for the spin down case. No significant variation in the charge density was observed with respect to the nD size.

Influence of the defect type
Moving on to the different defect centers considered here, we summarize in table 1 the average bond-lengths and bond-angles around the defect. These are measured on the impurity atoms with respect to their three closest carbon atoms of the nD host, A is the defect element and C 1 , C 2 , C 3 the three carbon atoms). The average bondlength is the mean value of the three distances r b (AC 1 ), r b (AC 2 ), r b (AC 3 ), and the three bond angles ϑ(AC 1 C 2 ), ϑ(AC 2 C 3 ), ϑ(AC 3 C 1 ), respectively. In case of the PN-defects, A=P and A=N are shown separately in table 1. Note, that for the PN + nn (PN + nnn ) defect, the distance between P and N at the center was found equal to 1.950 Å (2.764 Å ) and close to the surface 1.953 Å (2.761 Å ). As a first observation from this table, the NV − defect is more bound to its host nD, as its distance to the neighboring carbon atoms is the smallest of all defect cases. On the other side, the GeV 0 defect is the one which seems to repel the host C atoms. This holds when the defects are embedded into the center of the nD. Towards the other end, a defect close to the surface of the nD, the bonds of the defects to the neighboring C atoms slightly increase, while on average the respective bond-angles are either similar or roughly decrease. These trends indicate, that the bonding environment of a defect can be further varied if large nDs are considered. This would also be related to a change in their stability, thermal or optical, with size.
The electronic energy levels of the various defects at the center of the nD are given in figure 6. All defective nDs are compared to the pure hydrogenated nD of the same size without any defect. The numerical values of the defect level energies show large variations across different defect types (and the pure nD). The presence of the defect centers are first recognized through the fact that the (almost) six-fold degeneracy of the pure nD below the E LUMO pure level is broken. This occurs for all cases in their spin up configurations. In the spin down configuration, all defects apart from the NV − and the NiV − show the six-fold symmetry just below their Fermi level. In both spin states the NV − reveals levels closer to the E LUMO pure level, followed by the PN + nn , the NiV − , and the PN + nnn at the HOMO level. For the spin down state, the NV − levels are followed by the NiV − ones, while all other defect levels come as two degenerate groups more than 1 eV lower. This richness and complexity of the electronic energy levels of these defects underlines the choice of a different defect for different types of applications. The respective HOMO-LUMO energy gaps are summarized in figure 7. In some cases the band-gap for the spin up was lower than that of the spin down configuration, a trend which was reversed for other defects. Apart from the NV − center, in all other defect cases, much larger differences between the spin up and spin down states are observed. The smallest band-gaps are observed for the NiV − , while the larger ones are found for the SiV 0 and the GeV 0 centers in their spin up states. In the spin down state, the larger gaps were found for the PN + (nnn) . We next examine the distribution of the electron charge density around the defect. Contours of the electron number densities of all defects embedded at the center of the nD are depicted in figure 8. The NV − defect shows an increased electron density around the carbon atom closer to the vacancy, indicating a stronger local magnetization on that carbon site. The nitrogen atom has relaxed away from the vacancy, increasing its distance from 1.54 Å to 1.70 Å with a relative change of 0.16 Å (10.8%). As the Si and Ge atoms have the same valence configuration, the charge densities of the SiV 0 and GeV 0 defects display similar distributions. Both the Si and Ge atoms have relaxed towards and very close to the vacancy reducing their distance by 49.9% and 50.2%, respectively. The neighboring to the defect C atoms show an increased charge density at the side facing the Si/Ge atom, indicating the position with the highest local magnetization around the defect. This amount of charge density is higher than around other carbon atoms. Carbon is more electronegative compared to Si/Ge and accumulates more charge moving it away from the substituted elements. In the case of the NiV − defect, a clear charge accumulation further away from the nickel site is observed, with the a larger component found at one of the closest carbon atoms. In the nearest neighbor PN + defect the amount of accumulated charge around the nitrogen atom is the highest, while a very small contribution comes from the phosphorus atom. In the next nearest neighbor PN + defect, the largest local magnetization is found at the N and C sites close to P with a slight shift towards carbon. These findings

Defect type
Center Surface are physically intuitive considering the electronegativity order starting from nitrogen going to carbon, and to phosphorus. This can explain the accumulation of the charge around the intermediate carbon atom, as electrons move from the P atom to the intermediate C atom, and part of those from this C atom to the N atom, locally magnetizing those sites.
In order to provide a deeper understanding on the electron behavior in the defective nDs, we analyze the contrib utions of the defects on the electronic density of states of the nDs. The results are summarized in figure 9 and clearly reveal the states in the pure nD gap that the defect centers introduce reducing the band gap of the defective nDs. Due to its metallic nature the nickel impurity not only provides more states in the gap, but these also have much more pronounced peaks compared to the other impurities. A comparison of the spin up and spin down configurations reveals small shifts towards higher energies for the latter. As in all other results, the SiV 0 and GeV 0 show similar trends with a peak close to the HOMO arising from the carbon atom neighboring the defect. The PN + nn and PN + nnn defects provide states in the gap arising from the nD matrix and their C neighbors, while the P atoms contribute more than the N atoms. Next to the NiV − center, the NV − also shows rich intra-gap features. These arise from the carbon atoms neighboring the defect with a small contribution also from the N impurity. As opposed to the NiV − defect with more states close to the HOMO, the NV − provides states closer to the LUMO,  especially in the spin down configuration. The inset of the figure more clearly represents these states for the NV − , showing that the higher contribution close to the LUMO comes from the carbon neighbors of the defect, while the impurity N leads to less states closer to the HOMO. The shift of the spin down states closer to the LUMO is again evident. Accordingly, provided that enough energy is given, an electron can more easily hop via these intergap states to the LUMO affecting the optical and excitation properties of the nD.

Defect position in the nD
The electronic structure related to different NV − positions in the nD (see figure 3) is depicted in figure 10 and are compared to those of a pure (without a defect) hydrogenated nD of the same size. For all defect positions the defect levels can be clearly identified. For both spin configurations and all defect positions, the HOMO levels lie approximately 2 eV higher than the pure nD eigenvalues. The LUMO level of the spin up configurations is very close to the pure nD energies. The two degenerate LUMO levels for the spin down configuration move even closer  to the corresponding pure nD levels. The spin up LUMO state is closer to the pure nD LUMO than the spin down LUMO state indicating a localization of the spin up LUMO wavefunction at the nD surface. A comparison of the levels for all positions reveals that in all cases the defect levels are energetically very close and follow the pure nD trends apart from the band gap region. The six-fold symmetry of the bulk diamond HOMO levels seems to slightly break in the case of a pure nD. This can be justified by its very small size and the hydrogen termination. The defective nD levels for spin down below the HOMO level lie about 0.3 eV higher than those of the pure nD. This difference increases towards the band gap, but the defect levels remain always very close to each other despite the different defect positions. In the spin up configuration, above the LUMO the levels for a defect position at the nD center are on average about 76 meV higher than for a position closer to the surface. For the spin down configurations, the defect energies lie lower than the pure nD levels and the energy for the defect closer to the nD surface lies higher than all the others. This behavior is reversed at higher energies.
Finally, for all defects and their positions in the nD with a diameter of 2.4 nm, their electronic band-gaps are summarized in figure 11. The variety of different band gap values covered by the different defect centers is obvious. These span the range from about 0.5 eV to just below 4 eV. In the spin up configurations a decrease in the band gap is observed as the defect is placed closer to the surface for the NV − , PN + nn , and PN + nnn . An increase is found for the SiV 0 , GeV 0 , and NiV − . In the spin down configuration, all defects show a decrease in the band gap for shallower defects apart from the NiV − . This exception could be assigned to the metallic nature of this defect. Overall, the change in the band gap as the defect moves from the center to the surface of the nD is larger in the spin down configuration reaching about to about 10% for the PN + nn case. These variations are related to the differences in the occupations of the defect levels in the gap for the two spin components.

Conclusions
Using quantum-mechanical simulations, we have investigated in detail the influence of various defect centers on the structural and electronic properties of tiny quasi-spherical hydrogenated nanodiamonds. Let us note, that our aim is not to resolve the exact electronic structure of the nDs nor their dynamic evolution or thermal stability. At this level, the photostability of the defect is of no relevance. We use the defect as a means to introduce defect levels that shift the frontier orbitals of the nD and could be differently affected by the defect type. Accordingly, the influence of the type of the defect, its position within the nD, and the nD size were investigated in detail. Our study has provided a detailed comparison of in vacuo various defective nDs and could identify differences in their electronic structure and the respective electronic band gaps, which directly link the excitation efficiency of the defect (it is more difficult to move electrons across the gap when the HOMO and LUMO levels are far apart). However enhancing the optical properties of the nDs is possible by charging the defect or shifting the levels by other means. These points are indicated by our results for the NV − center at different positions. In the vacancy defects, the HOMO state is localized on the defect, while the LUMO state is more extended on the surface of the nD. In the Ni and PN defects, both HOMO and LUMO are localized around the defect center.
Overall, our aim was to show the variability in the electronic features of tiny defective nanodiamonds and the vast possibilities for tuning their properties. Many more aspects need to be further investigated, such as the interplay between symmetry, degeneracy in the electronic levels, and Jahn-Teller distortions would be of interest for many applications. Of further technological relevance would be the comparison of bulk diamond [70,71] to molecular-sized nDs and the influence of surface states or active sites on the energy-level splitting. Other factors include the amount of inherent (in the defect) or injected charge in the nD, the defect type, larger nD sizes variance or environmental conditions. In this respect, temperature, the interaction with other nanostructures or molecules, the presence of a solvent or the application of external fields can further enhance the potential of nanodiamonds by tailoring their properties and providing certain functionalities. In view of the high potential of defective nanodiamonds in novel applications, this is of extreme importance as it provides pathways for controlling the functionalities of materials.