{'fraction_non_alphanumeric': 0.0363003003003003, 'fraction_numerical': 0.02142942942942943, 'mean_word_length': 4.170289405042852, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photophysics of single nitrogen-vacancy centers in diamond nanocrystals

A study of the photophysical properties of nitrogen-vacancy (NV) color centers in diamond nanocrystals of size of 50~nm or below is carried out by means of second-order time-intensity photon correlation and cross-correlation measurements as a function of the excitation power for both pure charge states, neutral and negatively charged, as well as for the photochromic state, where the center switches between both states at any power. A dedicated three-level model implying a shelving level is developed to extract the relevant photophysical parameters coupling all three levels. Our analysis confirms the very existence of the shelving level for the neutral NV center. It is found that it plays a negligible role on the photophysics of this center, whereas it is responsible for an increasing photon bunching behavior of the negative NV center with increasing power. From the photophysical parameters, we infer a quantum efficiency for both centers, showing that it remains close to unity for the neutral center over the entire power range, whereas it drops with increasing power from near unity to approximately 0.5 for the negative center. The photophysics of the photochromic center reveals a rich phenomenology that is to a large extent dominated by that of the negative state, in agreement with the excess charge release of the negative center being much slower than the photon emission process.


I. INTRODUCTION
With the development of photonic quantum cryptography and quantum information processes, there is need for reliable and easy-to-use single photon sources. Such sources have been developed in recent years 1 , like single molecules, colloidal or epitaxial semiconductor quantum dots, and color-centers in diamond. [2][3][4] Here, we are interested in the latter, namely the NV center, formed by a substitutional nitrogen atom adjacent to a vacancy in the diamond lattice. NV centers have found numerous applications recently thanks to their unique physical properties such as excellent photostability 2,5-7 and long spin coherence times 8 as well as to improved control over both their production 9 and physical initialization protocols. 10 NV centers can be made available both in ultra-pure bulk diamond 8 and ultra-small crystals. 11 Applications range from high-sensitivity high-resolution magnetometry [12][13][14][15][16][17][18] , to fluorescence probing of biological processes 19,20 , solid-state quantum information processing 21,22 , spin optomechanics 23,24 , quantum optics [25][26][27] , nanophotonics [28][29][30][31] , and quantum plasmonics. [32][33][34] NV centers can take two different charge states with different spectral properties : the neutral center NV 0 , which has a zero-phonon line (ZPL) around 575 nm (2.16 eV), and the negatively charged center NV -, which has a ZPL around 637 nm (1.95 eV). 35 In addition to ZPLs, the fluorescence spectra of both centers exhibit a broad and intense vibronic band at lower energy. A single NV center, as a single-photon emitter, is characterized by a secondorder time-intensity correlation function that exhibits a photon antibunching dip at zero delay. In a first ap-proach, the photophysics of NV centers can be modeled by a two-level system with two photophysical parameters, the excitation rate, and the spontaneous emission rate. However, with increasing excitation power, the NV center, more particularly in the negative state NV -, can experience distinctive photon bunching at finite coincidence time, in addition to the expected antibunching at zero delay. This can be accounted for within a three-level system with additional photophysical parameters to describe photon decays to, or from, the additional shelving level. The aim of this paper is to give a detailed description of the intrinsic photophysics of single NV centers of both charge states in surface-purified 9 nanodiamonds (NDs) of size around 50 nm, or below, as a function of the illumination power from a CW laser. Understanding the intrinsic photophysics of NV centers is required before implementing small fluorescent NDs in a complex electromagnetic environment, such as practical single-photon devices, which will modify the photophysics. 28 It is also useful to the applications mentioned above as most of them exploit the single-photon emitter nature of isolated NVs. The statistics of both NV charge states has been studied previously in ND samples similar to those studied here and it was found to be size dependent, with a larger occurrence (> 80%) of the NVcenter over the entire size range from 20 to 80 nm. 9 In the present study, we further observed that most of the single NV centers being detected as neutral from their fluorescence spectrum at low illumination power (<0.5 mW) progressively gain with increasing power a photochromic character in that they also exhibit the NV -ZPL in addition to the NV 0 ZPL. Only a few of NV 0 s remain purely neutral over the entire power range, which we call pure NV 0 behavior. On the other hand, NV centers detected as negatively charged at low power remain so with increasing power, which we call pure NVbehavior. We first focus our attention on such pure NV 0 or NVcenters. In addition to such behaviors, we found that some rare NVs can see their charge switching between neutral and negative 36 already at the lowest excitation power. We also describe the photophysics of such a photochromic NV over the same power range as the non-photochromic centers and show how both photophysics can be linked. This allows us to find valuable information on the dynamics of photochromism. The paper is organized as follows. The experimental methods are described in Section II. The three-level model used to interpret the experiment is developed in detail in Section III. Section IV focuses on the experimental results and on the extraction of the various photophysical parameters as a function of the excitation power for a NV 0 and a NVcenter. Section V describes the photophysics of photochromism in a single NV. A summary is given in Section VI.

II. EXPERIMENTAL METHODS
Preparation of the ND sample was achieved following a procedure reported previously. 9,37 Commercial HPHT diamond nanocrystals are first irradiated using highenergy electrons, then annealed at 800˚C in vacuum to produce the fluorescent NV centers, and finally annealed in air at 550˚C to remove surface graphitic compounds. Colloidal dispersion in water and further sonication allow us to obtain a uniform solution of NDs. We consider here NDs with a typical size of 25 nm or 50 nm deposited on a fused silica substrate to minimize spurious fluorescence. 29 Single NV centers were optically addressed using standard confocal microscopy at room temperature. 29 A CW laser light (wavelength λ exc = 515 nm) falling within the absorption band of the NVs is used to excite the NV fluorescence. Excitation light was focused onto the sample using an oil immersion microscope objective with numerical aperture NA= 1.4. The NV fluorescence is collected through the same objective and is filtered from the remaining excitation, i.e., with wavelength below 532 nm, by a dichroic mirror and a high-pass filter. The collected fluorescence is subsequently sent either to a Hanbury-Brown and Twiss (HBT) intensity correlator (see below) or to a spectrometer. An example of fluorescence spectrum is shown in Fig. 1(b) for the NVcase. The ZPL corresponds to the resonant decay ( Fig. 1(a)) at λ = 637 nm while the wide fluorescence side-band is the phonon replica.
The second-order time-intensity correlation function contains the information on the classical versus quantum nature of light. It reads in the stationary regime g (2) (τ ) = P 2 (t + τ |t)/P 1 (t) (1) Figure 1: (Color online) (a) Schematic of the two-level system explaining the emission spectrum of the NV center. g : ground state, e : excited state. The several thinner horizontal lines represent the coupling with the phonon bath in the diamond matrix. The green arrow represents the excitation and the red ones the optical decay at resonance (full arrow) and out of resonance (dotted arrow) (b) Typical emission spectrum of an NVcenter in a 25 nm ND.
where P 2 (t + τ |t) = P 2 (τ |0) is the conditional probability to detect a photon at time t + τ knowing that another photon has been recorded at time t. This probability is normalized by the constant single photon detection rate P 1 (t) = P 1 (0). For a classical source of light g (2) (τ ) ≥ 1, whereas the observation of an anti-bunching g (2) (τ ) ≤ 1 is a clear signature of the quantum nature of light. [38][39][40] In particular, at zero delay g (2) (0) = 0 for a single photon emitter 40 , which means that the probability to detect simultaneously two photons vanishes. In practice, an HBT correlator ( Fig. 2(a)) 5,6 is used to measure g (2) . Here, the fluorescence of the NV center is sent on a beam splitter, which separates the signal in two equal parts sent to two avalanche photo-diodes (APDs) named APD1 and APD2. The APDs are connected to a time-correlated single-photon counting module to build histograms of delays between photon events detected by the upper « Start » APD1 and the lower « Stop » APD2. In order to avoid unwanted optical crosstalk between the APDs, a glass filter acting as a short-pass filter at 750 nm and a diaphragm are added in both branches. 41 In the standard configuration, no bandpass filter is added to the setup, in contrast with the configurations used to study photochromism (Section V), so that it can be used for both charge states of the NV (provided that the NVs do not experience charge conversion, which is the case for the selected NVs in the present work).  Fig. 2(b) shows a typical g (2) function measured for a single NV at low excitation power. The experimental data are compared with an equation stemming from a two-level model 5,29 : where r and γ are the excitation and spontaneous emission rates, respectively. Within this model g (2) (τ ) < 1, which means that the emitted light is non-classical at any delay τ . However, at higher excitation rate this simple model (Eq. 2) generally fails. This is particularly true for NVcenters. In this case, as shown in Fig. 2(c), the g (2) function includes a bunching g (2) (τ ) > 1 feature at finite delays, superimposed to the antibunching curve. This kind of correlation profile, which contradicts Eq. 2, calls for a third level that traps the electron, preventing subsequent emission of a photon for a certain time. [42][43][44][45][46][47][48][49] Therefore, for these delays, the correlation function is higher than one. Though well known for the NVcenter, the existence of a shelving level is less documented for its neutral counterpart but has been invoked theoretically 48 as well as in electron-paramagnetic resonance. 49 From now on, we will describe phenomenologically the color center dynamics by a three-level model and deduce the intrinsic photo-physical parameters of pure NVand NV 0 centers (Section IV) and of a photochromic center (section V). These NV centers have been selected from their fluorescence spectra.

A. Correlation function of a three-level system
In order to give a quantitative description of the threelevel model we first review briefly the interpretation of the g (2) function using Einstein's rate equations. Starting from Glauber quantum measurement theory 50 we have where I(t) is the quantum operator equivalent to the electromagnetic energy flow absorbed by a detector at time t and :: represents normal ordering. 38,39 Introducing the creation and annihilation photon operators a † (t) and a(t) in the Heisenberg representation, we obtain : For a two-level system with ground state g and excited state e, the creation operator a † (t) at time t is to a good approximation proportional to the rising transition operator |e, t g, t| 51 , we deduce where p(e, t + τ |g, t) is the conditional probability for the NV to be in the excited state e at time t + τ knowing that it was in the ground state g at time t. Like in Eq. 3 this probability is normalized to a single event probability, i.e., the probability p(e, t) for the NV to be in the excited state at the previous time t. The main interest of these equations is to link the probability of detection, as given by Eqs. 1 and 3, to the emission probabilities p(e, t + τ |g, t) and p(e, t) defined by the rate equations. Therefore, g (2) can be completely determined if the transient dynamics of the emitter is known.
In the context of the NV center, we must use a two-level system with a third metastable state to explain the bunching observed in the correlation measurements. 52 Since the previous calculations only considered the excited and ground states involved in the fluorescence process we will admit (see 38,54 for a discussion and justification) that these results still hold with a three-level system if we replace the excited and ground states by the levels 1 and 2 in the Jablonski diagram, respectively (see Fig. 3). Here we neglect the channels 1 to 3 and 3 to 2, because the system is not supposed to be excited at these transition energies, in contrast with previous models. 26,55 There, channel 3 to 2 was taken into account and channel 3 to 1 was neglected, because the quantum yield associated with the NV relaxation was supposed to be close to unity, while recent studies show that Q ≃ 0.6 − 0.7. 2,28,56-58 Therefore, we here take into account four channels. This leads to four unknown parameters : the excitation rate, k 12 , the spontaneous emission, k 21 , and the two parameters k 23 and k 31 of the additional decay paths involving the shelving level. k 21 being the only radiative channel, we can write the autocorrelation function as : where p 2 (t) is the population of state 2 at time t. p 2 (∞) represents the asymptotic limit of p 2 (t) when the transitory dynamics approaches the stationary regime. p 2 (t) can be obtained by solving the system of rate equations defining the three-level system presented in Fig. 3. With , the population of state i we can write the following set of equations : whereṗ i means time derivative of p i (t). The use of rate equations instead of Bloch equations is here fully justified since in ambient conditions, the coherence between levels is decaying very fast. 55 Eqs. 7 show that the system is necessarily in one of the states at any time. The steady-state analysis of these equations permits to find the explicit definition of the fluorescence rate R at which the system emits photons. 54 This definition allows us to explain the saturation behavior of the NV fluorescence, i.e., R tends towards a finite value for increasing excitation power. However, we are here looking for the timedependent analysis to find p 2 (τ ) as a function of the k ij coefficients. If we eliminate p 3 from Eq. 7, we obtain : The resolution of this pair of equations, with the initial conditions p 1 (0) = 1 and p 2 (0) = 0, leads to p 2 (τ ) = p(2, τ |1, 0) and therefore to the expression of g (2) (τ ). Important approximations can be made to obtain a simple expression. More precisely, in addition to neglecting channels 1 to 3 and 3 to 2, we suppose that : This is justified since, even if the third level is considered, the associated rates are supposedly very small compared to the singlet rates k 12 and k 21 . We will see latter that this is not really true for the NVcenter but that the results obtained are actually robust and keep their meaning even outside of their a priori validity range. Within the above mentioned approximations the second-order correlation function reads (see Appendix A for mathematical details) : where the parameters γ 1 , γ 2 and β are defined through the relations :

B. Determination of the kij coefficients
The aim of this sub-section is to determine the k ij coefficients for the specific measured NV centers. The fit to the g (2) (τ ) function allows us to determine γ 1 , γ 2 and β. From Eqs. 11-13 we deduce : We then have three equations for four unknown variables. A fourth equation is needed to solve the problem entirely.
In the experiment, we also have access to the radiation or fluorescence rate R, measured in s −1 . This rate is simply the average number of photons that the APDs collect per second. It represents, up to a multiplicative coefficient associated with the photon propagation in the setup, the probability for the system to be in level 2, multiplied by the transition probability k 21 to relax (supposedly by optical means) to the ground state. We have : where ξ is the collection efficiency of the system once the NV center has emitted its fluorescence. Note that the same formula was used in refs. 26,55 with a different definition of p 2 (+∞) because of a different Jablonski diagram. However, this has a physical meaning only if k 21 is associated with a pure radiative decay. Whereas the model of ref. 26,55 implied a unity quantum yield, the quantum yield in our approach is defined by In our phenomenological approach the third level is thus assumed to absorb all of the non-radiative transitions letting k 21 be a pure radiative decay. Finally we remind that the probability p 2 (+∞) needed in Eqs. 17 and 6 should be calculated in the asymptotic stationary regime, which can be obtained by canceling allṗ i in Eq.8 : .
Now four equations are at hand so that the system can be inverted to determine all coefficients. This calculation involves the numerical resolution through Cardano's algorithm 59 of a third-order polynomial as given in Appendix B.

IV. EXPERIMENTAL RESULTS AND EXTRACTION OF THE kij PARAMETERS
We here consider two representative examples of NDs, one hosting a single pure NVcenter, and the second one hosting a single pure NV 0 center. The NV --center ND is about 25 nm in diameter, whereas the one hosting the NV 0 center is about 50 nm in diameter. 60 The g (2) function was recorded for both NDs with different excitation powers P exc in order to study the evolution of the k ij coefficients. Fig. 4 depicts the g (2) function for the two NV centers and for excitation powers P exc ranging from 200 µW to 10 mW . The experimental curves are fitted with Eq. 10 taking into account the correction for the incoherent background light collected by the APDs. 5,6,29 The experimental antibunching dip does not drop to zero due to this incoherent background, which modifies Eq. 10 as exp. (τ ) = g (2) (τ )ρ 2 + 1 − ρ 2 , where ρ = S/(S + B) contains the signal S and background B contributions from the NV fluorescence and the spurious incoherent light, respectively. By recording the average intensity from the APD directly on the NV and at a location close to it, we experimentally determine ρ and the fit parameters of Fig. 4 as explained in refs. 5,6,29 . It is seen that the increase of P exc only induces a very small bunching for the NV 0 , contrary to what is observed for the NV -. However, the anti-bunching dip narrows with increasing power in both cases. Our observations are collected in Fig. 5, which shows the power evolution of the fit parameters γ 1 , γ 2 and β. The general trends seen in Fig. 4 are confirmed in Fig. 5 for both NVs since γ 1 , which is associated with the antibunching contribution, is clearly increasing with P exc , a fact which is reminiscent from Eq. 2. Furthermore, it is seen that γ 2 is also increasing significantly for the NVwhich is clear signature of the third energy level. As far as the β parameter is concerned, it remains at a constant value β ≃ 1 for the NV 0 , while it increases up to β ≃ 7 for the NV -. This behavior agrees with its definition from Eq. 10. Therefore, when β is very close to 1, there is no significant bunching, and when β increases with power, the bunching turns on. Now that the three parameters have been found, the k ij parameters can be traced back. This will be done in two steps.
A. kij parameters with constant k21 In Eq. 17 the collection efficiency ξ of the optical setup must be known precisely to extract the various k ij . However, since this can only be estimated, we will first calculate the k ij parameters by assuming that k 21 does not change with the excitation power. This hypothesis is intuitive because the parameter k 21 , i.e., the spontaneous emission rate, is supposed to be solely governed by the Fermi's golden rule, which in turn depends on the electromagnetic environment only. In order to determine the value of k 21 we observe that according to Eq. 14 we must have γ 1 = k 21 at zero excitation since in this case k 12 = 0. From the linear regression for γ 1 (Fig. 5) we deduce γ Thanks to Eqs. 14 to 16, we deduce the three other parameters as shown in Fig. 6 for both NVs. The first point to notice is that for both NVs, the k 12 parameter increases linearly with the pumping rate from zero to a value exceeding k 21 for P exc ≃ 2 mW . Moreover, for NV 0 , we see that k 31 increases, but keeps very small values compared to the other parameters, in agreement with the assumptions made in Eq. 9. However, for NV -, the same parameters are no longer negligible compared to the set of k 21 , k 12 values. They even overtake k 21 for P exc > 6 mW . However, we emphasize that assuming k 21 constant, if natural, is actually not fully demonstrated. In order to check how robust this hypothesis is, we will now try to approach the values and evolutions of the k ij parameters that were obtained here by modulating the collection efficiency ξ using Eq. 17.

B. Modulation of the collection efficiency ξ
Now we use the fourth equation Eq. 17 to calculate the k ij parameters. As already stated, we do not know exactly the ξ parameter but, as it turns out, slight variations in ξ can produce significant changes in k ij . To find the correct value of ξ, we adopt the following procedure. We let its value vary continuously and calculate the evolutions of the slope and Y-intercept of the linear regressions made with the obtained k 21 (P exc ) traces. The slope should vanish because k 21 is assumed to be constant, whereas the Y-intercept should reach the value obtained previously, i.e., the value of γ 1 at zero excitation. Therefore, we calculate the evolution of the k ij parameters as a function of the collection efficiency ξ. The results are shown in Fig. 7, where the constant horizontal curves depict the values to be reached by the slope and Y-intercept. For the NV 0 ( Fig. 7 (a)), there is indeed a value of ξ where the two parameters reach the assumed values (blue vertical line). This gives ξ (0) = 0.77×10 −3 , in agreement with a rough estimate of our setup collection efficiency taking into account the various optical compo-nents. However, for the NV - (Fig. 7 (b)), it is seen that the parameters do not reach exactly the previous values. Yet, there is an optimum ξ for which the k ij approach the previous values (blue vertical line). This corresponds to ξ (−) = 1.6×10 −3 , which differs only by a factor 2 from the NV 0 case. This appears reasonable since the measurements were not carried out the same day (optical alignments might be slightly different) and the collection efficiency ξ also depends on the unknown transition-dipole orientation in both NVs. Fig. 8 depicts the k ij (P exc ) curves deduced from the previous optimization. It is found that the order of magnitude of the coefficients is the same as for imposed k 21 . In particular, the relaxation rate k 31 is unchanged because it only depends on the g (2) fit parameters, Eq. 15. For the two NV centers, the excitation rate vanishes in the absence of any excitation power and linearly increases with P exc as it should. The main difference between both centers comes from the evolution of k 23 , k 31 and k 21 . Indeed, for NV 0 , k 23 and k 31 remain very small compared to k 21 , which is almost constant (1/k 31 ≈ 500 ns and 1/k 23 ≈ 1000 ns if P exc tends to zero). Therefore, the third level plays little role in the photodynamics of the NV 0 center. However, it is worth stressing that although very small, k 23 is not zero for NV 0 (errors bars are within the symbol size in Fig. 8 (a), right panel), which confirms the very existence of the third metastable level for this charge state. Regarding the g (2) curves, it implies that the optical channel 2 to 1 is favored, which prevents any significant bunching. Furthermore, for the NV 0 the narrowing of the antibunching dip is due to the increasing excitation rate as for a two-level system, i.e., Eq. 2 (see ref. 29 ). The analysis of the k ij parameters for the NVis more involved. Indeed, k 23 increases very quickly to reach the order of magnitude of k 21 , while the latter is increasing as well (at zero excitation power we have 1/k 21 ≃ 24 ns and 1/k 23 ≃ 500 ns). The increase of k 23 , associated with non-radiative transitions, actually explains the growth of the bunching feature on the g (2) curves. Although the physical justification of this finding is beyond our phenomenological treatment, it is likely that the variation of k 21 and k 23 with excitation power is due to a change in the local energy environment of the NV center at high power, in particular because the efficient coupling with the phonon bath in the diamond matrix is expected to be temperature sensitive. Increasing the excitation power could thus correspond to an increasing effective temperature, subsequently affecting the relaxation dynamics. It is worth pointing out that there is a limitation in the analysis done for the NVcenter. Indeed, the very fact that k 23 and {k 21 , k 12 } reach the same order of magnitude contradicts the hypothesis made in Eqs. 9-13. Actually, as already mentioned, the results obtained are much more robust that could be anticipated at first sight. This can be figured out by relaxing the constraint of Eq. 9 as done in the detailed calculation presented in Appendix C. The results obtained with the new rate-equation model are shown in Fig. 9 for the evolution of the k ij coefficients To complete the analysis we also computed the quantum yield evolution as given by Eq. 18 and compared with the evolution of the fluorescence rate. Fig. 10 confirms that within the Jablonski model sketched in Fig. 3 and in the considered excitation regime, the quantum yield of the NV 0 center is approximately constant Q ≃ 1, in agreement with the intuitive fact that the third level does not play a significant role in the dynamics. In contrast, the NVquantum yield decreases dramatically with increasing excitation power from a starting Q ≃ 1 to Q ≃ 0.5 at high power. This entails the fact that the NVdynamics is strongly dependent on the excitation power as discussed before.

V. PHOTOCHROMISM
Photochromism of NV centers has been reported in ensembles of NV centers in CVD diamond films under additional selective illumination 61 , with single NV centers in 90 nm NDs under femtosecond illumination, which results in the photo-ionization of the negative center to its neutral counterpart 35 , with ensembles of NV centers in type Ib bulk diamond at cryogenic temperatures under intense CW excitation 62 , and with a single center in natural type-IIa bulk diamond under CW illumination. 36 In this last report, a special scheme of cross-correlation photon measurements was applied in the emission band of both charge states to show that the collected fluorescence in the NV 0 and NVstates were correlated and originated form a single NV defect. Several studies have reported that charge conversion within the NVs critically depends on the illumination conditions. 57,63,64 A complete understanding of NV photochromism is lacking but a widespread view is that the optical excitation, either CW or transient, tunes the quasi Fermi level around a NV charge transition level, thereby inducing charge conversion. 61,65 This scenario has been reinforced recently by electrical manipulation of the charge state of NV ensembles and of single NVs by an electrolytic gate electrode used to tune the Fermi energy. 65 Depending on the sort of diamond studied, photochromism is thought to be favored by the presence of electron donor or acceptor defects, such as nitrogen, in the neighborhood of the NV center. 66 Recently, it was also found that resonant excitation of the NV 0 and NVstates in ultrapure synthetic IIa bulk diamond can induce reversible charge conversion in cryogenic conditions even at low power. 10 This was taken as evidence that the charge conversion process is intrinsic in this sort of diamond, not assisted by an electron donor or acceptor state. The goal of this section is to give additional information on NV photochromism detected in surface-purified NDs, 25 nm in size, subjected to a CW non-resonant excitation of increasing power. By comparing the behavior of a single photochromic center to that of non-photochromic centers in the same illumination conditions as described above, we gain valuable information on the rich photophysics of photochromism. For the purpose of studying photochromism, we use two additional configurations of the HBT correlator that differ only by the set of bandpass filters added in the interferometer branches. These configurations, called NV -/0 , and NV 0/respectively, are shown schematically in Fig. 11. In contrast to Fig. 2(a), these two configurations add selective bandpass filters in the interferometer branches. The NV -/0 configuration in   Fig. 11 (b)) uses a filter selective to the NV -(respectively NV 0 ) fluorescence in the start (respectively stop) branch. Therefore, in the NV -/0 configuration, a single photon emitted by a NVcenter and detected in APD1 gives the « Start » signal to the counting module, whereas a single NV 0 photon subsequently detected in APD2 produces the « Stop » signal. The NV 0/configuration works in just the complementary way. Note that our NV -/0 configuration is similar to the cross-correlation technique used in ref. 36 . These schemes turn out to be very powerful to study photochromism since cross correlation can be expected only if NVand NV 0 photons originate from the same defect center. Note that related techniques have also successfully been applied to identify various excitonic species emitted by single semiconductor quantum dots, see for instance. 67-71 With the above setup, we have located a few rare NDs hosting a single NV center that showed charge conversion at low excitation power already. In the following, we consider such a ND hosting a single photochromic NV   exp. function. The fit parameters γ1,γ2 and β using a 3 energy-level model are also given.

center.
The relevant spectra are shown in Fig. 12(a). It is found that both the NV 0 and NV -ZPLs are seen at any excitation power. The corresponding g (2) function measured in the standard configuration of the HBT (no bandpass filter) is shown in Fig. 12(b). It reveals a clear antibunching dip at zero delay. A precise analysis taking into account the background light (see Table 1) confirms this finding since at high excitation power we observe that the background B increases significantly with respect with the NV fluorescence signal : ρ = S/(S + B) ≃ 0.6 while at low power ρ ≃ 0.9 − 1. With this value and using the formula g (2) exp. (τ )/ρ 2 + 1 − 1/ρ 2 = g (2) (τ ) we deduce the actual value of g (2) (0) ≃ 0. 6,25,55,72 From these results, a natural interpretation for the observation of the NV 0 and NV -ZPLs together with NV uniqueness is that this particular ND is subjected to photochromism. In addition to the antibunching dip at zero delay, it is seen in Fig. 12 that g (2) exceeds 1 at longer delays. 55 In agreement with the previous sections, this is evidence for the presence of a trapping level and calls for a three-level description of the photochromic NV. The values of γ 1 , γ 2 and β parameters used for fitting the g exp. function agree qualitatively well with those obtained for the NVconsidered previously in agreement with the fact that the system is acting more like a three-level system. Although our study of NV photochromism is limited to a particular example we suggest that the dynamics of the system involves probably all energy levels of the NVand NV 0 with some possible hybridization. It could be interesting to know, whether or not, the third level involved in the photochromic case is identical in nature to the third level of the NV -. The role of the environment or of the radiation power 73 on the dynamics could be investigated in the future to clarify this point. Switching now to the cross-correlation regime we show in Fig. 13 the result obtained for P exc. = 3 mW (similar features, not shown here, were observed at other excitation powers). It is worth emphasizing that cross correlations, as described in e.g. refs. 68,69,74 , allow us to characterize the transitory dynamics between the two NV configurations. Using a formalism equivalent to the one leading to Eqs. 3-5 we indeed obtain in the NV -/0 configuration sketched in Fig. 11 (a) where p(e, NV 0 , t + τ |g, NV -, t) is the conditional probability for the NV to be in the excited level e of the NV 0 state at time t + τ knowing that it was in the ground energy level g of the NVcharged state at the previous time t. This also means that a first photon emitted by the NV 0 was detected at time t while a second photon emitted by the NVis detected at time t + τ . In a symmetrical way using the NV 0/configuration sketched in Fig. 11 (b) we get with similar definitions as previously but with the role of NV 0 and NVinverted. The experimental results corresponding to these two configurations are shown in Figs. 13(a) and (b), respectively. Here, it is also important to have τ ≥ 0 in the calculation leading to Eqs. 20 and 21 in order to have a clear physical understanding. However, the electronic delay ∆ = 100 ns included in the HBT correlator setup implies that sometimes even a photon emitted, say at time t 0 , by a NVis recorded after a second photon emitted later by the NV 0 state, i.e., at t 0 +τ . This corresponds to the 'negative' delay part of the graph, i.e., Fig. 13(a), which is actually associated with the inverse dynamics NV 0/-, i.e., Eq. 21. Here for clarity we did not subtract the delay from the abscises in Fig. 13. Therefore for τ > ∆ in Fig. 13(a) we have g (2) exp. (τ ) = g (2) −/0 (τ −∆)ρ 2 +1−ρ 2 while we have g 0/− (∆ − τ )ρ 2 + 1 − ρ 2 for 0 < τ < ∆ with ρ = S/(S + B) as previously. In the same way for Fig. 13(b) we have g −/0 (∆ − τ )ρ 2 + 1 − ρ 2 for 0 < τ < ∆. As it is clear from the definitions these cross-correlations should present some symmetries. In the present case we used a three energy-level fit, i.e. Eq. 10, for the theoretical functions, Eqs. 20 and 21. The parameters obtained to reproduce the data are the same for Fig. 13 (a) and (b) up to an inversion between the 'positive' and 'negative' delay for each graph. Therefore, we get for Fig. 13(a) γ 1 = 0.046 ns −1 , γ 2 = 0.01 ns −1 and β = 1.2 for τ < ∆, while we have γ 1 = 0.052 ns −1 , γ 2 = 0.01 ns −1 and β = 1.5 for τ > ∆. For the second cross-correlation curve the parameters are identical but the roles of τ > ∆ and τ < ∆ are inverted as it should be. We observe that these values are very close to each other and also from the one obtained in Table 1 at the same excitation power P exc. = 3 mW . This confirms that the system acts here mainly as a NVcenter. Interestingly, the last finding implies that after the emission of a photon in the spectral fluorescence band of the NVthe delayed emission of a second photon in the spectral band of the NV 0 , i.e., the conditional probability given by Eq. 20, is also characterized by the dynamics of the NVcontrary to the intuition. Such behavior was reported in ref. 36 for NV centers in bulk. In particular, in this paper the time dependence of the conversion NV 0 to NVprocess (and its inverse) was studied using pulse sequences. It was found that the relaxation from NVto NV 0 is a very slow process occurring with a decay time ≃ 1 µs. This agrees with our finding in Fig. 13 since, even if the full dynamics of the photochromic NV center is expected to depend on the energy levels of both the NVand NV 0 centers, it is clearly the NVcharacteristics which dominate during the transition associated with Eq. 20. A similar qualitative analysis can be done in the NV 0 to NVcenter conversion. The small dissymmetry between the two parts of the curves for positive and negative delays results from the presence of a small NV 0 contribution to the dynamics during the NV 0 to NVtransition which is absent in the NVto NV 0 conversion. Clearly, this complex charged/uncharged transition dynamics would deserve systematic studies in the future.

VI. SUMMARY
To summarize, we have experimentally studied the fluorescence photodynamics of NVand NV 0 centers in diamond nanocrystals of 50 nm size or below using HBT photon-correlation measurements as a function of the excitation power. The dynamics was theoretically modeled using Einstein's rate equations and the transition probability rates k ij entering the three-level model developed to analyze the data were deduced and used to infer a quantum efficiency to both charge states of the NV. It has been found that the shelving state, though present, plays a very small role on the neutral center in those small diamond crystals. The narrowing of the antibunching dip observed with increasing power for this center is a simple power effect that does not affect the near unity quantum efficiency. In contrast, the negative center experiences a distinctive photon bunching behavior at finite delay that increases with increasing power. This reflects the increasing role of the shelving state for this center, which in turn diminishes the quantum efficiency from near unity at low power to approximately 0.3 at high power. We have also studied the dynamics of a photochromic center. It reveals a rich phenomenology that is essentially dominated by the negative face of the center. In the future, it would be interesting to determine how the presence of, e.g., a plasmonic structure affects the dynamics of NV centers 28 in any charge state, in particular if the role of the shelving state in the NV 0 center as well as the neutral face of the photochromic center can be modified.

Annexe A
Here we briefly summarize the derivation of the main equations used in sections III and IV. We start from Eq. 8 written in a matrix form aṡ where For the sake of generality we here keep all transition coefficients k ij allowed by the three-level model. Such an Eq. A1 can be formally solved by defining a new vector q related to P by the equation P = e Mt · q. This leads to the new equationq = e −Mt J which has the general solution : where M −1 is the inverse matrix of M and where the initial condition corresponding to a system in the ground state at time t = 0 is q(0) = P(0) = 1 0 . Therefore, we deduce : In order to give an explicit form to the solution we need to diagonalize the M matrix. Therefore, we write The eigenvalues λ 1 and λ 2 of M are easily obtained from the secular equation det[M − λI] = 0. They reads : From Eq. A4 we easily obtain with In order to determine completely the solution we need to specify the transformation matrix Π = u 1 u 2 v 1 v 2 whose column vectors X i = u i v i with i = 1, 2 are solutions of the eigenvalue problem M · X i = λ i X i . These eingenvectors are determined up to an arbitrary normalization and here we choose Using Eq. A7 we therefore get : We can further simplify the solution since from Eq. A9 and the definition of M −1 we easily deduce After regrouping all terms we obtain p 2 (t) p 2 (∞) = − v 1 α v 1 α + v 2 ε e λ1t + ( v 1 α v 1 α + v 2 ε − 1)e λ2t + 1 which up to notations is equivalent to Eq. 10 of Section III if we write λ i = −γ i and β = v1α v1α+v2ε .

Annexe B
The results obtained in Appendix A are exact and no approximation were made for calculating the coefficients γ i and β. Now we will use the fact that {k 21 , k 12 } ≫ {k 23 , k 31 } to simplify and explicit these coefficients. First, we point out that we have which is Eq. 13. In order to solve the system of equations 11-13 we first eliminate k 23 from Eq. 12 using Eq. 13, i.e.
Inserting this result into Eq. 11 leads to which constitutes our first parameter solution. In order to determine the other parameters we insert the value obtained for k 31 into Eq. B5 and Eq. 11 to obtain k 23 = γ 1 γ 2 (β − 1) k 12 β k 21 = γ 1 − k 12 . (B7) In other words, all parameters are now expressed as a function of the excitation coefficient k 12 . To complete the solution we use Eqs. 17 and 19 for the single-photon rate R. After insertion of Eqs. B5 and B7 into Eq. 17 and some lengthly rearrangements we finally obtain k 3 12 − γ 1 k 2 12 + R ξ (γ 1 + X k 31 )k 12 + RX ξ = 0 (B8) with X = (β − 1)γ 1 γ 2 /β. The three roots of this cubic equation can be obtained numerically using Cardano's method.

Annexe C
We can generalize the analysis made in the previous Appendix B without any approximation. For this we first remark that we have the exact relations To complete the resolution of the system we use the exact relation Eq. B1. We finally obtain k 31 = γ 2 γ 1 /[γ 1 − γ 2 )β + γ 2 ], k 12 = B/2 − [B 2 − 4(C + A − Bk 31 )]/2, k 21 = C/k 12 , k 23 = B − k 12 − k 21 .

(C4)
The minus sign was chosen in the second equation for k 12 (i.e. solution of Eq. C3) since a Taylor expansion at low excitation power when γ 2 ≪ γ 1 , β ≃ 1 gives for the two roots In order to have the linear regime k 12 ∝ R/ξ we must therefore impose the minus sign.