Photophysics of quantum emitters in hexagonal boron-nitride nano-flakes

Quantum emitters in hexagonal boron nitride (hBN) have attracted significant interest due to their bright and narrowband photon emission even at room temperature. The wide-bandgap two-dimensional material incorporates crystal defects of yet-unknown configuration, introducing discrete energy levels with radiative transition frequencies in the visible spectral range. The commonly observed high brightness together with the moderate fluorescence lifetime indicates a high quantum efficiency, but the exact dynamics and the underlying energy level structure remain elusive. In this study we present a systematic and detailed analysis of the photon statistics recorded for several individual emitters. We extract the individual decay rates by modeling the second-order correlation functions using a set of rate equations based on an energy level scheme involving long-lived states. Our analysis clearly indicates excitation-power-dependent non-radiative couplings to at least two metastable levels and confirms a near unity quantum efficiency. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
As part of the growing family of two-dimensional van der Waals materials, hexagonal boron nitride (hBN) has received great interest as a building block for layered nano-electronic devices [1,2], mainly as a wide band gap substrate to facilitate for instance high performance of graphene electronic and opto-electronic properties [3,4]. With the discovery of optically active defect centers with high brightness and a narrow emission spectrum, hBN has emerged as a promising solid-state resource of single photons for applications in nano-photonics [5][6][7]. Single-photon emitters have been reported in single-layer, multi-layer and bulk hBN material [5,8], and emitters can be created by electron or ion irradiation and annealing [6,9]. The main emission peak, often referred to as the zero-phonon line (ZPL), is inhomogeneously distributed over a wide spectral range from around 500 nm up to 800 nm [10] and in comparison, the ZPL contribution to the total spectrum was found to be reduced in high purity crystals [11]. The bandwidth of the ZPL transition was recently shown to be nearly lifetime limited at cryogenic temperatures [12,13]. Strong dependence of the optical activity of a hBN emitter on the pump laser wavelength has been reported [14,15], and the wavelength dependence of the transition rates has been applied to resolve hBN emitters with super-resolution [16]. Optical absorption and emission polarization measurements of hBN emitters [12,17] suggest an indirect excitation mechanism depending on the excess energy of the excitation laser with respect to the ZPL [18]. Two types of defects have been identified and are denoted type I and II, where the former features a relatively broad ZPL and phonon side band (PSB) separated by roughly 160 meV, and the latter features a narrower ZPL and a suppressed PSB [6,19,20].
Fluorescence intermittency, commonly referred to as fluorescence blinking, and spectral diffusion of quantum emitters in hBN have been observed [11,14,19,21] and provide information about the emitter dynamics and their interaction with the local environment, e.g. surface charges.
Surface passivation of the substrate was shown to stabilize and narrow the emission line for type-II defects in multi-layer chemical-vapor-deposition grown hBN [19]. These observations suggest that a large range of time scales on the order of the interaction strength needs to be covered for characterization and that the defects are highly affected by their local environment due to the two-dimensional geometry of the host material and the resulting close proximity of the emitter to the surface.
Characteristic time scales for fluorescence intermittency may range from a few ns up to several seconds, indicating complex emitter dynamics as well as the presence of at least one additional state contributing to the level scheme. Very limited information about the emitter dynamics may be obtained from fluorescence count data: While the shortest time scale accessible is determined by the maximum photon count rate, binning over longer time scales averages the fluctuations. In contrast, the second-order correlation function g (2) (τ) contains dynamical information about the emitter level scheme spanning many orders of magnitude in time. Precise modeling of g (2) (τ) using an underlying level scheme may then yield the exact coupling strength between the emitter levels. This method has been used to study the photophysics of other systems such as quantum dots [22] and color-centers in diamond [23][24][25]. Measurements on hBN emitters show photon bunching with correlation times up to the ms-range, which clearly reveals the presence of long-lived metastable states [7,11,26], and their dependence on temperature has been investigated [10]. In those studies, the correlation functions have been empirically fitted with exponential functions, which allows for the extraction of the overall time scales but does not yield individual transition rates between the levels.
In this article we present a systematic study of the photon statistics emitted by individual defects in hBN nano-flakes. In particular, we focus on measurements of the second-order correlation function with delay times from sub-ns to the ms time scale at varying excitation power. We present data of four bright stable emitters and demonstrate the power dependence of their transition rates obtained by simulating the system dynamics with a set of rate equations. The extracted coupling strengths allow us to estimate the emitter quantum efficiency, which is found to be excitation-power dependent. Our experiments reveal complex dynamics with photon bunching up to the ms time scale and requiring a rate-equation model involving at least four levels. For some emitters and excitation powers we observe fluorescence intermittency on a timescale of seconds with distinct fluorescence-intensity levels. For these intensity levels we extract the individual correlation functions and demonstrate changes of distinct transition rates which we correlate with the observation of spectral diffusion.

Sample preparation and experimental setup
Our samples were prepared by drop casting 20 µL of an ethanol/water solution containing hBN nano-flakes (purchased at Graphene supermarket) onto a plasma cleaned quartz disc. The sample was then left to dry in ambient clean conditions before annealing for 30 min in a tube furnace at 850 • C with a controlled temperature ramp and cooling for several hours.
The sample was characterized in a confocal microscope with an immersion-oil objective (x100, NA=1.4). For excitation we used light from one of two 532-nm lasers focused to a diffraction-limited spot: a continuous-wave (cw) laser or a laser generating 4-ps-wide pulses with a repetition rate of 5 MHz. The detection setup consists of two avalanche photo diodes (APDs, Perkin Elmer, SPCM-AQR-14) separated by a symmetric beam-splitter forming a Hanbury Brown and Twiss (HBT) setup. During characterization of individual emitters, we used a notch filter at 532 nm and a combination of a long-and a short-pass filters in order to limit background signal and possible signal from nearby emitters. The following wavelength ranges have been used to characterize four selected emitters (denoted as emitters A-D): A: 550-785 nm, B: 647-700 nm, C: 647-785 nm, D: 600-785 nm.
The continuous and long-time acquisition of photon counting events from the APDs was carried out using a time-tagging module (PicoHarp 300, PicoQuant). From the raw data we subsequently calculated the second-order correlation function g (2) (τ) by numerically correlating each photon arrival event with all the other registered photons. We stress that the usually applied "start-stop" method gives a histogram that deviates significantly from the g (2) (τ) for times τ R −1 , where R is the photon detection rate [27].
The time resolution of our detection setup is ∼ 0.8 ns, which is the width of the instrumentresponse function (IRF). We measured the IRF by detecting lifetime of gold fluorescence, which is in the femtosecond time range. Photoluminescence spectra were recorded with a spectrum analyzer (Shamrock 500i, Andor) with a 532-nm notch filter before the spectrometer. All measurements were carried out at room temperature and air atmosphere.

Sample characterization
We report on the detailed characterization of in total four emitters, all from the same sample, which throughout this study are labeled emitters A-D. The confocal image presented in Fig. 1(a) shows a 10 × 10 µm 2 scan centered around the well-isolated emitter A. The normalized spectra of the four emitters are plotted in Fig. 1(b) with the x axis offset relative to the energy of the ZPL. We identify the ZPL as the highest in energy line in a range of periodic peaks corresponding to zero-, one-, and multi-phonon transitions. The ZPL is also the strongest line in agreement with previous works. We attribute any additional signal at higher energy to a background or distant emitters, as it is considerably less intense and separated by a random energy interval from the ZPL. We find that the one-phonon transition is red-shifted by 160.3 ± 4.1 meV (statistics taken for a total of 18 emitters including A-D) from the ZPL frequency E ZPL , which is in agreement with reports on type-I emitters [6,11,19]. The PSB magnitude of emitter C is roughly just ∼ 6% of the ZPL intensity which is also narrower compared to the other emitters [19]. The ZPL wavelength is indicated in the Fig. 1(b) legend and distributed over a ∼ 100 nm range. The mean count rateĪ for each of the emitters A-D as a function of cw excitation power P is presented in Fig. 1(c). The mean count rate was fitted with a saturation function I = I ∞ P P+P sat +bP+c, revealing the maximum count rate I ∞ for P − → ∞, the saturation excitation power P sat , a linear contribution b accounting for fluorescence background and a constant background c. In Table 1 we summarize these parameters for all four emitters and the total decay rate k tot obtained by fitting lifetime measurement data with a single-exponential function. All emitters have a significant brightness with I ∞ in the ∼ MHz range and relatively small background contributions bP sat at the saturation power P sat . A lower bound of 1 kHz is set for the constant background c to account for the dark-count level of the APDs and residual background noise from the setup. Throughout the investigation, emitter D bleached at an optical excitation power above 7 mW. Table 1. Summary of parameters for emitters A-D. The total decay rate k tot = 1/t LT is determined by low-power lifetime measurements, where t LT is the excited-state lifetime obtained by fitting measurement data with a single-exponential function. The maximum photon count rate I ∞ , the saturation power P sat , and the power-dependent background b are determined from fluorescence saturation measurements. The constant pump-independent background c was set to 1 kHz during fitting. Emitter

Second-order correlation function modelling
From the g (2) (τ) function measurements of all emitters we observe photon anti-bunching at vanishing delay times τ − → 0 and significant photon bunching on time scales between τ ≈ 10 −9 − 10 −3 s [e.g., see Fig. 3(a) for data of emitter A]. At times τ > 10 −3 s, the g (2) (τ) data flattens and approaches unity in the limit of τ − → ∞, demonstrating statistical independence of photons separated by large time intervals. The observation of photon bunching in g (2) (τ) over timescales up to the millisecond range reveals the presence of other transitions and levels than the radiative transition. Furthermore, examination of the tail of the correlation function after the bunching peak identifies clearly non-single-exponential decay, which implies that more than one metastable shelving state is involved in the dynamics of the emitter [28]. We therefore model the system with a set of four levels as indicated in Fig. 2, including a ground state (level 1), an optically excited state (level 2) radiatively coupled with k rad to the ground state, and two metastable shelving states (levels 3 and 4). With an optical excitation, the population is non-resonantly transferred from the ground state (level 1) to the optically excited state (level 2). The excited state may decay via a radiative transition back to level 1 and emit a photon or non-radiatively to one of the shelving states (levels 3. . . N). The shelving states subsequently relax to level 1. The diagram shows the possible transition pathways together with the corresponding transition rates. The manifold of vibrational sublevels of the ground and excited states (see details in text) are shown as gradients.  Figure 2) obtained from the fits of the measured correlation functions with Eq. (2). Solid lines are the linear fits to the power dependencies. The obtained linear trends are described by k 23 = 1.2(2) MHz P P sat + 2.0(2) MHz, k 31 = 1.07(3) MHz P P sat , k 24 = 5.4(1) MHz P P sat , k 41 = 6.0(2) MHz P P sat , and k exc = 211(7) MHz P P sat . The error bars show one standard deviation for the fitting parameters.
It is important to note that we did not obtain satisfying fits with small and randomly distributed residuals using a model with only one shelving state, which is in contrast to previous reports on emitters in hBN. We observed that for high optical excitation powers, addition of the second shelving state resulted in substantial (at least by half) decrease of the chi-square goodness-of-fit value, whereas the addition of further shelving states did not noticeably improve the goodness of fit (less than 2% change in the chi-square value). For lower optical excitation powers, the significance of adding the second shelving state was less pronounced, and for the lowest optical powers it did not play a role at all. This indicates that the dynamics related to the second shelving state becomes involved at higher optical excitation powers, which we analyze in more detail in section 3.
By means of optical excitation, population is transferred from the ground to the excited state with a rate k exc . The excited state may decay radiatively to the ground state or non-radiatively to any of the shelving levels at a rate k 2i , and the levels i may subsequently relax to the ground state at a rate k i1 . The actual ground and excited states consist of a manifold of vibrational states. Thus, under non-resonant excitation applied in our experiment, the emitter population is transferred from the ground state to a vibrational sublevel of the electronically excited state. However, at room temperature electron-phonon interaction leads to (i) quick vibrational relaxation due to loss of electron energy to phonons [29] and (ii) dephasing due to a modulation of the electron transition frequency by the lattice vibrations [30]. Both processes occur over timescales of 10 −14 − 10 −12 s which is confirmed by the spectral linewidth of the ZPL and the sharp initial onset of photo-counts [31] in the lifetime-measurement curve [see Fig. 5(b) later in the text]. The latter is almost instantaneous and fully determined by the IRF width of our setup. Slow intermediate relaxation would lead to a narrower ZPL linewidth or a pronounced slow rise of counts in the lifetime trace. It is interesting to note that the ZPL is the brightest line in spectra of all hBN emitters [ Fig. 1(b)]. From the Franck-Condon principle we infer that the vibrational potential of the emitter is nearly unchanged during electronic transitions, which is a strong indication of a high defect symmetry.
Since electron-phonon interactions are much faster than the luminescence lifetime of hBN emitters, transitions between vibrational sublevels and coherences between electronic levels may be neglected and hence the system dynamics can be summarized by the set of rate equations: where ρ i is the population of level i with the normalization i ρ i = 1. We assume that hBN emitters are single-photon quantum emitters and hence can only emit one photon at a time. To derive an expression for g (2) (τ) we also initially assume no background noise in the emission, which will be accounted for later. In this case, g (2) (τ) is related to the joint probability for the detection of a photon at time t and the detection of a subsequent photon at some later time t + τ. The detection of a photon at time t ensures that the emitter is prepared in its ground state (ρ 1 = 1). g (2) (τ) then follows the conditioned time evolution of time-averaged detected photocounts at time t + τ and hence the population of the excited state [28,32]. We can therefore use the solution of coupled differential Eq. (1) with initial condition ρ 1 (τ = 0) = 1 to obtain a model for the second-order correlation function via g (2) mod (τ) = ρ 2 (τ)/ρ 2 (τ − → ∞). In order to account for the background noise, we assume it is Poisson distributed and introduce a noise parameter σ = S/(S + B) as the ratio of signal level S to the total photon counts including background S + B. Then the measured correlation function can be written as [33]: The value of σ is found by fitting experimental data with Eq. (2) and is then compared with the saturation measurements (Table 1) to assess whether the experimentally detected noise level can explain the deviation from the expected condition g (2) (0) = 0 for a single emitter. Thus, the measured correlation function for a particular power is characterized by seven parameters: six transition rates k ij and the parameter σ. These are too many to describe the shape of the measured g (2) (τ) function, and we indeed found several distinct sets of parameters, which fitted experimental data equally well.
In order to decrease the number of fitting parameters, we jointly fit the entire data set of measured correlation functions at different excitation powers using the following assumptions: (i) the excitation rate k 12 is proportional to the excitation power P: k 12 = βP, (ii) the radiative rate k rad = k 21 is independent of power P and therefore is the same through all the experimental data sets, (iii) the emitter lifetime is measured in the regime of minimal excitation power and its value can be used to exclude the rate k rad from the fitting parameters for the lowest excitation power (and hence for all powers): 1/t LT = k tot = k rad + k 23 + k 24 . This leaves us with five fitting parameters {k 23 , k 31 , k 24 , k 41 , σ} for each data set and one joint parameter β for all data setsoverall approximately 5.1 parameters per data set. Since the parameter σ has low impact on the time dependence of g (2) (τ) and mostly determines its depth at τ = 0, there are only approximately four fitting parameters to describe the time variation of the measured correlation function. We thus exploit the correlations between measured g (2) (τ) functions at different excitation powers to minimize the number of fitting parameters and reliably fit the data using the four-level model.

Transition rates and photon bunching
A subset of measured g (2) (τ) functions of emitter A is presented in Fig. 3(a) for a range of excitation powers together with the curves obtained by fitting Eq. (2) to the data. We observe photon anti-bunching for all emitters over the entire range of excitation powers at a time scale τ 1 ns and photon bunching on time scales up to around 100 µs. One can see that Eq. (2) fits the experimental data very well except in a narrow region around τ = τ d corresponding to the zero time delay between channels, where the experimental data points form a rounded shape. The reason for this rounding is the limited time resolution of our setup (∼ 0.8 ns), which prevents resolving fully the narrow antibunching dip. The value of the correlation function at antibunching, g (2) (τ d ) ≈ 0.1 < 0.5, is the signature of single-photon emission. The deviation from an ideal anticorrelation [g (2) (τ d ) = 0] is fully consistent with the saturation measurement and the background counts extracted therefrom [ Fig. 1(c) and Table 1].
For all emitters we observe that bunching, quantified by the maximum g (2) (τ) value, decreases with increasing excitation power. As described in the previous section, we fit the g (2) (τ) functions of all emitters with a four-level model to obtain satisfying fitting results with homogeneously distributed residuals. This means that there are two metastable shelving states contributing to the dynamics of the emitters.
The rates extracted from the rate equation model for emitter A are summarized in Figs. 3(b) and (c) with the error bars illustrating one standard deviation. The solid lines are obtained by empirical linear fits to the data. All rates have a power dependence sufficiently well described by linear functions (R-squared > 0.98) and approach zero at vanishing excitation power, except k 23 with a finite value k 0 23 = 2.0 MHz in the limit P → 0. All scaling parameters are summarized in the caption of Fig. 3. Very similar results have been obtained for emitters B-D and are briefly presented in Appendix A. The observed linear dependencies of the deshelving rates k i1 with excitation power are in agreement with earlier observations on single defects in hBN [11] and diamond [34]. A model with linear scaling of also the shelving rates k 2i with excitation power has been successfully applied to describe the photophysics of SiV centers in diamond [24].

Fluorescence intermittence
For emitter A at an intermediate range of excitation powers between approximately 0.4 P sat and 1.5 P sat and after finalizing the characterization described in section 3.1, we observed telegraph-like fluctuations of the photon counts both under CW and pulsed laser excitation, clearly identifying bright and dark florescence levels. In the following, we will analyze the underlying dynamics related to this observation in more details using the model developed in section 2.3.
Three traces were recorded under CW excitation for a duration of 2 min each and for excitation powers of 0.44 P sat , 0.69 P sat , and 1.52 P sat . The fluorescence trace at 0.44 P sat is presented in Fig. 4(a). We processed the data by setting a threshold on the count rate and separated the counting events belonging to the bright and the dark fluorescence level, as indicated by the distinct blue (bright) and red (dark) color coding of the data. The g (2) (τ) functions were then calculated from experimental data for the bright and the dark fluorescence levels individually and each result fitted with Eq. (2).
The resulting g (2) (τ) traces for the bright and dark fluorescence levels are shown in Fig. 4(b) for the recording at 0.44 P sat . While the overall shape and the time scales are similar in both cases, photon bunching is significantly more pronounced for the dark fluorescence level. This behaviour is found for all three excitation powers we have investigated. Indeed, compared to the continuous acquisition (see Fig. 3), we observe an increase of k 23 [Fig. 4(c)], which explains the reduced counts recorded for the dark fluorescence level. Since k 31 remains unchanged, the increase in k 23 also accounts for the increased bunching level. On the other hand, an increase in k 41 is observed, which in contrast to k 23 has the effect of decreasing the bunching level. Evaluating the effect of the individual changes in the rates we conclude that the difference in bunching level between the dark and bright fluorescence levels is mainly due to the change in k 23 . One might argue that the difference in the background contribution for the dark and bright states, respectively, might cause changes of the bunching level. We found, however, that the dark state had always the same or a slightly higher background contribution than the bright state [for the measurements presented in Figs. 4(a) and (b), the respective values of σ in Eq. (2) are 0.95 and 0.96], which correspondingly decreases and increases the bunching level for the dark and bright state. Thus, without the background photocounts, the difference in bunching level would be even more pronounced.
The simultaneous acquisition of the fluorescence spectrum together with the fluorescence counts [Figs. 4(e)-(g)] proves a distinct correlation between fluorescence intermittency and spectral diffusion. This is expressed by a ≈ 1 nm red shift of the dark-level peak ZPL wavelength relative to the ZPL of the bright fluorescence level [ Fig. 4(g), two horizontal lines]. The timescale over which blinking is observed decreases with excitation power. Together with the spectral shifts this indicates that the process originates from site-specific and optically induced charge-state fluctuations, in agreement with similarly prepared samples on a SiO 2 substrate [14]. These fluctuations can be avoided by suitable surface passivation with for instance a thin layer of Al 2 O 3 [19].
We additionally recorded the time-resolved photon counts with pulsed laser excitation, in which case again from emitter A we observed telegraph-like blinking between a bright and a dark fluorescence level. The representative fluorescence time trace in Fig. 5(a) shows three distinct fluorescence levels, from which three separate photon-arrival-time histograms are obtained after post-processing the data [ Fig. 5(b)]. The lifetime of the respective fluorescence levels is well described by a time constant of 2.5 ± 0.1 ns apart from the brightest level, which contains an additional exponential component with a short time constant of 1.2 ± 0.1 ns. Considering the bandwidth of this approach, df = dτ/|τ| 2 ≈ 16 MHz, these results are qualitatively in agreement with the small changes of k 23 in the order of MHz obtained with CW excitation and recording of the second-order correlation function.

Internal quantum efficiency
The internal quantum efficiency (IQE), defined as the ratio of radiative decay rate to the total decay rate, is an important property of an emitter and characterizing its intrinsic brightness. It has been reported in numerous works that fluorescent emitters in hBN generally yield high photon count rates. Taking the lifetime into consideration, emitters hence should have a high IQE. Thus, a near-unity quantum efficiency was inferred from an emission rate of ∼ 4 × 10 6 counts/s at saturation [11]. In order to obtain a more quantitative estimate of the quantum efficiency, a calibration of the setup collection and detection efficiency has been made [5,15]. Then the IQE was determined by measuring the maximum count rate (at asymptotically infinite power), where unambiguously the count rate is proportional to the excited-state population. This yields values for the IQE of 0.65 for an emitter in hBN reported by Tran et al. [5], and 0.2-0.6 for one hBN emitter and 0.5-1 for another reported by Schell et al. [15]. Similar to our analysis, Tran et al. [5] fitted experimentally measured g (2) (τ) function at different excitation powers to extract radiative and non-radiative transition rates. From the reported rates we estimate the IQE for their emitter to be ∼ 97%.
With the pump-power-dependent relaxation rates obtained from the CW correlation-function measurements, we are also able to estimate the IQE of our emitters A-D. Assuming that the emitter dynamics are described by the level scheme shown in Fig. 2, the IQE can be calculated as η IQE = k rad /k tot = k 21 /(k 21 + k 23 + k 24 ) and is presented in Fig. 6. Since the shelving rates k 23 and k 24 are power dependent, the IQE also depends on excitation power. In agreement with the large saturation count rates observed, our estimation shows that IQE of all emitters is close to unity (above 95%) in the regime P P sat , but decreases with increasing excitation power. This power dependence explains a noticeable variation of the IQE values reported previously [5,11,15]. It has been argued for other single-photon sources (terrylene and dye molecules) [27,35] that power dependence of deshelving rates (in our case k 31 and k 41 ) results from optical excitation of higher states of the same multiplicity with a subsequent reverse intersystem crossing (relaxation) to the excited state. Fleury at al. [27] speculated that this process was responsible for remarkable photostability of terrylene at ambient conditions. This may as well be the case for fluorescent defects in hBN. We furthermore argue that shelving rates k 23 and k 24 depend on power due to thermally activated non-radiative processes. This has been experimentally demonstrated and studied for other color centers, e.g. nitrogen-vacancy centers in diamond [36]. It is important to note that the decrease of IQE with excitation power is an effect distinct from the saturation behaviour of a quantum emitter and is entirely due to non-radiative relaxation processes.

Conclusions
In conclusion, we report on the observation of bright single photon emitters in hBN nano-flakes at room temperature. For four emitters we carefully measure and evaulate the second-order correlation function in the regimes well below and above saturation power P sat . We observe photon anti-bunching at vanishing delay times proving single photon emission and photon bunching at long time scales up to the 100-µs range. The time scale and level of bunching decreases with power and indicates complex power-dependant emitter dynamics. Assuming a particular energy level scheme with corresponding transitions, we convincingly fit the experimentally measured second-order correlation functions and extract the transition rates of these emitters from the model. Our analysis shows that, in addition to optical ground and excited states, at least two additional levels contribute to the dynamics of these emitters. All emitters investigated here have similar optical spectra corresponding to type-I emitters in hBN reported previously in the literature.
Fluorescence intermittency with distinct fluorescence levels was observed for one emitter and the correlation function for each fluorescence level was individually characterized. Based on this analysis, we could associate the blinking behaviour mainly with changes in the internal quantum efficiency presumably caused by local and optically induced charge fluctuations. Under pulsed excitation at relatively high power, a similar blinking behaviour was observed, though no change in the total decay rate for the different brightness levels could be extracted from the lifetime measurements. Our observations indicate that the internal quantum efficiency of fluorescent defect centers in hBN is close to unity at low excitation power, but clearly decreases with power due to increased and power dependent non-radiative losses of the system. Fig. 7. Characterization of emitters B, C, and D, similar to analysis of emitter A presented in Fig. 3: Plots of transition rates obtained by fitting g (2) (τ) functions measured at different cw optical excitation powers with Eq. (2). Scaling of transition rates with excitation power is linearly approximated as follows: Emitter B: k 23 = 1.07(5) MHz P P sat , k 31 = 0.20(1) MHz P P sat , k 24 = 7.3(3) MHz P P sat , k 41 = 7.8(2) MHz P P sat ; Emitter C: k 23 = 7(2) MHz P P sat + 2.5(5) MHz, k 31 = 2.1(5) MHz P P sat , k 24 = 8.3(5) MHz P P sat , k 41 = 11.8(3) MHz P P sat ; Emitter D: k 23 = 0.4(1) MHz P P sat , k 31 = 0.5(1) MHz P P sat , k 24 = 3.0(3) MHz P P sat , k 41 = 5.2(5) MHz P P sat ;

A. Characterization of emitters B-D
Experimental data and fitting results for emitters B-D are presented in Fig. 7.

B. Saturation curves
The average count rate was extracted from the same time-tagged acquisition as were the g (2) (τ) functions. The data for the two detector channels were binned in 200-ms time slots for the full duration of the measurement. The data from two channels was added to obtain the combined count rate, from which the average count rate and standard deviation were calculated. The average count rate as a function of excitation power for all emitters is presented in Fig. 1(c) with error bars indicating one standard deviation. The curves are fitted with I = I ∞ P P+P sat + bP + c, where a lower limit of 1 kHz was imposed on the value c, which accounts for the combined dark count rate of both detectors. The resulting fit parameters are presented in Table 1.

C. Estimation of errors for transition rates k ij
In order to estimate the fitting errors for the rates k ij , we first calculate the Jacobian matrix J mn = ∂r m ∂x n , where r m is the vector of residuals and x n are the fitting parameters k ij . We then calculate the covariance matrix as K = σ 2 (J T J) −1 , with σ being the standard deviation of the residuals of the fit procedure. Each diagonal element of this matrix corresponds to the variance of the fitting parameter x n . The standard error of the rates is then found as the square root of the diagonal elements of K.