Probing the Optical Dynamics of Quantum Emitters in Hexagonal Boron Nitride

Hexagonal boron nitride is a van der Waals material that hosts visible-wavelength quantum emitters at room temperature. However, experimental identification of the quantum emitters' electronic structure is lacking, and key details of their charge and spin properties remain unknown. Here, we probe the optical dynamics of quantum emitters in hexagonal boron nitride using photon emission correlation spectroscopy. Several quantum emitters exhibit ideal single-photon emission with noise-limited photon antibunching, $g^{(2)}(0)=0$. The photoluminescence emission lineshapes are consistent with individual vibronic transitions. However, polarization-resolved excitation and emission suggests the role of multiple optical transitions, and photon emission correlation spectroscopy reveals complicated optical dynamics associated with excitation and relaxation through multiple electronic excited states. We compare the experimental results to quantitative optical dynamics simulations, develop electronic structure models that are consistent with the observations, and discuss the results in the context of ab initio theoretical calculations.

Even less is known about the visible QEs' optical dynamics.Optical dynamics arise from a QE's electronic structure together with radiative and nonradiative transitions between electronic states.State transitions can involve multiple processes including electron-phonon interactions, intersystem crossings between different spin manifolds, and ionization or recombination events.For QEs in h-BN, previous studies have reported photon bunching associated with metastable dark states [3,6,11], and yet the nature of these states and the transitions between them remains unclear.Some QEs exhibit magnetic-field-dependent modulation of their photoluminescence (PL) signal, consistent with a spin-dependent intersystem crossing, whereas others do not [11,14].An optically detected magnetic resonance signal was observed for a particular QE under excitation at 633 nm but not at 532 nm [13].Such observations present a complicated picture of the visible QEs, likely involving multiple defect classes (e.g., different chemical structures or charge states), strong local perturbations, and complex excitation and relaxation pathways.Improved understanding of the QEs' optical dynamics can resolve these mysteries.Such understanding is also a prerequisite to designing quantum control protocols that would facilitate their use in quantum technologies.Here, we use quantitative spectral, spatial, and temporal PL spectroscopy to investigate the optical dynamics of h-BN's QEs.
Photons emitted by a QE carry a wealth of information about its electronic structure and optical dynamics.For vibronic optical transitions, the photon energy and polarization distributions reflect the details of electron-phonon coupling and optical dipole selection rules, respectively.The QEs in h-BN generally exhibit linearly-polarized PL and strong electron-phonon coupling associated with a single vibronic transition [6], and yet other experimental and theoretical evi-dence points to the involvement of multiple excited states in the optical dynamics [10,28,43,46].Time-dependent measurements provide complementary information.The secondorder photon autocorrelation function, g (2) (τ ), is widely used to identify single-photon emitters.As a more general analytical tool, photon emission correlation spectroscopy (PECS) yields quantitative information about a QE's optical dynamics [47].Qualitatively, we distinguish between photon antibunching (g (2) (τ ) < 1) as a signature of non-classical light, with single-photon emission as a special case when g (2) (0) = 0, and photon bunching (g (2) (τ ) > 1 for τ = 0) as a signature of dark, metastable states accessed via nonradiative transitions.Quantitative measurements of g (2) (τ ) as a function of optical excitation power or wavelength can elucidate a QE's excitation and emission pathways as well as bunching mechanisms.
Prior observations of h-BN's visible QEs feature both bunching and antibunching signatures, although with some unusual, conflicting patterns.Some QEs respond to applied dc and ac magnetic fields in a manner consistent with spinmediated intersystem crossing transitions, whereas others do not [11,13,14].A lack of evidence for pure single-photon emission motivated a proposal that h-BN's QEs occur in pairs as "double defects" [48].In this work, we compare quantitative PL spectroscopy and PECS measurements of h-BN's QEs with theoretical simulations.We show that QEs in roomtemperature h-BN can exhibit pure single-photon emission, with g (2) (0) = 0 within experimental uncertainty.Furthermore, we find evidence for multiple electronic states connected by radiative and nonradiative transitions, with associated timescales spanning over five orders of magnitude.Comparing the experiments to theoretical proposals, we find that the boron dangling bond model provides a consistent, quantitative understanding of the observations for individual QEs as well as their heterogeneity.

II. RESULTS
The Results and Discussion are organized as follows.First, we report the basic optical characteristics of five well-isolated QEs across three samples, specifically including their PL spectra, PL saturation as a function of excitation power, polarization properties, and single-photon purity.Next, we investigate the QEs' optical dynamics using PECS as a function of excitation power and wavelength.We consider different models for the electronic level structure and simulate the corresponding optical dynamics.Finally, we compare the theoretical simulations to experimental observations and discuss the implications for understanding the QEs' electronic and chemical structure along with their spin and charge dynamics.

A. Photoluminescence characterization
We used a custom-built confocal microscope to study individual QEs in h-BN under ambient conditions.The h-BN bulk crystals were mechanically exfoliated into thin flakes us-ing a dry transfer process [49] and transferred to a SiO 2 /Si substrate patterned with circular trenches [6].Prior to the optical studies, the samples were annealed in a tube furnace at 850 °C in Ar atmosphere for 2 hours.This annealing process has been shown to brighten the emitters [35].The QEs are illuminated with either of two continuous-wave lasers operating at 532 nm and 592 nm wavelengths, where excitation power and polarization are controlled.To differentiate between excitation wavelengths in this work, data recorded under 532 nm (592 nm) excitation are plotted in green (orange) in the relevant figures.Some QEs disappeared during experiments, hence the set of measurements is not identical for each QE.See Materials and Methods for additional details on sample preparation, data acquisition, and analysis.
Figure 1 summarizes the PL characterization measurements.Each row corresponds to a particular QE (labeled A-E), and each column corresponds to a different experiment.The first column includes µ-PL images of each QE, acquired by scanning a fast steering mirror and recording the accumulated counts at each pixel.A two-dimensional Gaussian fit to each µ-PL image yields the background and signal levels for subsequent studies.The second column displays g (2) (τ ) measurements over short delay times, showing characteristic antibunching dips fit by an empirical model for a multi-level system (see Materials and Methods).The third column displays the steady-state PL signal as a function of excitation power.These data are fit using an empirical saturation model, where C is the background-subtracted, steady-state PL count rate, P is the optical excitation power, C λ sat is the saturation count rate for a particular excitation wavelength, λ, and P λ sat is the corresponding saturation power.The best-fit results are reported in Supplementary Table S2.
The fourth column of Fig. 1 presents PL emission spectra and polarization measurements.In each PL spectrum, the long pass filter cut-on wavelength is indicated as a vertical dotted line, and the excitation wavelength is a solid line.The inset to each PL spectra panel presents measurements of the QE's excitation and emission polarization properties.These data are acquired by varying the linear polarization of the excitation laser (colored circles) or by passing the PL through a linear polarizer placed in the collection path (black squares).For excitation polarization measurements, the linear polarizer in the collection path is removed.For emission polarization measurement, the excitation polarization is set to maximize the PL.At each polarization setting, we record the steady-state PL intensity as well as a background intensity from a spatial location offset ∼1 µm from the QE, which is subtracted to yield the PL signal.The order of the polarization angles is set randomly to minimize effects of drift and hysteresis.Solid curves are fits to the data using the model function where λ indicates the excitation wavelength, s indicates excitation (ex) or emission (em), A λ s is the amplitude, θ λ s is the polarization angle of maximum intensity, and B λ s is the offset.From the fit results, the visibility is calculated as where I max s and I min s are the maximum and minimum intensities, respectively.The misalignment between the excitation and emission polarization angles is The best-fit parameters are reported in Supplementary Table S2.

B. Spectral emission lineshapes
For spectra that are not cut off by the excitation filter (namely, the 532 nm excitation spectra for QEs A, C, and E, and the 592 nm spectrum for QE B), we find that the lineshapes are consistent with a Huang-Rhys model for a vibronic transition associated with a single zero-phonon line (ZPL).We use the analysis method described by Ref. [6] to fit the observed PL spectra using an empirical model in which the ZPL energy, ZPL width, Huang-Rhys factor, and vibronic coupling lineshape are free parameters; see Materials and Methods for additional details.The results are shown in Fig. 2. In the left column, we plot the normalized observed emission lineshape, L(∆E) ∝ S(∆E)/E 3 , where S(∆E) is the spectral intensity distribution as a function of the relative energy ∆E = E ZPL − E, with E denoting the photon emission energy and E ZPL denoting the ZPL energy.The factor 1/E 3 accounts for the photon energy dependence in spontaneous emission.Each solid curve is the result of a weighted leastsquares fit of the model to the experimental lineshapes.The right column of Fig. 2 shows the corresponding 1-phonon vibronic coupling lineshape for each fit.Best-fit parameters are reported in Supplementary Table S2.

C. Photon emission correlation spectroscopy
Temporal correlations between fluorescence photons reveal information about a QE's excitation and emission dynamics.In this work, we use PECS for two purposes: to verify the single-photon purity of the QEs and to probe their optical dynamics as a function of optical excitation rate.We calculate g (2) (τ ) from the photon arrival times acquired from two detectors in a Hanbury Brown and Twiss interferometer using a time-correlated single-photon counting module.For QEs in h-BN, the timescales over which antibunching and bunching occur can vary over at least 6 orders of magnitude [6,11,16].For this reason, we initially calculate and analyze g (2) (τ ) over a logarithmic scale spanning from 100 ps to 1 s.We fit the background-corrected data using a general empirical model for a QE's optical dynamics with a varying number of levels: Here, γ 1 is the antibunching rate, C 1 is the antibunching amplitude, γ i for i ≥ 2 are bunching rates, and C i for i ≥ 2 are the corresponding bunching amplitudes.We determine the number of resolvable timescales, n, by calculating and comparing the Akaike Information Criterion (AIC) and the reduced chi-squared statistic for each best-fit model.In optical dynamics models, an N -level system is characterized by N − 1 rates, corresponding to the nonzero eigenvalues of the generator matrix (see, e.g., Eq. ( 12) in Materials and Methods).Therefore, the inferred value of n places a lower limit on the number of electronic levels required to describe the observations, N ≥ n + 1.We extract the rates, amplitudes, and their corresponding uncertainty from these fits for comparisons with theoretical simulations.In order to assess the single-photon purity associated with the value of g (2) (0), we perform a subsequent analysis of g (2) (τ ) calculated over a linear scale of delay times, τ ∈ [− 20,20] ns.Examples of such data are shown in Fig. 1, along with constrained fits in which only the antibunching parameters γ 1 and C 1 are allowed to vary, and which account for the instrument response function associated with detector timing jitter.See Materials and Methods for further details.

D. Verifying single-photon emission
Any observation of sub-Poissonian statistics, g (2) (0) < 1, indicates the presence of quantized photon emission.The threshold g (2) (0) < 0.5 is often used to indicate singlephoton emission; however, a more precise interpretation is that a PL signal is dominated by a single-photon emitter when (0) FIG. 3. Single-photon emission characteristics.Zero-delay photon autocorrelation function as a function of optical excitation power.Orange (green) data correspond to excitation at 592 nm (532 nm).
Error bars represent one standard deviation derived from fits as described in the text.
g (2) (0) < 0.5 [47].An observation of g (2) (0) > 0 implies a non-zero probability of observing two detection events simultaneously, either due to background fluorescence, detection timing jitter, or the presence of multiple QEs.Studies of h-BN's QEs routinely report g (2) (0) < 0.5, however we are unaware of any prior room-temperature observations of pure single-photon emission with g (2) (0) = 0. Partially on the basis of such observations, Ref. [48] proposed that h-BN's QEs actually occur in pairs as double defects with parallel emission pathways.
We find that QEs in h-BN can indeed exhibit pure singlephoton emission at room temperature.Figure 3 shows g (2) (0) for each QE as a function of excitation power.These data are corrected for background fluorescence and detector timing jitter, as described in Materials and Methods.For QEs C, D and E, we observe g (2) (0) = 0 within the experimental uncertainty, particularly at low excitation powers.For QEs A and B, we observe g (2) (0) ∼ 0.1−0.2.The offset from zero could reflect a contribution from additional dim emitters, however we believe it is more likely to result from incomplete estimation of the background.For instance, QE B sits on an extended background feature whose contribution is not captured by our standard analysis method.For QE C, we attribute the increase in g (2) (0) as a function of excitation power to limitations in the instrument-response-function correction as the antibunching rate exceeds the detector timing resolution.

E. Probing the optical dynamics
Figure 4 summarizes the results of fitting the empirical model of Eq. ( 5) to PECS measurements as a function of optical excitation power.The figure includes the best-fit antibunching rate (top row) as well as the first two bunching rates and amplitudes (lower rows).The PECS data for QEs B, D, and E are best described by a three-timescale model (n = 3), whereas QE A exhibits four resolvable timescales (n = 4).For QE C, we resolve two or three timescales depending on the excitation power and wavelength.The best-fit results for QE A's third bunching component (γ 4 and C 4 ), as well as the antibunching amplitude (C 1 ) for all emitters are shown in Supplementary Figures S6 and S7, respectively.As in previous figures, colors indicate the optical excitation wavelength.The PL decay rate of QE A was directly measured to be 355 MHz using a picosecond pulsed laser (see Section S6 and Fig. S4 in the Supplementary Materials); this measurement is shown in Fig. 4(a) as a dashed black line.The PL lifetime measurement was only performed for QE A given the susceptibility of h-BN's QEs to disappear under pulsed excitation.
To fit these metadata, we consider the following empirical

models:
Model I (Linear) : Model III (Second-Order Saturation) : Model IV (Quadratic) : where P is the excitation power and the other variables are free parameters representing zero-power offset, R 0 , lowpower slope, m 0 , high power slope, m 1 , saturation value, R sat , and saturation power, P sat .Dotted curves in Fig. 4 show the best-fit results for the model listed in each panel; in each case, we select the model with the fewest free parameters that qualitatively fits the data.Best-fit parameters and uncertain-ties for each fit are reported in Supplementary Tables S5 and  S6.
The antibunching rate, γ 1 , exhibits a markedly nonlinear power dependence for QEs A, B, and C whereas the dependence appears to be linear for QEs D and E. However, we note that the power range in the data for QEs D and E might not be large enough for nonlinearities to emerge.For comparison, QEs B and C are excited with up to ∼ 4P λ sat whereas QE D is excited with up to ∼ 2P λ sat (see Supplementary Table S2).The zero-power antibunching-rate offset (R 0 ) for QEs B-E is clearly nonzero, whereas the fits using Model II for QE A are poorly constrained, yielding R 0 = 0 ± 261 MHz and R 0 = 0 ± 167 MHz for green and orange excitation, respectively.The antibunching amplitudes (see Supplementary Figure S7) for all QEs show a nonlinear saturation dependence on excitation power with an expected convergence to C 1 ∼ 1 at zero excitation power.
The bunching dynamics exhibit significant quantitative and qualitative variations across emitters.The fastest bunching rate, γ 2 , scales linearly with excitation power and has a nonzero offset for QEs A, D and E, whereas it exhibits saturation behavior and zero offset for QEs B and C. The magnitudes of γ 2 range from several kilohertz (QEs A, B, and D) up to several megahertz or faster (QEs C and E).The slower bunching rate, γ 3 , exhibits the largest qualitative variation across emitters, including linear (QE D), quadratic (QEs A and E), and saturation models (QEs B and C).Only QE D exhibits clear evidence for a non-zero offset for γ 3 .The magnitudes of γ 3 are typically in the kilohertz range, with the exception of QE C, whose γ 3 increases beyond 1 MHz at high powers.The bunching amplitudes primarily depend nonlinearly on excitation power, except for C 2 of QEs B, C (green excitation) and E, which scale linearly with excitation power.All of the bunching-amplitude fits are consistent with zero offset, except for QE E, where small residual offsets (R 0 < 0.1) likely reflect minor systematic errors in the analysis or inaccuracies of the empirical models.For QE A, we restrict the meta-analysis of bunching parameters to the orange-excitation data, which extend to higher excitation power.However, we note that the green-excitation bunching parameters generally track the data for orange excitation.

F. Electronic model and optical dynamics simulations
We find that the key features observed in Fig. 4 can be understood using the four-level electronic model shown in Fig. 5(a).Figure 5 summarizes the results of optical dynamics simulations for this model.Given a set of transition rates for the model, we simulate g (2) (τ ) including the effects of timing resolution and shot noise (e.g., Fig. 5(b)), and we subsequently fit the simulated data using the empirical model of Eq. ( 5) with n = 2 to extract the antibunching and bunching parameters, as shown in Figs.5(c)-(e).For reasons explained later in this section, the simulated data were best described by an n = 2 model despite having three eigenvalues.See Materials and Methods for more information regarding the simulations.
The four-level model consists of a ground state (level 1), an excited radiative state (level 2), a higher-lying excited state (level 3) and a nonradiative metastable state (level 4).We consider two optical excitation pathways from the ground state to excited states 2 or 3, represented by the rates Γ 12 and Γ 13 , respectively.The magnitudes of these two rates depend on the corresponding optical cross-sections for absorption at the excitation wavelength.A difference in cross section can result from the difference in electric dipole matrix elements between the different electronic states, the atomic configuration coordinate overlap for vibronic transitions, or both of these factors.For the simulations in Fig. 5, we set Γ 12 = 0, since we are particularly interested in the situation where Γ 12 /Γ 13 1, such that the dynamics feature indirect excitation of the radiative state 2 via nonradiative relaxation from excited state 3, at a rate κ 32 .This was informed by the nonlinear powerscaling of γ 1 for QEs A, B and C. In the Supplementary Materials, we report simulations over a range of settings where Γ 12 /Γ 13 ∈ [0, 2], with qualitatively similar results (see Supplementary Figure S8).
In addition to the indirect excitation pathway formed by states 1, 2, and 3, optical excitation results in population and relaxation of metastable state 4 via nonradiative transitions with rates κ 24 and κ 41 .We consider two types of nonradiative transition mechanism for the metastable state: spontaneous and optically pumped.Spontaneous transition rates are independent of the optical excitation rate (in this case, Γ 13 ), whereas optically pumped transition rates scale linearly with Γ 13 .In this model, the optically pumped transition rates κ 24 and κ 41 can approximate more complicated processes; for example, they could involve re-pumping from levels 2 → 3 or from levels 4 → 3 with subsequent nonradiative relaxation (see Supplementary Figure S10), or they could involve transient population of additional levels.Their approximation as individual pumped transitions remains accurate as long as optical pumping remains the rate-limiting step.The key observable difference between spontaneous and optically pumped transitions manifests in the excitation power dependence of the corresponding bunching rate (Fig. 5d); the bunching rate for spontaneous transitions features a non-zero zero-power offset and saturates at high power, whereas the bunching rate for optically pumped transitions has zero offset and scales nearly linearly with power, even as the corresponding bunching amplitude (Fig. 5e) saturates.
For both bunching mechanisms, the simulated data were best described by only a single bunching level (n = 2 in Eq. 5) despite the fact that there should be 3 eigenvalues that describe this system.The reason for this is that the indirect excitation and emission process through levels 1, 2, and 3 can lead to two of the eigenvalues being complex.These eigenvalues have the largest real values and are responsible for the antibunching dynamics.When we include practical limitations on timing resolution and signal-to-noise ratio at short delay times, the fit cannot distinguish these two values, and the goodness-of-fit analysis prefers a single real rate that approximates the true model.The result is an effective antibunching rate that scales nonlinearly with increasing excitation rate.This effect persists even when a direct transition from state 1 → 2 is included.We performed simulations varying the ratio Γ 12 /Γ 13 , and observed qualitatively similar results (see Supplementary Figure S8).

A. Photoluminescence, spectral, and polarization properties
Our experiments provide clear evidence that visible QEs in h-BN occur as isolated point defects with emission originating from a single, dominant optical transition.The QEs are spatially resolved in high signal-to-background µ-PL images.They exhibit PL saturation, high polarization visibility in emission, and emission spectra consistent with individual vibronic transitions.Most convincingly, several emitters exhibit pure single-photon emission, with g (2) (0) = 0 within small experimental uncertainty, as shown in Fig. 3.This finding contrasts with previous suggestions that h-BN's QEs occur in pairs [48].We do not contend, however, that such pairing cannot occur.On the contrary, in the course of our exper-FIG.5. Electronic level structure and optical dynamics simulations.(a) An investigative four-level electronic model consisting of the ground state (level 1), radiative state (level 2), excited state (level 3), and metastable state (level 4).Orange arrows represent excitation pathways (with rates Γ12 and Γ13), the wavy red arrow represents radiative emission (with rate Γ21), and dotted black arrows represent nonradiative transitions (with rates κ32, κ24 and κ41).(b) Simulated g (2) (τ ) for Γ12 = 0, Γ13 = 84 MHz, Γ21 = 300 MHz, κ32 = 600 MHz, κ24 = 60 kHz, and κ41 = 30 kHz.Error bars represent simulated photon shot noise.(c-e) Best-fit parameters γ1, γ2, and C2 determined by fitting simulated g (2) (τ ) data using Eq. 5 with n = 2.The results are plotted as a function of Γ13/Γ21, where Γ21 = 300 MHz is a fixed parameter.
iments we observed multiple instances of spatially isolated emitters with high polarization visibility, and yet g (2) (0) is substantially larger than zero.We have focused here on emitters showing the highest likelihood of being single defects.
Qualitatively, the QEs' room-temperature PL spectra are similar to those reported elsewhere in the literature [6,19,27,28].The analysis shown in Fig. 2 indicates that the PL spectra are consistent with individual vibronic transitions between two optical-dipole-coupled electronic states.The one-phonon lineshapes for QEs A, B, C, and E are all qualitatively similar despite the fact that QEs A and B feature ZPL energies near 1.9 eV, compared to 2.1 eV for QEs C and E. All four emitters exhibit strong coupling to low-energy phonons ( 50 meV) as well as to higher-energy phonons (150-200 meV) that are typically associated with longitudinal optical modes in bulk h-BN [4].Coupling to low-energy phonons is a key feature in determining the asymmetric shape of the dominant emission peak [27,50].The ZPL corresponds to the transition from the lowest vibrational level of the initial (excited) state to the lowest vibrational level of the final (ground) state.When low-energy phonons are involved, transitions can occur from the lowest vibrational level of the initial state to the first vibrational level of the final state, showing up on the low-energy side of the ZPL and leading to the asymmetric shape.Failure to account for low-energy phonons in interpreting experimental spectra leads to an underestimation of the Huang-Rhys factor, S HR , which quantifies the strength of the vibronic coupling and is a crucial parameter for comparing with theoretical calculations.Our model captures the asymmetric spectral shape.However, the precise details of the low-energy phonon-coupling lineshape become correlated with the ZPL width (assumed to be Lorentzian) and S HR when fitting the model to experimental data.We account for these correlations by performing the fits using varied constraints on the low-energy phonon coupling motivated by scaling considerations.We follow the method described in Ref. [6], in order to estimate uncertainties on S HR and the ZPL linewidth.Overall, we find that the best-fit ZPL linewidths are narrower or comparable to those reported in the literature for off-resonant excitation of h-BN's QEs at room temperature [10,51,52], and the values of S HR are somewhat higher.We consider comparisons to theoretical proposals in detail later; we note here that the ZPL energies and S HR values closely match the calculated properties of boron dangling bonds [43,53].
The QE's polarization-resolved PL excitation and emission characteristics (Fig. 1) begin to reveal more complicated features of their optical dynamics.Both QEs A and C exhibit linearly polarized emission with nearly complete visibility, again consistent with emission through a single optical dipole transition.For QE A, the PL intensity varies as a function of excitation polarization angle in a manner consistent with excitation through a single optical dipole, with high visibility and an angle aligned with the emission dipole, independent of excitation wavelength (532 nm or 592 nm).In the case of QE C, the emission polarization visibility and dipole angle is similarly independent of the excitation wavelength.However, QE C's excitation polarization dependence varies dramatically as a function of excitation wavelength; the excitation dipole is aligned with the emission under 592 nm excitation, but misaligned under 532 nm excitation with substantially reduced visibility.Emission polarization data are not available for QEs B, D, and E, but the excitation polarization measurements are qualitatively similar to those for QEs A and C. All three QEs show polarized absorption with varying degrees of visibility.
The heterogeneous polarization responses are consistent with previous observations for QEs in h-BN [6,18,27,28].In particular, Ref. [28] studied the variation of polarization visibility and alignment between excitation and emission as a function of the energy difference between the excitation photon energy and the ZPL photon energy, ∆E.They observed that the excitation and emission dipoles are aligned (∆θ = 0) when ∆E 200 meV, whereas if ∆E 200 meV, ∆θ can take any value.Our observations are consistent with this empirical finding.For QE A, the excitation and emission angles are aligned despite relatively large energy differences, ∆E 592 = 169 meV and ∆E 532 = 405 meV, for 592 nm and 532 nm excitation, respectively.For QE C, the dipoles are aligned for excitation at 592 nm (∆E 592 = 10 meV; ∆θ 592 = 0.0 ± 1.4°) but misaligned at 532 nm (∆E 532 = 247 meV; ∆θ 532 = 46.5 ± 1.6°).
Misalignment between absorption and emission dipoles is expected if the optical dynamics involve multiple excited states.Whereas the invariance of the PL polarization, visibility, and spectral shape to exctiation energy implies that PL emission occurs through a single optical transition, offresonant optical pumping can involve transient excitation of higher-lying excited states through transitions with different optical dipole orientations, which subsequently relax to the radiative state as shown in Fig. 5(a).Depending on energy level arrangement and the vibronic copuling strengths, a single excitation laser can drive both transitions between states 1 → 2 and 1 → 3. The excitation polarization dependence will then reflect a superposition of two optical dipole transitions, with an orientation and visibility determined by the underlying dipole transition orientations and their relative optical cross section.To test this hypothesis, we performed a simultaneous fit of QE C's emission and excitation polarization data under 532 nm excitation assuming a single shared dipole for excitation and emission via states 1 ↔ 2 together with a second dipole for excitation via 1 → 3 (see Section S4 and Fig. S2 in the Supplementary Materials).We find that the data are consistent with such a model, in which the dipole projection for transition 1 → 3 is misaligned from that of transition 1 ↔ 2 by 63 ± 1°, and the ratio of excitation cross sections is Γ 12 /Γ 13 ∼ 0.5.
In interpreting these results, we note that the observation of highly polarized emission implies the presence of at least one symmetry axis for the underlying electronic states.For most defect models under consideration, symmetry allows for optical dipole transitions aligned perpendicular to the h-BN plane (along z) or within the plane either parallel or perpendicular to the defect's symmetry axis (along x or y).Hence, the observation of dipoles misaligned by ∼60°seems surprising.However, since our polarization-resolved experiments are primarily sensitive to the projection of the dipole perpendicular to the microscope's optical axis, it is possible that sample misalignment or local distortions of the defect that tend to tilt the z axis could explain the observations.Alternatively, our model of two superposed excitation dipoles might not capture all salient features of the excitation process; more than two transitions might be involved, and yet the superposition of any number of dipole absorption patterns will ultimately yield a polarization dependence consistent with Eq. (2).

B. Optical dynamics
PECS experiments reveal key details regarding the nature of the QEs' excited states and optical dynamics.The PECS results summarized in Fig. 4 resolve individual dynamical processes, their associated timescales, and their dependence on optical excitation power.All QEs feature three or more timescales in their autocorrelation spectra, which implies that the optical dynamics involve at least four electronic levels.In addition to antibunching on nanosecond timescales, all QEs exhibit bunching with two or more resolvable timescales that are orders of magnitude longer (typically microseconds to milliseconds).These bunching timescales are broadly consistent with past observations [3,4,6,14,26,54], and they indicate the role of metastable dark states in the optical dynamics.Here, we emphasize and discuss two key features of the PECS measurements in Fig. 4: the nonlinear power dependence of the antibunching rate, γ 1 , that is clearly observed for QEs A, B, and C; and the heterogeneous behavior of the bunching rates and amplitudes, which feature qualitatively diverse power-dependent variations.
For a QE featuring a direct optical transition between a ground state and a radiative excited state, the antibunching rate scales linearly as a function of optical excitation power, with a zero-power offset corresponding to the QE's spontaneous emission rate.This is the case even for QEs that also feature metastable charge and spin states, such as the NV center in diamond [55].In Supplementary Figure S5, we present measurements of the antibunching rate of single NV centers in nanodiamonds as a function of excitation power; the results show clear linear scaling and a zero-power offset for γ 1 consistent with the expected optical lifetime.The PECS observations of h-BN's QEs in Fig. 4 defy this expectation.The power-scaling of γ 1 for QEs A, B, and C is clearly sublinear, with a saturation behavior (Model II or Model III) characterized by a steep slope at low power tapering off to a shallow slope at high power.Moreover, the γ 1 measurements for QE A are all less or equal to the measured spontaneous decay rate (dashed line in the upper-left panel of Fig. 4), whereas γ 1 always exceeds the spontaneous rate for a direct optical transition.The zero-power offset for γ 1 in QEs A and B is consistent with zero but poorly constrained due to the steep low-power slope; the offset is non-zero for QEs C, D, and E. QEs D and E exhibit linear power-scaling of γ 1 , however the range of available powers is smaller than for the other emit-ters, and we cannot rule out a saturation behavior at higher power.Previous studies of QEs in h-BN have revealed hints of power-independent antibunching rates [54] and nonlinear power scaling [3,26], however these observations were never satisfactorily explained.
The antibunching rate's nonlinear power dependence can be understood in the context of an indirect excitation mechanism, as illustrated in Fig. 5(a), where optical excitation leads to the population of multiple states: levels 2 and 3, with competing rates Γ 12 and Γ 13 .Indirect population of the radiative state (level 2) through such a mechanism creates a rate-limiting step (3 → 2) to the optical emission pathway (2 → 1) that leads to nonlinear scaling of the observed antibunching rate, as shown in Fig. 5(c).The rate-limiting nature of this process is intuitively obvious in the limit where Γ 12 /Γ 13 1.However, we find that the nonlinear saturation behavior remains qualitatively consistent independent of the pumping-rate ratio across a wide range of simulation settings where Γ 12 /Γ 13 ∈ [0, 2] (see Supplementary Figure S8).In the regime Γ 12 /Γ 13 1, the population of level 2 is still mostly determined by the indirect excitation pathway through level 3, with rate κ 32 , and the dominant antibunching rate saturates to a value close to κ 32 + Γ 21 .In the regime where Γ 12 /Γ 13 > 1, two underlying rates in the dynamical system are associated with the antibunching dip.As discussed previously, the eigenvalues associated with these rates can be real or complex depending on the relative magnitudes of transitions in the system.However, the fast rate associated with the direct population of level 2 and the subtle signatures of complex eigenvalues on the shape of the antibunching dip turn out not to be detectable when we include realistic assumptions for the experimental limits on timing resolution and shot noise.Instead, we observe a single effective antibunching rate γ 1 that exhibits nonlinear saturation similar to the slow rate.
The bunching dynamics observed in Fig. 4 can also be understood within our optical dynamics models by including metastable shelving states.In Fig. 5, the key qualitative difference between spontaneous population of the metastable state(s) (e.g., spin-dependent intersystem crossings) and optically pumped transitions (e.g., ionization/recombination) appears in the power-scaling and zero-power offset of the associated bunching rate.Spontaneous transitions are characterized by a rate with a non-zero offset that tends to saturate with increasing pumping power, whereas optically pumped transitions have zero offset and increase quasi-linearly.Previous studies considering the power-scaling of bunching rates for QEs in h-BN nano-flakes and exfoliated flakes have proposed similar optically-pumped models [3,26].Similar behavior has also been observed in color centers such as the silicon-vacancy center in diamond, attributed to power-dependent de-shelving from higher lying states to the metastable state [56].However, the heterogeneity and complexity of these processes for QEs in h-BN, both regarding the number of levels and the type of transitions, have not been considered before.
We observe both qualitative bunching behaviors in the data of Fig. 4, with several QEs exhibiting multiple bunching levels that apparently have different transition mechanisms.In some cases, individual bunching rates exhibit power scalings with features of both phenomena; for example, γ 2 for QEs A, B, D, and E appears to have a non-zero offset and yet increase linearly with power.This could indicate that the associated state can be populated both spontaneously and through an optically pumped pathway.Our simulations support this intuitive reasoning.For example, Supplementary Figure S9 shows the results of simulations of the same four-level system as in Fig. 5(a), but with rates chosen to reproduce the observations for QE A from Fig. 4. We indeed find that a combination of spontaneous and optically pumped transitions to the metastable state (κ 24 and κ 41 ) yields a bunching rate γ 2 with a non-zero offset that scales linearly with pumping power.Moreover, setting κ 32 < Γ 21 creates a situation where the spontaneous emission rate exceeds the observed antibunching rate, Γ 21 > γ 1 , over a wide range of pumping power, in agreement with our observations.The quantitative magnitudes of γ 1 , γ 2 , and the bunching amplitude, C 2 , are also reproduced by the simulations.This highlights the versatility of optical dynamics simulations as a valuable tool to recreate or predict optical dynamics based on complex combinations of radiative and nonradiative processes.To fully capture the observed dynamics of any particular QE, including the additional bunching rates γ 3 and γ 4 (where applicable), more metastable states are required in the simulations.We further note that the number of observed bunching timescales represents a lower limit on the number of metastable states, and hence some states could actually represent multiplets associated with different spin manifolds.Even with those caveats, these observations present the opportunity for quantitative comparisons with theoretical predictions.

C. Consistency with theoretical proposals
Several defect structures have been proposed as the origin of visible-wavelength single-photon emission in h-BN, including the boron dangling bond (DB) [43], [46], and V B C N [34].The negatively-charged boron vacancy, V − B , has been suggested to give rise to an optically detected magnetic resonance signal observed for emitter ensembles [12], however V − B has a ZPL of ∼1.7 eV and couples more strongly to phonons (S HR ∼ 3.5) [57], producing a PL band between 800-900 nm that does not overlap with the emitters considered here.Early studies highlighted V N N B as the potential origin of visible QEs [2], but recent calculations show that the coupling to phonons is substantially larger than observations [46].More recently, V B C N has been proposed based on the observation that carbon is correlated with the emission signal, but the calculated PL spectrum [34] does not match our observations.The V B C N calculations also predict a single, linearly-polarized absorption dipole, which is inconsistent with our measurements.The calculated PL spectrum and strain dependence of V N C B [46] are in reasonable agreement with the our observations.However, the optical transition for V N C B occurs in the triplet channel, while the calculated ground state is a singlet; the authors did not propose a mechanism through which the triplet channel is populated quickly enough to give rise to the optical emission they considered.
The boron DB is predicted to possess an optical transition at 2.06 eV with a Huang-Rhys factor of 2.3 [43], which is in close agreement with the values observed in this study.In addition, the variations in ZPL and S HR for the observed emitters can be explained by out-of-plane distortions [53].The ground state of the boron DB is a singlet, and the predicted existence of a triplet excited state can explain the presence of level 4 in Fig. 5(a).Another important feature of the boron DB model is the proximity of the states to h-BN's conduction band [43]; this allows electrons to be optically excited directly into the conduction band, depending on the excitation energy, explaining the misalignment of the absorptive and emissive dipole when the excitation energy is increased.Other proposed models do not provide an explanation for the misalignment.For instance, in the case of V N C B the optical transition occurs in the neutral charge state, and for the excitation energies considered here, photoionization will not occur [58].
Within the boron DB model, we would interpret level 3 in Fig. 5(a) as the conduction band and κ 32 as the nonradiative capture rate.To support this interpretation, we have estimated the relevant capture rate κ 32 of a photoionized electron from the conduction-band minimum into the DB excited state [level 2 in Fig. 5(a)].This capture rate is a product of a capture coefficient and the density of electrons in the conduction band.A first-principles calculation yields a capture coefficient of 4 × 10 −7 cm 3 s −1 (see Supplementary Materials, Section S9).The density of electrons is estimated based on the thermal velocity of the photoionized electron (∼ 10 5 m s −1 ) and a typical electron energy relaxation time of ∼ 1 ps [59].In the time it takes the electron to relax to the conduction-band minimum, it can thus travel ∼ 100 nm; this distance corresponds to an effective electron density of 2.4 × 10 14 cm −3 .Multiplying this value with the calculated capture coefficient gives a rate of κ 32 ∼ 100 MHz, in compelling agreement with the observed saturation antibunching rates of γ 1 ∼ 300-800 MHz for QEs A, B, and C. Our calculations also show that capture into the excited state is favored over capture into the ground state by more than 5 orders of magnitude, justifying the neglect of κ 31 in the general model of Fig. 5(a).
The inclusion of photoionization allows us to further rationalize the heterogeneity in bunching behavior of the observed emitters: the photoionized electron is not necessarily re-captured at the same QE, but may instead be captured by a neighboring defect, leaving the QE in a nonfluorescent, ionized configuration that likely requires optical excitation of additional free electrons to restore emission via subsequent electron capture.This process would be represented in Fig. 5(a) by an optically pumped transition, where level 4 represents an ionized state of the QE.The emitters may therefore be highly sensitive to the local defect environment.Unlike other proposed defect models, we conclude that the boron DB model is thus capable of explaining numerous aspects of the experimental observations, lending support to this proposed microscopic structure.

IV. CONCLUSION
The observations in this work reveal that h-BN's QEs have intricate electronic level structures and complex optical dynamics including multiple charge or spin manifolds.Our proposed electronic-structure models complement previous reports [3,26,54] and explain the quantitative features of our observations.In particular, the models explain the observation of nonlinear power-scaling of the antibunching rate as well as heterogeneous magnitudes and power-scaling behavior of multiple bunching rates.Whereas past reports have lacked consensus on mechanisms to explain the observed optical dynamics of h-BN's QEs, and many posited chemical and electronic structure models have failed to adequately explain the heterogeneous observations, we show that the boron dangling bond model is remarkably consistent with experiments, especially accounting for the role of local distortions, photoionization, electron capture, and the QEs' heterogeneous local defect environment.Future experiments should be designed to investigate these details, for example time-domain studies of transients associated with charge and spin dynamics, and temperature-and excitation-energy-dependent variations of the PL lineshape, vibronic spectrum, and polarizationdependent excitation cross section.Combined with theoretical models, such experiments can resolve the underlying transition rates and resolve the disparate influences of the QEs' intrinsic properties with those of their local environments.The observation of pure single-photon emission with g (2) (0) = 0 resolves earlier questions about h-BN's QEs [48], affirming their potential for use in photonic quantum technologies.More generally, we hope that the approach and techniques presented in this work -especially the quantitative use of PECS -present a model to formulate optical dynamics models for QEs in any material platform [47,60].Our models can be adapted to account for recent observations of magnetic-field-dependent optical dynamics [11] and optically detected magnetic resonance [13,14] in h-BN.Subsequently, they can be used to design protocols for initialization, control, and readout of quantum-coherent spin states for quantum information processing and quantum sensing.

A. Sample preparation
We used a confocal microscope (see Supplementary Materials) to isolate individual QEs in h-BN under ambient conditions.The h-BN samples were sourced from two batches (purchased ∼2 years apart) of bulk, single crystals from HQ Graphene.Each batch consisted of roughly 20 different individual crystals.The bulk crystals were mechanically exfoliated using a dry transfer process [49] resulting in thin (≤100 nm) and large area (∼10 µm) flakes of h-BN.The exfoliated flakes were transferred to a SiO 2 /Si substrate with micro-fabricated circular trenches 4 µm to 8 µm in diameter and 5 µm deep [6] using a dry transfer process.Table S1 in Supplementary Materials highlights the crystal from which the h-BN flake under study came.
Prior to the optical studies, the exfoliated h-BN samples were cleaned with a soft O 2 plasma (Anatech SCE 106 Barrel Asher, 50 W of power, 50 sccm O 2 flow rate) for 5 minutes to remove polymer residues resulting from the transfer process.The samples were then annealed in a tube furnace at 850 °C in low flow Ar atmosphere for 2 hours.Annealing h-BN has been found to brighten the emitters [35].While not the focus of study here, annealing for longer time (2 hours vs commonly used 30 minutes) appears to improve emitter stability.One sample was also exposed to a focused ion beam chamber operated in scanning electron mode (FEI Strata DB235 FIB SEM) but was not directly exposed to the electron beam.Supplementary Table S1 summarizes the three samples investigated, the QEs studied in each sample and the annealing treatment received by each sample.

B. Experimental details
Supplementary Figure S1 depicts a simplified schematic of the room-temperature confocal microscope used to measure the QEs.There are two available excitation sources: a 532 nm (green) cw laser (Coherent, Compass 315M-150) and a 592 nm (orange) cw laser (MPB Communications, VF-P-200-592).The power and polarization of each excitation path can be independently selected.Reported power values are measured just prior to the objective.In addition, a shutter completely blanks the excitation source when imaging is not in use to mitigate unnecessary light exposure.The excitation paths are combined with the collection path using a long pass (LP) dichroic mirror (Semrock, BrightLine FF560-FDi01 for green and Semrock, BrightLine FF640-FDi01 for orange).The LP dichroic cut-off is 560 nm for green excitation and 640 nm for orange excitation.A fixed half-wave plate in each of the excitation paths corrects for the birefringence induced by the dichroic mirrors.The co-aligned excitation and collection paths are sent through a 4f lens system with a fast steering mirror (Optics in Motion, OIM101) and a 0.9 NA 100x objective (Olympus, MPI Plan Fluor) at the image planes.This allows for the collection of wide-field, rastered, micro-photoluminescence (µ-PL) images.The objective is mounted on a stage system for changing the field of view.
The collection path consists of a linear polarizer (Thorlabs, WP25M-VIS) for measuring the emission polarization as well as a wide-band variable retarder (Meadowlark, LRC-100) which compensates for the birefringence induced by the dichroic.A LP filter specific to the excitation color fully extinguishes any scattered excitation light and the Raman signal.The cut-on wavelengths are 578 nm (Semrock, BLP01-568R-25) and 650 nm (Semrock, BLP01-635R-25) for green and orange, respectively.The filtered light is focused onto the core of a 50 µm core multi-mode fiber (Thorlabs, M42L01) acting as a pinhole.The output of the fiber is connected to a fiber switch (DiCon, MEMS 1x2 Switch Module) which can switch the collected emission to either a 50:50 visible fiber splitter (Thorlabs, FCMM50-50A-FC) or a spectrometer (Princeton Instruments, IsoPlane160 and Pixis 100 CCD).The outputs of the fiber splitter are sent to two identical single-photon counting modules (SPCM, Laser Components, Count T-100) resulting in a Hanbury Brown Twiss interferometer.The outputs of the SPCMs are either measured by a data acquisition card (National Instruments, DAQ6323) for general-purpose counting or a time-correlated single-photon counting (TCSPC) module (PicoQuant, PicoHarp 300) for recording the photon time-ofarrival information with a full system resolution of ∼350 ps.

C. Spectra analysis
The PL spectra are collected as multiple exposures and averaged after correcting for dark counts, cosmic rays and wavelength-dependent photon collection efficiency.The PL spectra are measured as a function of wavelength, λ and binned to determine spectral distribution function, S(λ).To analyze the vibronic coupling, the measured spectra must be converted to a form suitable for analysis with the general theory of electron-phonon coupling in three dimensional crystals [61,62].To do this, the spectral probability distribution function is obtained through where h is Planck's constant, c is the speed of light, and E is the photon energy.The emission lineshape, L(E) is derived from S(E) as which accounts for the photon-energy dependence of spontaneous emission.The emission lineshape is then fit following the method described in [6].From the fit, the following free parameters are determined: the ZPL energy, E ZPL , the ZPL Lortentzian linewidth, Γ ZPL , the Huang-Rhys factor, S HR , and the one-phonon vibronic coupling lineshape, approximated as an interpolated vector of values spanning the phonon spectrum in h-BN.The Debye-Waller factor, w DW , can be calculated from w DW = e −SHR .

D. Second-order photon autocorrelation function
For a given QE, all autocorrelation measurements are performed with the excitation polarization set at the angle of maximum excitation and the collection path has the polarizer removed.Due to the varying QE brightness, which affects the signal-to-noise ratio of the antibunching signal, measurements are integrated for 10 s to 140 min with repositioning occurring every 2 min.All errors from the fitting denote one standard deviation.
Due to timing jitter in the single photon counting modules (SPCMs) introducing systematic artifacts at short delay times, the data are analyzed in two stages: logarithmic and linear scales.First, the the autocorrelation data are binned over a log scale for visualizing the dynamics over 9 orders of magnitude in time (0.1 ns to 1 s) corrected for background [63] (see Supplementary Materials), and then fit by multiple instances of Eq. 5 with n = [2,5].The best fit, and corresponding n, is then determined by calculating the AIC and comparing the reduced chi-squared statistic.This method determines the number of bunching levels and their rates and amplitudes that best explain the observations.The QE's autocorrelation data are then binned over a linear scale that contains the antibunching features (τ ≤ 30 ns).To account for the timing jitter in the SPCMs, the instrument response function (IRF) is found by measuring the autocorrelation signal of an attenuated picosecond pulsed laser sent through the HBT interferometer and binned over the same linear scale as the QE.A convolution of the IRF with a modified Eq. 5 is fit to the backgroundcorrected data, given by where C B (τ ) is the total bunching contribution found from the log scale analysis (first step) and only C 1 and γ 1 are allowed to vary.The autocorrelation at zero delay is then given by which is used to determine the purity of single-photon emission from the QE.

E. Electronic level structure simulations
A four-level optical rate equation is used to model aspects of the observed autocorrelation data.The model is defined as where P is a vector of state populations P i and G is a generator matrix describing the transition rates and is given by where Γ 13 is the excitation rate, Γ 21 is the radiative emission rate, and κ ij are nonradiative rates that are either fixed or a proportion of the excitation rate.The g (2) (τ ) can be found by calculating the time evolution of the radiative state, P 2 given the system started in state P 1 and normalizing by the steady state population of P 2 .This is given by where τ = t 2 − t 1 .The differential equation (Eq.11) given the initial state is solved in MATLAB using the function ode15s.Timing resolution limitations and shot noise are added to the simulated autocorrelation function to best recreate the measurements.To model timing resolution, the simulated data are only analyzed for t 0 ≥ 0.5 ns.To include shot noise, a standard deviation, σ 0 , is set for the first delay time.
Assuming shot noise, this standard deviation is converted to mean number of photons as The log-scale processing results in the average number of photon correlations detected in each bin increasing linearly with the delay time, From this, a simulated, noisy g (2) (τ ) is calculated as where For the pumped bunching, κ 24 = 6 kHz/MHz × Γ 13 and κ 41 = 3 kHz/MHz × Γ 13 .

F. Theoretical calculations
We perform first-principles density-functional theory calculations as implemented in the VASP code [64,65].We utilize the hybrid functional of Heyd, Scuseria, and Ernzerhof [66,67] to ensure accurate energetics, electronic structure, and atomic geometries.The fraction of Hartree-Fock exchange is set to 40%, consistent with previous studies [43,53].A plane-wave basis with projector augmented-wave potentials [68] is used, and the energy cutoff for the basis is set to 520 eV.
The boron dangling bond is modeled in a 240-atom supercell with volume 2110 Å3 within periodic boundary conditions [69].A single, special k-point (0.25, 0.25, 0.25) is used to sample the Brillouin zone.Lattice vectors are held fixed while the atomic forces are relaxed to below 0.01 eV/ Å.To calculate the nonradiative capture coefficient, we utilize the formalism of Ref. [70] implemented in the Nonrad code [71].  Excitation (θ ex ) and emission (θ em ) dipole orientation, excitation (V ex ) and emission (V em ) visibility.

Simultaneous Fit of Emission and Excitation Polarization of QE C for 532 nm Excitation
The emission and excitation polarization of QE C is independently modeled using Eq. 2 and 3 to estimate its optical dipole orientation and visibility.Further, the misalignment of the emission and excitation dipole is estimated from the independent fit results using Eq. 4.However, the low visibility of excitation polarization for 532 nm excitation suggests superposition of multiple excitation dipoles which means the excitation pathways could be via multiple dipoles, one of them oriented along the emission dipole.To check this hypothesis, the following empirical equations are simultaneously fit to the emission and excitation polarization data: where A is the normalized amplitude, B is the normalized offset, θ 1 is the emission dipole orientation and orientation of one of the two excitation dipoles and θ 2 is the second excitation dipole.
Figure S2 shows the normalized emission and excitation polarization data and the resultant simultaneous fits.Table S3 presents the fit result.The emission and first excitation dipole is oriented 111.22 o .This agrees with the independent fit to the emission polarization using Eq. 2 which gives dipole orientation of 111.09 o (Table S2).This proves the hypothesis that the excitation polarization is a superposition of two excitation dipoles, one aligned with the emission dipole and the other misaligned by ∼60 o .Collectively, the two excitation dipoles result in an effective dipole with low visibility.Figure S3 illustrates the second-order autocorrelation function calculated and analyzed over logarithmic time scale for each emitter and fit using an empirical model (as discussed in the main text), g (2) (τ ) = 1 − C 1 e −|t|/τ1 + n i=2 C i e −|t|/τi .The resultant fit parameters are summarized in Table S4.

S5.2 Background Correction
Background correction is done to account for the background and incoherent light detected along with the signal which can affect the autocorrelation function.The following two background correction techniques are used: 1. Recording background from an offset spot: The autocorrelation data was acquired from a background spot, same as the emitter.The background spot is an offset spot, ∼1 µm from the emitter which seems to emulate the true background.
The background data was acquired for the equivalent time as the emitter, with all other experimental conditions such as excitation power kept same.From the background data, the average background count rate was determined.This technique was applied to QEs A and B.
2. Recording background from a 2D scan: Instead of recording background data of an offset spot, X and Y line (µ-PL) scan of the emitter were acquired.Since the emitter is tracked during the acquisitions, the µ-PL line scans along X and Y pass through the center of the emitter.A 2D Gaussian fit to the line scans provides the background and signal of the spatially isolated emitter.This technique speeds up the data acquisition by a factor of 2 since the data from an offset spot is no longer needed to be acquired.This technique better approximates the background since the estimation is done right around the emitter, instead of an offset spot.This technique was applied to QEs C, D and E, and adopted as future autocorrelation background correction technique.
Using the background and signal acquired, following background correction equations from Ref. [1] are used to determine background-corrected autocorrelation function: where g (2) (τ ) is determined from the QE and g bkgd (τ ) is the background-corrected autocorrelation function.

S5.3 Akaike Information Criterion
For a given set of data, the Akaike Information Criterion (AIC) is a quantitative method to determine the relative quality of a collection of statistical models.Using the AIC, the relative quality of each of the models compared can be estimated.It can be applied to any data and fitting routines.Thus, it can be used to determine which model best fits to a given data set.However, it does not determine absolute quality.Using the AIC, the relative quality of Eq. 5 for different n is determined.The likelihood of a model (n = [2,5]) explaining the actual data is determined by comparing the AIC of each model.
The AIC is defined as: where p is the number of fit parameters in the model.Here, Gauss log likelihood estimation is used.The Gauss log likelihood is defined as: where −2 ln (L G ) is the Gauss log likelihood, N is the number of data points, c i is the measured value of the i th data point, m a is the model function, m a (x i ) is the model predicted value of the i th data point, x i is the residual of the i th data point and σ i is the standard deviation of the i th data point.Thus, The likelihood of a model to explain the actual data from a collection of models is determined by calculating the weight of the model, defined as: where w i is the weight of the i th model, AIC min is the minimum AIC value among all models and AIC i is the AIC value of the i th model.The model with highest w best explains the actual data, amongst the collection of models.However, for relatively close values of two w i , a simpler model could be selected.Thus, the AIC is used to determine relative quality of n for n = [2,5].

S6. Lifetime Measurement (QE A)
Figure S4 shows time-resolved PL of the QE A. Time-resolved PL is acquired by using a pulsed laser (NKT Photonics, Fianium Whitelase) with the excitation wavelength centered ∼580 nm and 40 MHz pulse rate.The PL is recorded in histogram mode with the TCSPC module (PicoQuant, PicoHarp 300).The lifetime is obtained by fitting the convolution of the IRF and an empirical model (A × exp (−t/τ ) + B) to the data, given as IRF * (A × exp (−t/τ ) + B).The IRF is obtained in same configuration as the time-resolved PL of the QE.The lifetime is measured as 2.82±0.004ns (355 MHz).

S7.1 Nitrogen Vacancy Center in Diamond
The antibunching rate as a function of excitation power is typically linear, where the y-intercept provides the inverse of the radiative lifetime and the slope relates the pumping power to excitation rate.This is the case for the NV center in diamond.The autocorrelation function of single NVs in multiple nanodiamonds is calculated to verify linear dependence of antibunching rate on excitation power.Figure S5 shows the antibunching rate, γ 1 as a function of excitation power of NV center in two nanodiamonds.The nanodiamonds were dropcasted on a silicon wafer and probed for single NVs.The nanodiamond sample was studied in the same setup as the h-BN samples discussed in the main text.The data acquisition and analysis is as discussed in the main text (no IRF correction was performed).The data is acquired using green (532 nm) excitation.To check for linear dependence, the rates are fit using Eq.6a.The y-intercept is calculated to be 51 MHz (NV in ND1) and 67 MHz (NV in ND2), which agrees with that reported in the literature [2].

S7.3 Antibunching Amplitude
Figure S7(a-e) presents the excitation power dependence of the antibunching amplitude, C 1 of the QEs. Figure S7(f) presents the simulated antibunching amplitude for spontaneous and pumped transition mechanisms.The dashed lines are fits using the first order saturation model (Eq.6b).The fit results are shown in Table S6.  Figure S9 shows the result of simulations to recreate the observed optical dynamics of QE A. The same model discussed in the main text is used in the simulations (Fig. 5(a)).The radiative and nonradiative rates are chosen such that the resultant dynamics best recreates the observed optical dynamics of QE A (1 st column of Fig. 4).There are two important qualitative features of simulating QE A that differ from the simulations presented in Fig. 5 of the main text.The first is that when κ 32 < Γ 21 , we find that the observed antibunching rate, γ 1 , is less than the spontaneous emission rate, Γ 21 , over a wide range of power settings.Here, Γ 21 = 300 MHz and κ 32 = 60 MHz.The second important feature of simulating QE A is that the nonradiative transition mechanism involves both the spontaneous and optically pumped components.The nonradiative rates κ 24 and κ 41 are given the following form:

S7.4: Simulations
where κ ij,0 is the spontaneous emission rate, and β ij is a scaling factor for the optically pumped transition, giving the corresponding transition rate at saturation when Γ 13 = Γ 21 .Here, we set κ 24,0 = 24 kHz, κ 41,0 = 18 kHz, β 24 = 9 kHz, and β 41 = 3 kHz.The combination of spontaneous and pumped transition quantiatively reproduces the non-zero offset of the bunching rate and the quasi-linear power scaling.We calculate the nonradiative capture coefficient C n for the capture of an electron from the conduction band into the boron dangling bond using the formalism implemented in the Nonrad code [3].We will focus on the ground state [level 1 in Fig. 5(a)] and the optically active excited state [level 2 in Fig. 5(a)] of the dangling bond, which are separated by 2.06 eV [4].In equilibrium, the dangling bond is in the negative charge state and is occupied by two electrons.When the excitation energy is sufficiently large, an electron can be excited into the conduction band and the dangling bond is photoionized, changing the charge state from negative to neutral [process Γ 13 in Fig. 5(a)].Subsequent re-capture of this electron returns the dangling bond to the negative charge state.We consider two potential scenarios for this nonradiative process mediated by electron-phonon coupling: (1) The electron is captured directly into the ground state of the dangling bond, with rate κ 31 [not depicted in Fig. 5(a)], or (2) the electron is captured into the excited state of the dangling bond [κ 32 in Fig. 5(a)].Process (2) puts the defect in the optically active excited state, from which a photon can then be emitted, with an emissive dipole unaligned with the absorptive dipole.
To evaluate the nonradiative capture coefficients, we extract several parameters from our density-functional theory calculations: the transition energy, the mass-weighted root-mean-square difference in atomic geometries ∆Q, the phonon frequencies in the initial (i) and final (f ) states Ω i/f , and the electron-phonon coupling matrix element W if .The transition level for capture into excited state, which is used to determine the transition energy, is above the conduction-band minimum, while the single-particle states are in the gap.For the purposes of our capture coefficient evaluation, we shift the transition energy to be consistent with the 200 meV difference observed experimentally [5]; we verified that the conclusions are insensitive to the choice of the energy shift.The degeneracy factor in the nonradiative rate [3] is set to 1 since the dangling bond does not possess any configurational degeneracy.A scaling factor that accounts for charged defect interactions [see Sec.III.E. of Ref. [6]] is not necessary in this case because capture occurs in the neutral charge state and the electron-phonon coupling is evaluated in the neutral charge state.
At room temperature, we calculate C n for capture into the ground state to be 1.2 × 10 −12 cm 3 s −1 and into the excited state to be 1.2 × 10 −7 cm 3 s −1 .We can thus safely assume that capture into the excited state will dominate.These capture coefficients are larger than typical radiative capture coefficients, which are on the order of 10 −13 -10 −14 cm 3 s −1 [7], justifying our implicit assumption of nonradiative rather than radiative capture.
Previous work has already demonstrated that out-of-plane distortions are important for understanding both the symmetry [4] and transition rates [8] of the boron dangling bond.Here we include the effect of out-of-plane distortions on the capture coefficient, following the methodology of Ref. [8].In short, a plane neighboring the dangling bond is bent to create a "bubble", with the height of the bubble being referred to as the applied distortion h.The bubble induces an out-of-plane distortion in the dangling bond and allows us to study its effect on the capture coefficient.The influence of the out-of-plane distortion on the calculated parameters is shown in Fig. S11.For comparison purposes, we will use the average at room temperature 4 × 10 −7 cm 3 s −1 as a representative value for capture into the excited state.The lines are a quadratic fit to the calculated parameters and are intended to guide the eye.The calculated (e) electron capture coefficient at 10K (solid), 300K (dashed), 600K (dashed dotted), and 900K (dotted).The parameters for capture into the ground state are shown in blue and the excited state are shown in orange.

FIG. 1 .
FIG. 1. Photoluminescence characterization.In all panels, data plotted in green (orange) were acquired under 532 nm (592 nm) excitation.(Column 1) µ-PL images of the QEs (circled), acquired under 592 nm (QEs A-D) or 532 nm (QE E) excitation.Scale bars denote 1 µm.(Column 2) Second-order photon autocorrelation function (colored points), fit using an empirical model discussed in the text (black curve).Error bars represent Poissonian uncertainties based on the photon counts in each bin.(Column 3) Steady-state, background-subtracted PL intensity as a function of excitation power (points), fit using an empirical saturation model discussed in the text (solid curves).Error bars represent one standard deviation based on three measurement repeats.(Column 4) PL spectra and polarization data.Vertical colored lines represent the excitation laser wavelengths, and black dotted lines indicate cut-on wavelengths for long-pass optical filters in the collection path.Insets: PL intensity as a function of linear excitation polarization angle (colored circles) or filtered by linear polarization angle in emission (black squares).Solid curves are fits to the data using an empirical model discussed in the text.

FIG. 2 .
FIG. 2. PL emission lineshapes.(Left Column)Observed spectral emission lineshapes (black points) and fits according to Huang-Rhys theory (blue curve).Experimental uncertainties are comparable to the size of the data points.(Right Column) One-phonon vibronic coupling lineshape corresponding to the fits at left.

FIG. 4 .
FIG.4.Photon emission correlation spectroscopy.Excitation power and wavelength dependence of (Row 1) the antibunching rate (γ1) denoted as circles, (Row 2) the bunching rate (γ2) denoted as circles, (Row 3) the bunching amplitude (C2) denoted as circles, (Row 4) the bunching rate (γ3) denoted as squares and (Row 5) the bunching amplitude (C3) denoted as squares, of the five QEs presented in each column.The black dashed line (QE A) represents the lifetime.Orange (green) data correspond to excitation at 592 nm (532 nm).The error bars represent one standard deviation.The dotted lines are fits to the metadata as discussed in the text.
Figure S1.Experiment Setup.A simplified version of the room temperature optical setup showing the essential optical and electronic components used to probe the QEs in h-BN.The green dashed line represents the 532 nm (green) excitation path which can be switched to 592 nm (orange) excitation.

2 Figure S2 .
Figure S2.Quantum emitter C: emission and excitation polarization with green excitation.The normalized emission polarization is shown with blue circles and normalized excitation polarization is shown with orange circles.The green and yellow curves are simultaneous fits to the normalized emission and excitation polarization, respectively.

Figure S4 .
Figure S4.Quantum emitter A: lifetime measurement.Time-resolved PL (light orange circles) and fit (dark orange circle) to the data.The error bars represent one standard deviation.

Figure S5 .
Figure S5.Antibunching rate of NV centers in nanodiamonds.The antibunching rate (denoted by circles and squares) is measured as a function of excitation power.The lines (dashed and dotted) are fits to the rates.The error bars represent one standard deviation.

4 Figure S6 .
Figure S6.Quantum emitter A: γ4 and C4 for orange excitation.The error bars denote one standard deviation.

Figure
FigureS8presents result of simulations for Γ 12 /Γ 13 = [0, 2], for the case of spontaneous transition via level 4. As highlighted in the main text, the result of simulations for different Γ 12 /Γ 13 are qualitatively similar.

Figure S7 .
Figure S7.Antibunching amplitudes.(a-e) Antibunching amplitudes (circles) of QEs A to E. The error bars represent one standard deviation.The dashed lines are fits using an empirical model discussed in the text.(f) Simulated antibunching amplitudes as a function of excitation rate for spontaneous (circles) and pumped (squares) transition mechanism discussed in the main text.

2 Figure S8 . 2 Figure S9 .
Figure S8.Simulation of optical dynamics for different Γ12/Γ13.Antibunching and bunching parameters resulting from simulation of the model discussed in main text for Γ12 as a factor of Γ13 for spontaneous transition.Simultaneous excitation to level 2 and level 3 takes place at different rates which is ratio of the two rates.

Figure S11 .
Figure S11.Calculated parameters and capture coefficient for capture of an electron from the conduction band into the neutral boron dangling bond, as a function of applied distortion h.Calculated (a) transition energy, (b) mass-weighted root-mean-square difference in atomic geometries, (c) initial (up triangle) and final (down triangle) phonon frequencies, and (d) electron-phonon coupling matrix elements.The lines are a quadratic fit to the calculated parameters and are intended to guide the eye.The calculated (e) electron capture coefficient at 10K (solid), 300K (dashed), 600K (dashed dotted), and 900K (dotted).The parameters for capture into the ground state are shown in blue and the excited state are shown in orange.

Table S1 .
The Samples, Quantum Emitters and Sample Treatments.

Table S3 .
Result of Simultaneous Fit of QE C's Emission and Excitation Polarization.

Table S4 .
Logarithmic Time Scale Autocorrelation Function Fit Parameters.