Bright $\mathrm{\textit{ab-initio}}$ photoluminescence of NV+ in diamond

The positively charged nitrogen vacancy (NV+) centre in diamond has been traditionally treated as a dark state due to the experimental lack of an optical signature. Recent computational studies have shown that it is possible for the NV+ defect to have an excited state transition equivalent to that of the negatively charged (NV-) centre, but no PL predictions have been reported so far. We report the first $\mathrm{\textit{ab-initio}}$ calculation showing that the NV+ center presents quantum emission, with zero phonon line at 765 nm and a non-zero transition dipole moment, approximately 4x smaller than the transition dipole moment of NV-. We calculate the energy levels of the multielectron states under time-dependent density functional theory (singlet and triplet E states), and using our recently developed frequency cutoff method, we predict the full PL spectrum. Our results suggest that this state cannot be considered intrinsically 'dark' and charge specific quenching mechanisms should be investigated as the cause of the lack of optical activity in experimental characterizations.

The positively charged nitrogen vacancy (NV + ) centre in diamond has been traditionally treated as a dark state due to the experimental lack of an optical signature. Recent computational studies have shown that it is possible for the NV + defect to have an excited state transition equivalent to that of the negatively charged (NV − ) centre, but no PL predictions have been reported so far. We report the first ab-initio calculation showing that the NV + center presents quantum emission, with zero phonon line at 765 nm and a non-zero transition dipole moment, approximately one quarter of the transition dipole moment of NV − . We calculate the energy levels of the multielectron states under time-dependent density functional theory (singlet and triplet E states), and, using our recently developed frequency cutoff method, we predict the full PL spectrum. Our results suggest that this state cannot be considered intrinsically 'dark' and charge specific quenching mechanisms should be investigated as the cause of the lack of optical activity in experimental characterizations.
The nitrogen-vacancy centre in diamond is an attractive platform for quantum information processing [1], quantum sensing [2] and for storing quantum information [3], and has been widely studied theoretically [4].
The study of the NV centre in diamond has been dominated by the negatively and neutrally (NV 0 ) charged states. This is due to the relative stability these charge states have in typical situations where the NV centre is found, for example in nanodiamonds in solution or substrate, or for diamond in bulk. The NV − is a promising candidate for singly addressable spin due to its excellent opto-mechanical properties and has demonstrated spin initialization and readout with a long spin coherence time.
Recently, the NV + charged state has attracted interest as a potential way to access long-lived nuclear spin states for spin coherence storage due to its absence of electronic spin [5]. Furthermore, there is interest in studying alternate charge states of NV − as efficient readout of a single spin state, long term classical data storage, and high resolution microscopy [6][7][8]. However, these studies typically use NV 0 as it can be created with photoionization and easily identified optically.
Optically, NV 0 and NV − are very easy to see under typical experimental conditions via the optical absorption/emission spectrum due to a distinctive zero phonon line and phonon side band. In comparison, the NV + state is not thermodynamically favourable. Specifically, the NV + state is only stable for Fermi levels lower than approximately 1 eV above the VBM, however N-doped diamond has a Fermi level of around 2 eV above the valence band maximum (VBM) [9]. In order to fabricate stable NV + experimentally, it is necessary to reduce the Fermi level.
Experimental work towards creating the NV + state has successfully implemented both band bending, which increases the bands and relatively reduces the Fermi level, with functionalization of the diamond surface chemically [10], and directly controlling the Fermi level with in plane gates [5,11], as well as a combination of both [12,13]. These experimental papers reached a Fermi level that theoretically should stabilize NV + but only detected background optical emission. This is problematic, as the exact charge transition level varies by ∼0.5 eV in theoretical estimates [9,12,14,15], so this technique cannot guarantee the existence of the NV + charge state. More concretely, Pfender et al. [5] have measured an increased amplitude Rabi Oscillations of the nuclear spin, consistent with computational predictions that NV + should have no total electron spin [16].
The lack of an optical signature makes the NV + difficult to identify and characterize. While hyperfine constants have been calculated [16] suggesting potential avenues of magnetic identification, no current optical identification is known and this state is often considered dark due to the lack of optical activity. This lack of optical activity is considered characteristic of the NV + state, however there exist optical transitions in the NV + that are dipole and spin allowed [16]. The experimentally observed dark state suggests there may be external quenching of the photoluminescence (PL) due to non-radiative transitions [5,16], and in [13] it is theorised that the PL is quenched by proximity to the two-dimensional hole accumulation layer created by functionalizing the surface that quenches NV + selectively.
Calculating photoluminescence is generally limited by the computational complexity of quantum chemistry. A solid state approach uses Delta self consistent field (∆-SCF) with the supercell method which has found success with the vibrational lineshape of the PL of a defect arXiv:2104.04256v2 [cond-mat.mes-hall] 25 Aug 2021 in bulk, but is only theoretically valid for excited states with different symmetry from the ground state [17]. Cluster calculations allow us to use time-dependent density functional theory (TD-DFT) to obtain accurate ZPL and transition dipoles [18][19][20] but the result does not apply to solid state systems. Recently, we have demonstrated a method based on cluster calculation with frequency cutoff that can accurately predict the ZPL and phonon sidebands of solid state systems [21].
In this work we report the first ab-initio photoluminescence calculation showing that the NV + state is optically active and has a characteristic emission and absorption spectra.
We perform TD-DFT calculation of the energy levels of the multielectron states (singlet and triplet states), and find that the transition dipole moments of the lowest singlet-singlet transition is non-zero. Using our recently developed frequency cutoff method [21], we predict the full PL spectrum. This suggest that this state cannot be considered 'dark' but there might be another effect preventing experimental observation of the optical signature.
Our methodology is based on the frequency cutoff approach, recently reported in [21]. We use their nanodiamond shape, exchange-correlation functional and basis set. In brief, we calculate the photoluminescence (PL) spectra of a defect in bulk by applying a frequency cutoff to cluster calculations. The PL of the cluster is given under the displaced harmonic oscillator model. The only inputs for this model are the relaxed ground state, found with DFT; the relaxed excited state, found with TD-DFT; and the normal modes given by finite-difference under DFT of the ground state. This gives us the PL of the defect in a nanodiamond. To recover the solid state spectrum, this is followed by a TD-DFT calculation while constraining the outer layer of the cluster, which allows us to identify normal modes that occur on the surface, unique to the cluster. We eliminate this region of normal modes with a cutoff to recover the bulk PL spectrum.
To study the NV + center, we construct a 1 nm diameter nanodiamond by adding to the vacancy the four nearest-neighbour carbon atoms in their bulk position and removing atoms from the outer layer to achieve C 3v symmetry. We protonate the outer carbons to have consistent sp3 hybridisation and only CH and CH2 functional groups, which gives the chemical formula: C 197 NH 140 .
Next, the ground state is relaxed from the bulk positions and the electron density and associated ground state energy under DFT in TURBOMOLE [22,23] are converged for the stationary atoms. Throughout the process the symmetry is fixed to C 3v . Consistent with previous calculations with NV − [21], we use def2-SV(P) [24] as a basis set with PBE0 [25,26] as the exchange-correlation functional. This allows explicit comparisons with NV − calculations.
FIG. 1. The ground state structure of NV + . Carbon atoms are shown in grey, and nitrogen atoms in blue (only the atoms around the defect are shown). The table reports the distance inÅ for atoms in the defect for NV + and NV − . Due to the reduced negative charge in the defect, the positive charged state has larger separations around the vacancy compared to the negatively charged state. Z axis for C3v geometry is shown and corresponds to the direction from the centre of the vacancy to the nitrogen atom. Figure 1 shows the ground state structure of NV + as well as the bond lengths and distances between the relevant atoms in the vacancy. For NV + , the lack of the electron and consequent decreased negative charge within the vacancy compared to NV − increases the electrostatic repulsion between nuclei, which results in a larger space around the vacancy. The conformal change in comparison to NV − will lead to different coupling of vibrational modes, and a distinct vibronic spectrum.
The excited state is simulated under Linear Response Time-Dependent DFT using the TURBOMOLE software [27,28] with the same basis set and functional. For the excited states, the symmetry is fixed to C s . NV + was initially relaxed under C 1 symmetry and found to adopt C s symmetry. Furthermore C s symmetry was found to give accurate excited state properties for NV − in the harmonic approximation [29].
As the NV + is diamagnetic, all electrons are paired by spin. This allows for DFT to be calculated with the closed shell Kohn-Sham equation [30] rather than unrestricted spin resolved calculations as with NV − . The single electron orbitals are shown in Figure 2(a). Similar to NV − shown in Figure 2 Kohn-Sham single electron orbital energies for NV − (a) and NV + (b). The multielectron energy ladder for NV + singlet (left) and triplet (right) are shown in c). Part a) shows alpha (left) and beta (right) electron orbitals separately since electrons are unpaired in NV −Ḟ ishhooks represent electron occupation. All multielectron states below 2 eV are intraband transitions; all higher states up to 4 eV are transitions from the valence band into the |e state. Our proposed bright transition is shown as a black arrow in b) and c). The lowest singlet valence band transition occurs at 3.09 eV and is shown as a cyan arrow. It should be noted the virtual Kohn-Sham orbitals do not correspond to actual energies if they were occupied, so the bandgap cannot be read from a) or b) and excitation energies must be read from c). bands shown are from virtual Kohn-Sham orbitals which may not accurately reflect the actual conduction band energies. The transition energies are calculated in the multielectron description.
We calculate the multielectron state energy ladder and report it in Figure 2(c), the ground state is a singlet state with A 1 symmetry. For the excited states within the bandgap, there are two singlets states with A' and A" symmetry as well as a triplet E state. The singlet states are degenerate on absorption (with E symmetry), but have different adiabatic transitions due to the reduction of symmetry from C 3v to C s in the excited state geometry. All higher excited states involve exciting electrons from the valence band to the empty e levels. The energy levels are reported in Figure 4. As TD-DFT is used, the associated transition dipoles are calculated as reported in Figure 4.
To access the singlet-triplet transition, it is necessary to use the Tamm-Dancoff Approximation (TDA) to avoid instabilities. The TDA ignores contribution by "de-excitation" terms i.e. the coupling from virtual to occupied orbitals, and has been shown to be valid for determining energies of single-excitation states, particularly for optical spectra [31,32]. Comparison of TDA transition energies for the singlet state are shown in supplementary material section I.
To obtain all the transition energies and strengths, we analyze the combined relaxed ground and excited states under the displaced harmonic oscillator model. The results, reported in Figure 4, show that, in contrast to common conclusions reported in the literature, the transition dipole moments, µ, are not zero and about a quarter of the value of NV − , which implies that the NV + is not a dark state. Furthermore, the zero phonon lines (ZPL) for the NV + states are lower in energy but under the phonon sideband of the NV − . This could explain why, to date, experimentally photon emission from this state has not been observed [9][10][11][12][13].
We consider convergence with nanodiamond size. Unfortunately, all energetics are highly dependent on the shape of the nanodiamond and, furthermore, larger nanodiamond sizes can form graphite on the surface [33]. Figure 3 shows emission values for select sizes where the cluster retained diamond geometry (majority of bond angles are 109.7 • ). The ZPL can be seen to flatten after 121 carbon atoms and oscillate by around 0.1 eV around their central values.
We now analyse the vibrational properties of the NV + defect. Vibrational modes for the ground state are calculated with finite-difference under DFT with SNF software [34]. We can then calculate the Partial Huang-Rhys (PHR) factors for each vibrational mode and each transition, which show how each vibrational mode couples with a transition. The results are shown in Figure 5 The difference between unconstrained and constrained spectra is plotted in Figure 5(e) and (f). Here we can see that at low frequencies (on the left of the green shaded area) the presence of large positive peaks indicate that the coupling to vibrational modes are present in the nanodiamond but not in bulk, whereas at high frequencies (on the right of the green shaded area) the vibrational modes equally couple to the transition for both bulk and nanodiamond. The green shaded area is therefore identified as a region separating the unique properties of the nanodiamond from those in common with bulk. Eliminating the vibrational modes with frequencies below this cutoff region should return the bulk behaviour. The lower boundary of the cutoff region is at 36.8 meV, and the upper boundary is at 49.5 meV. Choosing a cutoff value within this region has been shown to reproduce a simulated PL consistent with experimental measurements [21].
The large peaks in the PHR factors of NV − have been shown to be directly responsible for the main features of the PL spectrum. We compare the PHR factors for NV + reported in Figure 5 to those reported for NV − in [21]. NV − exhibits a large peak in the low frequency region, at 30 meV which is reduced under constrained calculation. Similarly the NV + has a large peak at 25 meV. In the higher frequencies, NV − has a large peak unaffected by constraining, centered at 66 meV, whereas NV + has two peaks, one at 60 meV and one at 80 meV associated with the A' and A" transitions respectively. It can also be seen that there is non-negligible coupling at frequencies up to 175 meV, compared to NV − for which it becomes negligible above 100 meV.
To generate the PL spectra, we apply the cutoff at the center of the region (43.2 meV) to the unconstrained nanodiamond calculations. The modified PHR values are used in VIBES to get a vibrational lineshape, which we center on the ZPL from Figure 4 to get the PL. Figure 6 shows the PL spectra for the two transitions A' and A". The intensity of A" is higher than A' as expected from the larger transition dipole moment as shown in Figure 4, but both lower than NV − by a factor of around twenty. The photon emission rate is calculated as the Einstein A coefficient [35]. The ZPL is at 1.60±0.1 and 1.62 ± 0.1 eV for A' and A" respectively, which are both lower energy than 1.95 eV of NV − and within the phonon sideband of NV − . Experimentally, the emission spectrum should match the A' curve (dashed red line) as it is more energetically favourable for the excited state to deform to the A' configuration (longer wavelength). It has been proposed in previous literature that the defect charge transition from NV 0 to NV + at 1 eV above the VBM indicates photoionisation may occur instead of bright transitions of higher energy [16]. The defect charge transition state is calculated by assuming the defect is coupled to an electron bath with some chemical potential. The transition state is defined as the chemical potential of the bath where the defect formation energy for two charge states are identical. This corresponds to the energy difference between the two charge states when the chemical potential of the bath is at the CBM. The defect transition energies can be considered the chemical potential of the bath required for 50% population of each charge state. This model adequately explains thermodynamics of charge states. Experimentally, when the Fermi level (the chemical potential of the bath) is tuned to cross the defect transition energies, the popula-tion of one charge state changes to the next as indicated by change in photoluminescence [5].
However, investigation into NV + /NV 0 extrinsic optical dynamics shows radiative transitions that cannot be explained under this model [37]. There exist both radiative and non-radiative transitions between charge states. Non-radiative transitions consist of electrons tunnelling between defect centres and electron donors/acceptors and can be controlled by isolating the defect. Radiative transitions involve electrons entering or exiting the defect via the bulk bands and require explicit calculation of defect to band transitions. Explicit calculation of NV 0 /NV − radiative photoionisation transition energies has recently been performed by Razinkovas [38] using defect to band transitions. They calculate the energy difference between NV − and NV 0 with one electron excited to the conduction band. They found the energy difference of the ground state of NV − to NV 0 with one CBM electron gave the experimental threshold for photoionisation of 2.64 eV.
Our equivalent transition is between the ground state  of NV + with NV 0 with one fewer electron in the valence band. Under TD-DFT we can simulate this by relaxing the excited state with one electron excited from the valence band, which is shown by the cyan arrow in Figure 2(c). We perform this calculation in C s geometry, Similar to the transition to the 1 E state, however now we optimise for the second excited A' and A" states. The energies are 3.09 and 3.14 eV respectively, which are higher than our proposed bright ZPL of 1.7 eV. Illumination by 1.7 eV light should be able to drive the bright transition within the defect without inducing this charge conversion transition.
There is also a transition from the excited 1 E state to a possible NV 0 state, by first exciting from |a 1 to |e , and then from the valence band to |a 1 as shown in Figure 2(b). This transition has an energy of 1.3 eV and may result in charge state conversion. However, since it requires two excitations, it will scale second order with power, and the timescales will be much longer than transitions within a charge state as is the case for the equivalent second order NV − /NV 0 transition [37]. Therefore, it should be theoretically possible to see PL without resulting in immediate photoionisation. For further details see supplementary information section II.
As NV + measurements often occur under electric field [5,[11][12][13], we consider the effect of the electric field on the energy of the defect. Hauf et al. [11] experimentally observed NV + under an electric field strength of 100 kV/cm. To understand the effect of the electric field on the ground state energy, we apply a linear electric field and perform a relaxation under DFT in TURBOMOLE, to get the new ground state energy and associated deformation for those values. The calculation results are reported in Table I, where we see that applying an electric field of 100 kV/cm shifts the homo-lumo gap by only 7 meV, which suggests that the expected effect on the PL spectra and photoionisation is small.
We have demonstrated that the NV + state contains an optically active transition between its singlet ground and singlet excited levels. We have studied the emission both electronically and vibrationally, calculating the TD-DFT energy levels for the excited triplet state and found that the dipole transition moment is non-zero. Finally, we have predicted the characteristic photoluminescence spectra and shown that the effects of the electric field are negligible. These results are consistent across three sizes of nanodiamond. While our results appear to be in con-tradiction with recent experimental observations of NV + concluding that it is a dark state, this could be accounted for by a charge dependent quenching mechanism.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

SUPPLEMENTARY MATERIAL I: Tamm-Dancoff Approximation
The entire methodology was performed for NV + in nanodiamonds of three different sizes: C 121 NH 100 , C 145 NH 100 and C 197 NH 140 . In the main text we showed only results for the largest nanodiamond, and in this section we report values for the other sizes. Table S1 shows the transition energies and strengths for all three sizes of nanodiamond. It can be observed that the Tamm-Dancoff Approximation (TDA) overestimates the singlet transition energies by 10% for all sizes. The dipole moments of C 121 NH 100 and C 145 NH 100 diamond are overestimated by an order of magnitude. For singlet states, we use TD-DFT without TDA, however for the triplet states we use TDA due to triplet instabilities. The singlettriplet transition energy is 0.79 ± 0.08 eV, where we estimate the error, from the singlet-singlet result, to be about 10%. Singlet-triplet transitions are dipole forbidden, so there is no need to calculate the transition dipole moments for the triplet state.
There is a second triplet state transition at 2.80 eV that corresponds to exciting from the valence band into the defect in a triplet E state. There is no radiative transition from the ground state into this state. The transition between triplet ground state and this valence band state is approximately 2 eV. Figure S1 shows the multielectron energy levels and subsequent charge state transition between NV − , NV 0 , and NV + . There are intrinsic dynamics within a charge state and extrinsic dynamics that allow transitions between charge states [37]. Non-radiative transitions involve tunnelling of an electron between the NV centre and a nearby well, typically an N defect or a functional group on the surface and vary according to proximity to the defect [37]. These are shown as dotted arrows in Figure S1. Radiative transitions, however, involve excitations of the bulk diamond band structure and are due to the defect itself. Assuming isolation of the defect, under illumination, charge transitions will occur primarily through radiative transitions.

SUPPLEMENTARY MATERIAL II: Charge State Conversion
A radiative excitation from NV − to NV 0 occurs through the conduction band. An electron from NV − can be excited into the conduction band leaving the defect in the NV 0 charge state. The electron in the conduction band can then relax into another well such as a defect or surface functional group [37]. The energy to excite an electron from the ground state is 2.6 eV and can be seen experimentally [37]. These are shown as solid arrows in Figure S1.
We propose a similar mechanism but involving the valence band for the transitions from NV + to NV 0 . In our case, we assume an isolated NV + defect with no electrons in the conduction band under illumination. The optical transitions available are given as solid arrows in  S1. Diagram of transitions between the charge states of the NV centre in diamond. The transitions between the NV 0 and NV − , shown inside the black border, are reproduced from Doherty et al. [37] with permission from Elsevier. The radiative charge state transition from NV − to NV 0 occurs first through excitation to the conduction band. In addition, we propose that for NV + , assuming no electrons in the conduction band, the only available radiative mechanism occurs through an excitation from the valence band ("V.B."). Figure S1. We propose that under illumination, an electron can be promoted from the valence band into the defect. The hole in the valence band can then be filled elsewhere in the crystal similar to how the conduction band electron is assumed to relax elsewhere in the NV − to NV 0 transition [37]. Our proposed intrinsic dynamics involve a bright transition between the ground 1 A 1 state and 1 E singlet states at 1.7 eV. We also see a triplet 3 E state in-between the singlet states.
We also propose bright transitions where a valence band electron is promoted to a defect level. These are shown as a grey box labelled "V.B." in Figure S1. For more details, see main text Figure 2(b) and 2(c). The transition from the ground state to the lowest band state is 3.08 eV. We propose this corresponds to the equivalent NV − transition at 2.6 eV.
The second order transition requires exciting the defect to the singlet E state before exciting into the valence band. The second excitation should take 1.38 eV. It is possible for this transition to occur under illumination of 1.7 eV light, however, since it is second order, it should occur at a longer timescale than the intrinsic dynamics [37] and, more importantly, scale as the square of input power. Therefore, this effect will be suppressed at low power. Calculating the transition dipole moment between two excited states is outside of the scope of this paper, so the actual rates cannot currently be compared.
If the defect is in its triplet ground state, there is also an optical transition of around 2.0 eV. This, similarly, will occur at long timescales since NV + must change spin states.