Inverse-designed photon extractors for optically addressable defect qubits

Solid-state defect qubit systems with spin-photon interfaces show great promise for quantum information and metrology applications. Photon collection efficiency, however, presents a major challenge for defect qubits in high refractive index host materials. Inverse-design optimization of photonic devices enables unprecedented flexibility in tailoring critical parameters of a spin-photon interface including spectral response, photon polarization and collection mode. Further, the design process can incorporate additional constraints, such as fabrication tolerance and material processing limitations. Here we design and demonstrate a compact hybrid gallium phosphide on diamond inverse-design planar dielectric structure coupled to single near-surface nitrogen-vacancy centers formed by implantation and annealing. We observe device operation near the theoretical limit and measure up to a 14-fold broadband enhancement in photon extraction efficiency. We expect that such inverse-designed devices will enable realization of scalable arrays of single-photon emitters, rapid characterization of new quantum emitters, sensing and efficient heralded entanglement schemes.


Introduction
Optically addressable defect qubits in materials like diamond 1, 2 and silicon carbide 3 are promising for the realization of a wide range of quantum technologies including singlephoton generation, quantum metrology, and quantum information protocols 4,5 . In all of these applications, photon collection efficiency is a key figure of merit. The photoluminescence (PL) detection efficiency is intrinsically limited by the non-directional emission of PL and total internal reflection between the high-index host material and its low-index environment. To enable high PL collection efficiency, photonic structures such as solid-immersion lenses 6 , waveguides 7, 8 and microcavities [9][10][11][12] have been utilized. Nanophotonic devices are particularly attractive for their small footprint and potential for scalable integration. However, optimizing nanophotonic structures for efficient defect integration is often nontrivial and bespoke due to the unique constraints imposed by the targeted application for a specific defect system. For high-sensitivity quantum metrology with nitrogen vacancy (NV) centers in diamond, photonic structures should be optimized for efficient broadband extraction of PL from NVs located a few nm from the sensing surface. In contrast, for quantum information applications, photonic coupling to the sharp zero-phonon line (ZPL) of deep NV centers (> 100 nm from the surface) and operation under the high-cooperativity (Purcell-enhanced) regime is preferable.
Here we present a flexible inverse-design optimization framework that can reconcile a wide range of design constraints and generate planar dielectric gallium phosphide(GaP)-on-diamond photonic structures. For a dipole located 100 nm from the surface and oriented perpendicular to the NV axis, an optimized 1.5 µm × 1.5 µm device is calculated to provide a 17-fold average PL enhancement of the free-space PL collection. Recent work 13 suggests that the performance of such an inverse-designed extractor is close to the theoretical limit. Orders of magnitude larger enhancement factors are found for defects positioned closer to the surface. Additionally, during optimization, we constrain patterning to the dielectric layer in order to minimize perturbations of the defect environment.
Experimentally, we fabricate the optimized GaP-ondiamond photon extractors and observe efficient PL collection from shallow (≈ 100 nm deep) single NV centers created by nitrogen ion implantation and vacuum annealing. These versatile devices exhibit a broadband PL enhancement for wavelengths in the measured range of 575 nm to 750 nm. We observe up to 14-fold enhancement of free-space PL collection for device-coupled single NVs. Extensive optical characterization of the NV centers performed both before and after fabrication provides insight into changes in the local NV environment. Post-fabrication sample treatment substantially improves device NV stability and points to important fabrication/design considerations for future defect qubit-photonics integration.

Design of inverse-designed photon extractors
The photon extractors are designed for robust enhancement of ZPL photon collection from near-surface NV centers, and modeled via a 3D finite-difference frequency-domain (FDFD) solver, with the frequency set at the negativelycharged NV − ZPL of 637 nm. Topology optimization 14 was used to design the device within a design region situated directly on top of the diamond substrate with dimensions 1.5 µm × 1.5 µm × 0.25 µm. A harmonic dipole source representing the NV center is situated 100 nm below the diamond- The optimization objective for any given dipole polarization is the net Poynting flux through a square collection surface ∂ F of sidelength 1.5 µm, situated 0.4 µm above the top surface of the device. To make the device robust with respect to uncertainty in the NV center polarization state, a minimax formulation was used so as to maximize the minimum Poynting flux contribution from dipoles withx,ŷ, andẑ polarizations.
The optimization problem is formulated as follows: where S = 1 2 Re{E * × H} is the Poynting vector, ε i the permittivity at the ith location within the design region, and the superscript j ∈ [x, y, z] runs over different polarizations of the Figure 2. a. Illustration of material stack and measurement configuration. Implanted NV centers coupled to grating devices are excited with off-resonant 532 nm (resonant 637nm) laser source and the ZPL at 637 nm (PSB emission from 650 to 800 nm) is collected for photon-antibunching (PLE) studies. b. Comparison of simulated (grey) and measured (blue) broadband enhancement for inverse extractor devices coupled to NV ensemble (Sample A). The red line indicates the average simulated spectrum and the red envelope is one standard deviation around the average. Variation in individual spectra for both simulation and measurement stems from randomness of NV positions relative to any given device on top of the ensemble sample. c. Spectra showing broadband NV PL enhancement for NV ensemble (Sample A) d. Spectra comparing ZPL emission of device-coupled and a non-device (bare) single NV (Sample B). e. Saturation counts from devices coupled to single NVs. Inset: Histogram of measured devices and the observed ZPL enhancement (Sample B). f. g 2 (τ) measurement of device coupled to single NV under 532 nm 300 µW excitation (Sample B). time harmonic dipole source. t is an auxiliary variable representing the minimum of the flux for the different source polarizations, introduced so the optimization objective and constraints are all differentiable with respect to the degree of freedom. In principle, topology optimization allows the relative permittivity of each spatial discretization voxel within the design region to be a degree of freedom; to accommodate fabrication using electron beam lithography, the actual degrees of freedom form a 2D grid and represent a top-down view of the device. Density and binarization filters 15 were used to restrict the minimum feature size to around 50 nm, well within the capabilities of electron beam lithography. The adjoint method 16,17 is used to efficiently calculate the gradient for the design objective with respect to every degree of freedom.
A schematic representation and xy cross-section of the final design are shown in Fig. 1a. Fig. 1b shows the LDOS and collection flux enhancement spectra of the device for different dipole orientations, withx andŷ orientations producing nearly the same spectrum due to the resulting near mirror device symmetry. A few salient features of the optimized design are worth commenting upon. To begin with, the extractor yields orders of magnitude larger flux enhancements forẑ compared to in-plane polarized dipoles. Such a vast difference in relative improvement follows from the fact that the bulk of the radiation from aẑ polarized dipole underneath a bare diamond interface undergoes total internal reflection, suppressing its radiation by more than a order of magnitude compared to that of an in-plane polarized dipole; the main role of the extractor is therefore to out-couple such non-radiative modes.
Further, device performance is robust against spatial displacement of the NV center; the enhancement factor decreases by less than 50% for spatial displacements smaller than 70 nm from the original location of the NV center (Supplemental Information SI.6).
It is noted that while the extractor leads to large enhance-3/8 ments in the collected flux, it does not significantly increase the net radiation from the NV center, quantified by the local density of states (LDOS); hence, most of the observed improvement comes from a higher collection efficiency as opposed to fluorescence rate enhancement via the Purcell effect. Fig. 1d, e explores this apparent trade-off by comparing the relative LDOS and collection enhancements of inverse-designed devices targetingẑ polarized dipoles at various depths. For shallow dipoles, optimizing for larger LDOS also significantly enhances the collected flux, though to a lesser degree than what is achieved by directly optimizing for collected flux. However, the disconnect between these two design objectives is observed to grow rapidly with increasing dipole depth: beyond a depth of 75 nm, devices designed to maximize the collected flux maintain a flux enhancement ratio of about 100, while their LDOS enhancement is close to one. Notably, beyond a dipole depth of 150 nm, even LDOS optimized devices are only able to achieve LDOS enhancements of roughly two, reflecting the rapidly decaying access of the device to the dipole's near field. Devices that combine the benefits of Purcell enhancement with increased collection efficiency through interference are therefore only needed for NV centers sufficiently close ( 100 nm) to the interface.

GaP photon extractors coupled to implanted NVs
We fabricate photon extractors (Fig. 1c) on two samples, one for ensemble-averaged measurements (Sample A) and the other for single-emitter characterization (Sample B). Sample A is a high pressure high temperature synthesised diamond (Element Six, N < 200 ppm, B < 0.1 ppm) implanted with 14 N accelerated to 20 keV and vacuum annealed at 800 • C. Here, annealing forms NVs primarily by vacancy recombination with native nitrogen. The resulting NV distribution is dictated by the vacancy diffusion profile 18 , yielding a dense layer (≈100 NVs per 800 nm excitation spot diameter) of NV centers ≤100 nm below the diamond surface.
Sample B is a chemical vapor deposition diamond (Element Six, electronic grade, N < 5 ppb, B < 1 ppb) implanted with 15 N accelerated to 85 keV and vacuum annealed at 1200 • C (see Methods). During annealing, NV centers are formed primarily by vacancy recombination with the implanted nitrogen, yielding a thin layer (≈ 3 NVs per excitation spot) of single NVs 100 nm ± 20 nm from the surface.
After NV-formation, a 250 nm thick gallium phosphide (GaP) membrane is transferred to each sample via a wet lift-off process 11,12,19 . Electron beam lithography and subsequent plasma RIE etching of the GaP layer forms the 1.5 µm × 1.5 µm photon extractors. The small footprint was chosen for compatibility with on-chip electrodes enabling optical frequency tuning 12,20 . Over 100 000 devices are fabricated in multiple arrays on the 2 mm × 2 mm diamond substrate. An array of fabricated devices and a false color SEM image highlighting the material stack is shown in Fig. 1c. On average, each feature is measured to be within 10% of the design, with near vertical (θ = 88 • ) GaP sidewalls. See Methods, Supplemental Information SI.2 for more details on fabrication.

Enhanced photon extraction
To study NV PL enhancement from the fabricated devices, we use spatially-resolved photoluminescence spectroscopy under off-resonant 532 nm excitation (Fig. 2a). First, sample A is utilized for broadband characterization of ensemble-averaged enhancement over the entire NV spectrum. The room temperature spectra from devices are normalized to a non-device ensemble on the same sample. We observe an average sixfold enhancement that is relatively flat over 575 to 750 nm and matches well to the theoretical spectrum-averaged four-fold enhancement (Fig. 2b, c).
Having demonstrated broadband NV PL enhancement, we move on to single-NV coupled devices on sample B and characterize the enhancement of the NV zero-phonon line (ZPL) emission. Identification of single-NV coupled devices is performed by comparing the low temperature (8 K) NV − ZPL spectra under 532 nm excitation of devices to nearby nondevice NVs (Fig. 2c). Enhanced ZPL collection rates ( Fig. 2d) were exhibited in 74 of 480 tested devices. The distribution of observed enhancement can be seen in Fig. 2e. For devices on sample B, the placement of NVs with respect to the optical mode is random; coupled device yield can be significantly improved by targeted implantation or pattern alignment to registered defect centers 21 during fabrication.
We verify that the observed enhanced PL emission corresponds to a single-emitter using an autocorrelation measurement on the NV − ZPL at 637 nm. Fig. 2f shows the normalized coincidences under continuous off-resonant excitation for a device-coupled NV. The dip in coincidences at 0 time delay is 0.12 < 0.5 and verifies single-photon emission. Autocorrelation measurements on the eight brightest devices confirm single-NV coupling in six of the devices. In addition to the dip at zero time delay, significant bunching is observed in all devices on the 100 ns timescale. In modeling the autocorrelation curve (red line, Fig. 2f), we find it is necessary to include both the NV − singlet state and the NV 0 charge state to reproduce the magnitude and timescale of the bunching for a series of power-dependent measurements. Consideration of charge state dynamics is a critical for NV device performance as discussed further below. Consistent with theoretical expectations (Fig. 1b, dashed lines), excited state lifetime reduction (due to an increased LDOS) is not necessary to obtain good agreement between the experimental data and model. Details of the modelling procedure and optimized rates are provided in the Supplemental Information SI.5.

Charge state and spectral stability
Single defect qubit devices have requirements beyond photon collection efficiency including charge state stability and spectral homogeneity and stability. These properties should be preserved during device integration. Our autocorrelation model suggests that rapid charge-state conversion occurs between the neutral (NV 0 ) and negatively-charged (NV − ) states. For both sensing and quantum information applications, the NV − charge state is required, hence we need to minimize ionization into the NV 0 charge state. Every NV center, device or non-device, has emission at both NV 0 and NV − ZPL transitions, with the ratio of the ZPL intensities (Fig. 3a) determined by the local Fermi-level 22, 23 and excitation intensity 24 . We also observe a broadening in the inhomogeneous NV − ZPL distribution (Fig. 3b) for device-coupled NVs vs nearby nondevice NVs. This static inhomogeneous broadening may be due to local variation in the strain environment 25,26 . Identical photon emission from different coupled defects are essential for photon-mediated defect qubit entanglement schemes 27,28 .
The observed static inhomogeneity can be bridged by Stark tuning 12,20 or quantum frequency conversion techniques 29,30 . (More information on characterization of the local strain environment for individual NV centers is provided in Supplemen-tal Information SI.4) High resolution photoluminescence excitation (PLE) spectroscopy gives us further insight into the temporal spectral stability of individual defects. In PLE measurements, a narrowband tunable laser is scanned across the NV − ZPL while collecting the NV − phonon-sideband PL (650 nm to 800 nm). From the PLE spectra we obtain the average NV − optical linewidth as well as the scan-to-scan variation in the ZPL frequency, indicating NV spectral diffusion. During PLE we can observe a loss of the NV − PL signal due to ionization to the NV 0 state. When PL is lost, we apply a short 532 nm repump pulse (0.1 s) between scans to reinitialize into the NV − charge state. Hence, the interval between repump pulses is another indicator of the NV − charge state stability.
For a deep, single NV center incorporated during growth, NV − is the preferred charge state at low excitation power 5/8 (NV − /NV 0 =3.4 at 60 µW of 532 nm excitation). An average linewidth of 44 MHz is observed (Fig. 3d.i). Minimal spectral diffusion is observed between repump pulses, with the NV − ZPL frequency exhibiting a standard deviation of 48 MHz, It is a challenge to demonstrate this level of optical stability with device-coupled implanted near-surface NVs. Pre-fabrication, the implanted single NV centers exhibit an average optical linewidth of < 100 MHz (Fig. 3c, green). Between repump pulses, the standard deviation of the NV − ZPL frequency is ∼100 MHz, indicating low spectral diffusion. (Fig. 3d.ii). The NV − charge state remains stable, no ionization is observed for multiple scans over ten minutes of measurement time. This data compares favourably with the grown-in NVs.
Post-fabrication, for both device and non-device implanted NV centers, the preferred charge state at low excitation power is NV 0 . In addition to the low NV − /NV 0 ratio (Fig. 3a dashed  lines), we observe broadening of the average single NV − linewidth (Fig. 3c, blue), rapid ionization and large spectral diffusion. We suspect that although our design avoids etching into the diamond, the GaP photonics plasma etch (Ar/Cl 2 /N 2 ) modifies the surface termination 31, 32 of the diamond and introduces new surface charge traps. Encouraged by prior studies 22, 33 , we performed a post-fabrication oxygen anneal of the sample at 400 • C (Supplemental Information SI.3). From Fig. 3a (solid lines), we see that the surface treatment more than doubles the NV − /NV 0 ratio. For the non-device NV centers, we recover the pre-fabrication average NV − linewidth (Fig. 3c, red), but the NVs retain the larger spectral diffusion (Fig. 3d.iii). Post-oxygen annealing, for device-coupled NV centers (Fig. 3d.iV), we measure an average device-coupled NV linewidth of 844 MHz (1.5 GHz pre-O 2 annealing). This linewidth is ≈ 8 times larger than non-device NVs. Further studies are needed to determine if this fast spectral diffusion can be attributed to the GaP-diamond interface (charge traps) or laser-induced due to modification of the excitation intensity profile. One promising avenue for further reduction in fabrication-induced NV instability is to perform the GaP device transfer post-fabrication via a stamping process 34 . Our inverse-design optimization framework can accommodate additional constraints such as interconnected support structures required by the stamp-transfer fabrication technique.

Discussion and conclusion
Compared with related devices in the literature such as solid immersion lenses 6,35 , diamond nanowires 36 , bullseye gratings 37 , diamond metalenses 38 and parabolic reflectors 39 , our approach enables compact devices and avoids directly etching the diamond via harsh processes that degrade the optical properties of the coupled defect qubits. Other designs based on hybrid metal-dielectric gratings / plasmonic resonances may achieve Purcell enhancement along with high directivity, but they require the NV to be either very close (few nm) to the surface 40 or embedded in the device within a nanodiamond 41-43 , both extremely challenging environments for the realization of high charge stability and spectral purity required for quan-tum information applications. Our design addresses efficient photon extraction given the constraint of an NV, formed by implantation, embedded in a low-perturbation environment within a bulk diamond sample. As near-surface defect qubit engineering advances, our inverse-design platform can readily be used for even higher enhancement of photon emission from shallower defects and provides the design flexibility for integration with emerging technologies such as stamp-transferred devices 34 .
During the preparation of this manuscript we were made aware of a contemporaneous theoretical proposal by Wambold et al. 44 . We see the two schemes as complementary, both showcasing the ability of nanophotonic inverse design for handling non-trivial design objectives and constraints. Our design is geared towards quantum information applications, while Wambold et al. is focused on quantum metrology.

Fabrication details
In Sample B, we use electronic grade CVD diamond samples (ElementSix). Based on Refs. [45,46] We perform preimplantation etching to remove the highly strained diamond surface layer (∼5µm deep) by two step reactive-ion etching (RIE) [Ar/Cl 2 , O 2 ]. Near-surface NV centers were created by 15N+ ion implantation (beam energy=85 keV, beam fluence= 3×10 10 /cm 2 , CuttingEdge Ions), followed by a twostep anneal. A 2-hour, 1200 • C annealing step was performed under vacuum (1e-7 mbar). A subsequent 4-hour 460 • C anneal was performed under constant oxygen flow in order to oxygen-terminate the surface and improve stability of the negatively charged state of the near-surface NV centers. For both samples A and B, a 250 nm thick GaP membrane was transferred to the diamond via epitaxial liftoff in hydrofluoric acid and Van der Waals bonding. We spin-on hydrogen silsesquioxane (HSQ) resist and pattern the devices by electron-beam lithography. RIE with an Oxford PlasmaLab-100 ICP etcher was then used to etch the devices (3.0 mTorr, 1.0/6.0/3.0 sccm Cl2/Ar/N2, 235 V dc bias). More fabrication details are provided in the Supplemental Information SI.2. Post-fabrication oxygen anneals were performed for 4-hours at 400 • C. The lower O 2 annealing temperature was used to prevent damaged to the fabricated GaP structures. We observed that annealing at 425 • C caused irreversible damage to the GaP layer destroying the photon extractors.

Measurement details
Measurements on sample A (B) were performed at 298K (8K) with a custom built confocal microscope. For low temperature measurements the samples are placed in a closed-cycle He cryostat (Janis SHI-4XG-5-M). The off-resonance NV measurements are performed by optically exciting with a 532 nm diode-pumped solid-state laser (RGBLase LLC FB-532) at powers between 0.1-5 mW with a spot diameter of ∼800 nm. For the photon autocorrelation measurements, the NV − ZPL photons are filtered through a spectrometer 6/8 (Princeton Acton 2750, 1200g grating). NV ZPL photons are measured with avalanche photodiodes (free-space: Excelitas SPCM-AQRH-16). A IDQuantiq 801 TDC is used to record photon arrival times. For NV − resonant excitation measurements, a Newfocus velocity tunable 637nm laser is utilized. When required 2.88 GHz sidebands can be added to the velocity laser using a fiber-coupled EOM (EOspace PM-0K1-PFA-637). The filtered NV − PSB photons are measured with fiber-coupled photodiodes (Excelitas SPCM-AQ4C). The reported linewidths in Fig. 3(c, d) are calculated by fitting each scan to a Lorentzian and then averaging all the fitted linewidths. More measurement details are provided in the Supplemental Information SI.1.

2
The NV centers are optically excited with a 532 nm diode-pumped solid-state laser (RGBLase LLC FB-532) at powers between 0.1-5 mW with a spot diameter of ∼800 nm. A polarizing beamsplitter with an automated half-wave plate is used in the excitation path to preferentially excite a given NV orientation. For off-resonant measurements, a combination of dichroic (Semrock 532LP) and 532 nm razor edge filter is used to attenuate the reflected excitation signal. Spectra are taken with a Princeton Acton 2750 spectrometer with either a 300, 1200 or 1800 groove grating. For resonant excitation studies a New Focus Velocity tunable 637 nm external-cavity diode laser reflected from a 90/10 T/R beam sampler is used in the excitation path. The phonon-sideband (PSB) emission is filtered in the collection path by a combination of 532 nm razor edge and a 660-800 nm bandpass filter (Semrock FF01-731/137-25). The NV − PSB photons are measured with a fiber coupled avalanche photodiode (Excelitas SPCM-AQ4C). Photoluminesence from Sample A is attenuated with an optical density filter (Newport-530-OD1) to avoid detector saturation.
On sample B, post-fabrication, background fluorescence is observed under 532 nm excitation from the resist (HSQ) and GaP structure. To mitigate this issue only the ZPL is utilized for off-resonant optical characterization. For the photon autocorrelation measurements, the NV − ZPL photons are filtered through the spectrometer with the 1200 groove grating. The filtered photons are then incident on a 50/50 plate beam splitter, with the two output modes measured by free-space coupled low dark count (<30 cts/s) avalanche photodiodes (Excelitas SPCM-AQRH-16). Two additional filters (Semrock 750SP, 631/25BP) are introduced to reduce spurious counts from detector afterpulsing. The photons counts are recorded on two independent channels by a timing board (ID Quantique 801 TDC).
During resonant excitation measurements, if the NV center under study is in a high strain environment (evident by the splitting of the E x and E y excited state transitions under off-resonant excitation), sidebands at 2.88 GHz are added to the resonant scanning laser using a fiber-coupled electro-optic modulator (EOSPACE PM-0K1-PFA-637). Electron beam lithography (100 kV, 2 nA, JEOL-6300) is utilized to pattern the designed photon extractor to a thin (150 nm) HSQ resist layer on top of the GaP-on-diamond stack. However, charging can have a severe impact on the patterning resolution because of the non-conductive diamond substrate even with the thin GaP membrane. To mitigate charging, we spin a water soluble polymer (4.5% PSS + 1% TritonX) over the HSQ resist layer and sputter ∼ 5 nm of Au+Pd metal. This metal layer readily lifts-off during development of the HSQ resist. Multiple arrays of devices are patterned with the design feature size and (x,y) aspect ratio are varying in 12.5 nm increments across the array to compensate for fabrication variations.The fabricated devices were observed under scanning electron microscopy for process validation and determination of device yield. It can be seen in Fig. SI.2 that on average each feature is within 10% of the designed dimensions.

SI.3. OXYGEN ANNEALING DEVICES POST FABRICATION
Post fabrication and after measurements presented in Fig.2, both samples A and B were annealed at 400 • C for 4 hours under constant oxygen flow. On sample A, we observe an increase in overall NV PL emission (Fig. SI.3a). The ensemble coupled devices show a similar trend. After oxygen annealing sample B, we do not observe any notable changes to PL intensity. The observed ZPL intensity for a device-coupled single NV (Fig. SI.3b) is similar to pre-anneal. As discussed in the main text, we see a significant improvement in NV − /NV 0 ratio, increased charge state stability and reduced spectral diffusion. Further annealing of sample A at 425 • C caused irreversible damage to the GaP layer destroying the photon extractors.  The NV − excited state structure is sensitive to the local strain or electric field perturbations [1][2][3]. Once the strain perturbation exceeds the spin-orbit splittings, additional transverse strain increases the energy splitting between the Ex and Ey orbital branches whereas a perturbation parallel to the NV axis shifts the entire excited state manifold. We can observe this splitting under off-resonant excitation (532 nm) (Fig. SI.4a) in low temperature (T = 8 K) spectra. Both device-coupled and non-device NV centers exhibit varying magnitudes of splitting, from < 8.2 GHz (resolution limited by spectrometer) to ∼90 GHz. The two transitions (E x , E y ) exhibit near-orthogonal excitation polarization dependence (Fig. SI.4a, inset). We confirm that the observed (E x , E y ) pairs correspond to single NV centers by performing autocorrelation measurements (Fig. SI.4b).
The observed implanted NV (E x , E y ) splitting does not preclude their use for identical singlephoton generation as it is within the capabilities demonstrated by Stark tuning [2,4]. However, the high strain regime can be detrimental to quantum information schemes that rely on high fidelity single-shot qubit spin readout [5]. In a high transverse strain environment, the probability of a spin-flip during an optical transition increases, placing the NV − center in an effective dark state. This is illustrated in our resonant excitation measurements. To perform resonant excitation studies, we have to effectively mix the ground state m s = 0 and m s = ±1 populations. We achieve this by adding 2.88 GHz sidebands to our resonant scanning laser to simultaneously drive both spin states (Sec. S.1).   (2) (τ ) data. b. g (2) (τ ) for 120 µW non-resonant excitation. c. g (2) (τ ) for 300 µW off-resonant excitation. d. g (2) (τ ) for 600 µW non-resonant excitation. Note the increase in bunching peaks with increase in excitation intensity.

SI.5. MODEL FOR AUTOCORRELATION (g (2) (t)) MEASUREMENTS
The normalized autocorrelation measurements on the NV − ZPL photons emitted shown in Fig. SI.5(b-d) provide definitive evidence that we are observing enhancement from a single quantum emitter, yet also reveal significant photon bunching. Power dependent photon bunching is a well-documented effect in NV centers, however the number of coincidences we observed at delays around 20 ns for the 300 µW and 600 µW data is larger than typically observed in literature. We observe similar bunching in implanted, non-device coupled NV centers in the same sample, but under higher excitation intensity. This similarity suggests that the devices enhance the excitation field at the defect, but are not responsible for the presence of the bunching. Photon bunching is detrimental for many applications since it arises from non-cycling transitions shelving the NV center into inaccessible, or dark, states.
The connection between the photon statistics and the NV center's state is given by the following where p(e, τ |g, 0) is the probability of finding the NV center in the excited state at time delay τ given that it was in the ground state at 0 time delay, and p(e, ∞) is the steady state population in the excited state under continuous excitation.
We model the NV center as a 6 level system in order to account for the NV − ground states, excited states, singlet states, and NV 0 states as shown in Fig. SI.5(a). In our labelling scheme, manifolds 1 and 3 corresponds to the m s = 0 ground and excited states respectively. Manifolds 2 and 4 correspond to the m s = ±1 ground and excited states. Manifold 5 contains the singlet states, and manifold 6 consists of the NV 0 states. The probability of being in level i is p i and the rate from level i to j is given by k ij , and P is the optical excitation power. Coherences between the different energy levels are assumed to decay on time scales faster than the excitation/relaxation rates due to the non-resonant pumping, which simplifies the model to following set of 6 rate equations.
2) One subtlety is that each level of our model corresponds with a manifold of NV states. This has implications in our treatment of off-resonant excitation into the excited state, the multiplicity of the m s = ±1 level, and the NV 0 internal dynamics. We address the issue of off-resonant excitation by taking the excitation rate to depend linearly on power, while the photon emission rate is power independent. We account for the multiplicity of the m s = ±1 states by properly normalizing the branching ratios. For example, the recombination from the NV 0 manifold into m s = ±1 is taken to be twice the recombination into m s = 0. Lastly, we ignore the internal dynamics of NV 0 , and treat recombination as a two-photon process that in the weak excitation limit should scale as the power squared.
The initial condition is determined by the spin polarization, p 1 /(p 1 + p 2 ), and is determined by forcing our solution to be self consistent. In other words p 1 (0)/(p 1 (0) + p 2 (0)) = p 1 (∞)/(p 1 (∞) + p 2 (∞)). We find that the equilibrium spin polarization approaches a maximum value of 0.8 at low power and decreases with power.
Older models of NV dynamics ignore ionization and recombination, and focus on the role of the NV − singlet state and spin non-conserving transitions [6]. We are unable to reproduce the power dependent bunching we observe with the rates from these models. Our model is more similar to recent literature models which account for both the long-lived NV − singlet state, as well as the NV 0 charge state, but we ignore internal dynamics within the singlet and neutral states [7,8] as we are able to obtain satisfactory agreement simulataneously fitting the three power-dependent g (2) curves with the simplified model (Fig. SI.5. We begin the optimization with rates taken directly from these models and optimize over each parameter. This is necessary, since the rates can vary significantly between different NV centers (Table SI.1). We find the final optimized rates are in reasonable agreement with reported values, with the only disagreement being in lower shelving rates into the singlet state, and lower deshelving rates out of the singlet states. Nevertheless, it is clear that ionization and recombination are important factors in the emission properties of an NV center, and should be considered in a variety of applications.