The positively charged nitrogen vacancy (NV$+$) center 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$\u2212$) center, but no photoluminescence (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$\u2212$. We calculate the energy levels of the multielectron states under the 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.

## I. INTRODUCTION

The nitrogen-vacancy center 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 center 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 center is found, for example, in nanodiamonds in a solution or a substrate or for a diamond in bulk. The NV$\u2212$ is a promising candidate for singly addressable spin due to its excellent opto-mechanical properties and has demonstrated spin initialization and readout with 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 the absence of electronic spin.^{5} Furthermore, there is interest in studying alternate charge states of NV$\u2212$ as an efficient readout of a single spin state, long term classical data storage, and high resolution microscopy.^{6–8} However, these studies typically use NV$0$ as it can be created with photoionization and easily identified optically.

Optically, NV$0$ and NV$\u2212$ are very easy to see under typical experimental conditions via the optical absorption/emission spectrum due to a distinctive zero phonon line and phonon sideband. In comparison, the NV$+$ state is not thermodynamically favorable under these conditions. Specifically, the NV$+$ state is only stable for Fermi levels lower than approximately 1 eV above the valence band maximum (VBM), where 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 toward creating the NV$+$ state has successfully implemented both band bending, which increases the band energies 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 $\u223c0.5$ eV in theoretical estimates,^{9,14,15} so this technique cannot guarantee the existence of the NV$+$ charge state. More concretely, Pfender *et al.*^{5} have measured 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 experimental 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, it has been predicted that 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 Screyvogel *et al.*,^{13} it is theorized that 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 the difference of self consistent field energies ($\Delta $-SCF) with the supercell method that has found success with the vibrational lineshape of the PL of a defect in bulk but is only accurate under the assumption that the excited state has different symmetry from the ground state.^{17} Cluster calculations allow us to use time-dependent density functional theory (TD-DFT) to obtain accurate zero phonon lines (ZPLs) and transition dipoles,^{18–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 are non-zero. Using our recently developed frequency cutoff method,^{21} we predict the full PL spectrum. This suggests that NV$+$ cannot be considered “dark” but there might be another effect preventing the experimental observation of the optical signature.

## II. METHODOLOGY

Our methodology is based on the frequency cutoff approach, recently reported in Karim *et al.*^{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-neighbor carbon atoms in their bulk position and removing atoms from the outer layer to achieve C$3v$ symmetry. We protonate the outer carbons and ensure that our structure has consistent sp3 hybridization and only CH and CH$2$ 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 the NV$\u2212$ charge state,^{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$\u2212$ calculations. However, the NV$+$ is diamagnetic, and this allows for DFT to be calculated with the closed shell Kohn–Sham equation^{27} rather than unrestricted spin resolved calculations as with NV$\u2212$.

The excited state is simulated under linear response time-dependent DFT using the TURBOMOLE software^{28,29} 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 for all excited states. Furthermore, C$s$ symmetry was found to give accurate excited state properties for NV$\u2212$ in the harmonic approximation.^{30} This procedure was performed for four excited states: the two lowest A$\u2032$ and the two lowest A$\u2033$ states.

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, Sec. I.

Vibrational modes for the ground state are calculated with finite-difference under DFT with the SNF software.^{33} 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 vibrational modes unique to the nanodiamond but not bulk NV$+$can be found by applying the same constraints present in the bulk by calculating the relaxed excited state while constraining the outer layer of carbons in the nanodiamond, as detailed in Karim *et al.*^{21} This new excited state and the ground state vibrational modes are used in VIBES^{34} to compute the constrained PHR factors

The difference between constrained and unconstrained PHR factors allows us to separate the modes into regions. The low frequency region consists of modes that are suppressed by the constrained calculation. These typically oscillate the entire nanodiamond due to finite-size effects. After this is a region where the modes should be relatively unaffected (or enhanced) by constraint. Between these regions, we can define a cutoff frequency. We remove all unconstrained PHR factors below this cutoff frequency and use the resulting factors to generate a PL spectrum under the displaced harmonic oscillator model in a modified version of VIBES (see Karim *et al.*^{21} for details).

## III. RESULTS AND DISCUSSION

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$\u2212$ increases the electrostatic repulsion between nuclei, which results in a larger space around the vacancy. The conformal change in comparison to NV$\u2212$ will lead to different coupling of vibrational modes and a distinct vibronic spectrum.

The single electron orbitals are shown in Fig. 2(b). Similar to NV$\u2212$ shown in Fig. 2(a), there are clear conduction and valence bands and importantly, there are also intraband levels local to the defect center. These results are consistent with and extend previous calculations reported in the literature.^{16} The conduction bands shown are from virtual Kohn–Sham orbitals that 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 Fig. 2(c), the ground state is a singlet state with A$1$ symmetry. For the excited states within the bandgap, there are two singlet states with A$\u2032$ and A$\u2033$ symmetry as well as a triplet E state. The singlet states 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 Fig. 4. As TD-DFT is used, the associated transition dipoles are calculated as reported in Fig. 4.

To obtain all the transition energies and oscillator strengths, we analyze the combined relaxed ground and excited state geometries under the displaced harmonic oscillator model. The results, reported in Fig. 4, show that, in contrast to common conclusions reported in the literature, the transition dipole moments, $\mu $, are not zero and about a quarter of the value of NV$\u2212$, which implies that the NV$+$ is not a dark state. Furthermore, the ZPL for the NV$+$ states are lower in energy but under the phonon sideband of the NV$\u2212$. This could explain why, to date, experimentally photon emission from this state has not been observed.^{9–13}

We consider convergence with the nanodiamond size. Unfortunately, all energies are highly dependent on the shape of the nanodiamond, for instance, larger nanodiamond sizes form graphite on the surface.^{37} Figure 3 shows emission values for selected sizes where the cluster retained diamond geometry (majority of bond angles are $109.5\xb0$). The solid lines show the transition within the defect and can be seen to flatten after 121 carbon atoms with fluctuations less than 0.1 eV. The transition from the valence band (Fig. 3 dotted line) similarly flattens with fluctuation less than 0.2 eV. The fluctuations follow the same trend with the size, for instance increasing for the C$145$NH$100$ nanodiamond, which supports our theory that this variation is due to shape dependence. The transition energies may be overestimated due to quantum confinement effects; however, our largest nanodiamond has been previously shown to recreate the ZPL using TD-DFT for the NV$\u2212$charge state^{21} and has a diameter of 1.37 nm, which has been shown to be large enough for vertical transitions in NV$\u2212$ to converge^{38} and is larger than 1 nm diameter nanodiamonds where quantum confinement effects are significant.^{39} Similarly, while the transitions between defect and bands are more susceptible to quantum confinement effects, the vertical transitions between an NV$\u2212$ defect and the conduction band did not continue to decrease beyond 1.4 nm diameter sized nanodiamond.^{38} We are confident our calculations are sufficiently large to predict transition energies within error of shape dependence.

We now analyze the vibrational properties of the NV$+$ defect. The PHR factors from unconstrained relaxation are shown in Fig. 5(a) for A$\u2032$ and (b) for A$\u2033$, respectively. The PHR factors when the outer carbons are fixed to ground positions are shown in Figs. 5(c) and 5(d) for A$\u2032$ and A$\u2033$ transitions, respectively. The difference between unconstrained and constrained spectra is plotted in Figs. 5(e) and 5(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 the 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 the bulk. Eliminating the vibrational modes with frequencies below this cutoff region should return the bulk behavior. 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 of NV$\u2212$ consistent with experimental measurements.^{21}

The large peaks in the PHR factors of NV$\u2212$ have been shown to be directly responsible for the main features of the PL spectrum. We compare the PHR factors for NV$+$ reported in Fig. 5 to those reported for the NV$\u2212$ charge state.^{21} NV$\u2212$ 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$\u2212$ 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$\u2032$ and A$\u2033$ transitions, respectively. It can also be seen that there is non-negligible coupling at frequencies up to 175 meV, compared to NV$\u2212$ for which it becomes negligible above 100 meV.

The mode with the largest contribution in the 60 meV region has energy of 56.5 meV. This mode has $e$ symmetry and describes the oscillation of the carbon atoms with dangling bonds in the vacancy. This corresponds to the oscillation of C–C (across) in Fig. 1. For the 80 meV region, the mode with the largest contribution to both transitions occurs at 87.9 meV. It has $a1$ symmetry and consists of the nitrogen atom oscillating along the z-axis (as annotated in Fig. 1) with minimal movement from neighboring atoms in the defect. Similarly, the mode at 88.4 meV has $a1$ symmetry and consists of the nitrogen moving along the z-axis but now the three carbons with dangling bonds in the vacancy radially displace from the z-axis. This corresponds to the oscillation of C–N (across) and C–C (across) in Fig. 1. Finally, the largest contribution at higher energies is at 149 meV. This mode has $a1$ symmetry and involves oscillation of the neighboring carbons to the nitrogen radially outward corresponding to the oscillation of C–N in Fig. 1.

To generate PL spectra, we apply the cutoff at the center of the region (43.2 meV) to the unconstrained nanodiamond calculations. After a cutoff is applied, the total Huang–Rhys constants can be found by adding the squared PHR factors. The total HR constants are 2.524 and 2.030 for A$\u2032$ and A$\u2033$ transitions, respectively, compared to approximately 3 for NV$\u2212$.^{40} The modified PHR values are used in VIBES to get a vibrational lineshape, which we center on the ZPL from Fig. 4 to get the PL.

Figure 6 shows the PL spectra for the two transitions A$\u2032$ and A.” The intensity of A$\u2033$ is higher than A$\u2032$ as expected from the larger transition dipole moment as shown in Fig. 4, but both lower than NV$\u2212$by a factor of around 20. The photon emission rate is calculated as the Einstein A coefficient.^{35} The ZPL is at $1.60\xb10.1$ and $1.62\xb10.1$ eV for A$\u2032$ and A$\u2033$, respectively, which are both lower energy than $1.95$ eV of NV$\u2212$ and within the phonon sideband of NV$\u2212$. Experimentally, the emission spectrum should match the A$\u2032$ curve (dashed red line) as it is more energetically favorable for the excited state to deform to the A$\u2032$ configuration (longer wavelength).

It has been proposed in the previous literature that the defect charge transition from NV$0$ to NV$+$ at 1 eV above the VBM indicates photoionization may occur instead of bright transitions of higher energy.^{16} The defect charge transition state is calculated by assuming that 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 VBM. 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 population 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.^{41} There exist both radiative and non-radiative transitions between charge states. Non-radiative transitions consist of electrons tunneling between defect centers 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$\u2212$ radiative photoionization transition energies has recently been performed by Razinkovas *et al.*^{42} using defect to band transitions. They calculate the energy difference between NV$\u2212$ and NV$0$ with one electron excited to the conduction band. They found the energy difference of the ground state of NV$\u2212$ to NV$0$ with one CBM electron gave the experimental threshold for photoionization 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 Fig. 2(c). We perform this calculation in C$s$ geometry, similar to the transition to the $1$E state; however, now we optimize for the second excited A$\u2032$ and A$\u2033$ states. The energies are 3.09 and 3.14 eV, respectively, which are higher than our proposed bright ZPL of 1.6 eV. Illumination by 1.6 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 $|a1\u27e9$ to $|e\u27e9$, and then from the valence band to $|a1\u27e9$ as shown in Fig. 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$\u2212$/NV$0$ transition.^{41} Therefore, it should be theoretically possible to see PL without resulting in immediate photoionization. For further details, see supplementary material, Sec. II.

As NV$+$ measurements often occur under the electric field,^{5,11–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 photoionization is small.

. | Gap (eV) . | Ground energy (eV) . | ||
---|---|---|---|---|

Electric field (kV/cm) . | Value . | Diff . | Value . | Diff . |

−100 | 2.4645 | 0.0071 | −207 661.967 | −0.004 |

100 | 2.4787 | 0.0071 | −207 661.975 | 0.004 |

0 | 2.4716 | … | −207 661.971 | 0 |

. | Gap (eV) . | Ground energy (eV) . | ||
---|---|---|---|---|

Electric field (kV/cm) . | Value . | Diff . | Value . | Diff . |

−100 | 2.4645 | 0.0071 | −207 661.967 | −0.004 |

100 | 2.4787 | 0.0071 | −207 661.975 | 0.004 |

0 | 2.4716 | … | −207 661.971 | 0 |

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 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 contradiction with recent experimental observations of NV$+$ concluding that it is a dark state, this could be accounted for by a charge dependent quenching mechanism.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the comparison of transition energies under the Tamm–Dancoff approximation for nanodiamonds of sizes C$121$NH$100$, C$145$NH$100$, and C$197$NH$140$ and for a detailed explanation of our proposed transition between NV$+$ and NV$0$.

## ACKNOWLEDGMENTS

We thank Cristian Bonato, Alastair Stacey, and Brett Johnson for useful discussions. A.P. acknowledges RMIT University Vice-Chancellor’s Senior Research Fellowship, ARC DECRA Fellowship (No. DE140101700), and Google Faculty Research Award. This work was supported by the Australian Government through the Australian Research Council under the Centre of Excellence scheme (Nos. CE170100012 and CE170100026). It was also supported by computational resources provided by the Australian Government through the National Computational Infrastructure Facility and the Pawsey Supercomputer Centre.

## DATA AVAILABILITY

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.