Multidimensional spectroscopy and imaging of defects in synthetic diamond: excitation-emission-lifetime luminescence measurements with multiexponential fitting and phasor analysis

We report the application of phasor analysis and nonlinear iterative fitting to complex spatial and spectroscopic luminescence decay data obtained from multidimensional microscopy of a CVD diamond grown on a HPHT substrate. This spectral and lifetime-resolved analysis enabled spatial mapping of variations in concentrations of nitrogen vacancy (NV) defects in both charge states and the quenching of NV− defects, as well as the identification of SiV− luminescence. These imaging and spectroscopic modalities may be important for reliable fabrication of quantum devices based on such defects in diamond, which will require well-defined and characterised quantum electronic properties.


Introduction
Diamond is a mechanically strong and optically transparent host for a variety of complex luminescent defects, such as the nitrogen vacancy (NV) [1] and the silicon vacancy (SiV) [2]. These 'colour centres' have been extensively studied, because of their potential for quantum sensing, communication and computing applications. Manufacturing pure diamond, or diamond with well-controlled defects, is nontrivial as unwanted impurities can be introduced during growth. It is therefore important to understand the interactions between defects to be able to fabricate diamond-based quantum devices.
Diamond colour centre emission is complex, arising from a range of emitting defects [3]. Through quenching mechanisms of these defects, the quantum electronic properties such as luminescence lifetime can be impacted [4][5][6][7]. Spectrallyresolved imaging can be used to distinguish and map emitting defects, but at room temperature, many defects spectrally overlap due to broad phonon side bands [1,5]. This makes separation non-trivial and can obscure weak emission. Lifetimeresolved imaging can provide contrast between spatial regions on the diamonds that may contain differing concentrations of spectrally overlapping defects [8]. Thus, inhomogeneities in diamond samples, e.g. due to a change in growth conditions, can potentially be identified using a combination of spectral and lifetime imaging. For convenience we refer to luminescence lifetime imaging as FLIM (fluorescence lifetime imaging) [9] since we are utilising tools developed for fluorescence studies.
Complex time-resolved luminescence decay profiles are typically analysed using nonlinear iterative fitting techniques-usually to a multiexponential decay model. However, the assumption of a specific decay model may not always be justified and can be computationally intensive for large data sets. Phasor analysis offers a fast, fit-free method to analyse complex luminescence decay profiles and can be used to segment luminescence lifetime images based on the phasor components. It can also be used to determine the lifetimes of a biexponential decay if they are assumed to be constant across the image or region. Iterative fitting of complex decay profiles is still useful for quantification of population fractions of different components but can be sensitive to initial guesses for the component lifetimes. Here, we use phasor analysis to provide the initial guesses for iterative fitting, which can increase the speed and accuracy of the fitting. Reasonable initial guesses are particularly important when there are more than two spectrally overlapping lifetime components, or when little is known about the emitters present in a sample.
Phasor analysis of fluorescence decay data is a graphical approach utilising polar coordinates, which are often plotted on the complex plane. The application of this approach to analyse complex fluorescence lifetime data was first explored by Clayton et al [10]. and Redford et al [11]. using frequency domain lifetime measurements, and then applied to time-correlated single photon counting (TCSPC) fluorescence decay data by Digman et al [12]. Phasor analysis has been applied to fluorescence resonant energy transfer (FRET) experiments read out by FLIM, where emitting donor fluorophores have a reduced fluorescence lifetime based on their proximity to spectrally overlapping acceptor fluorophores [13]. To the best of our knowledge, phasor lifetime analysis has not previously been applied to diamond defect luminescence-for which time-resolved studies are frequently undertaken using the Hanbury Brown and Twiss method.
Spectrally-resolved lifetime measurements of diamond have previously been undertaken [8,14], where they have provided useful contrast between different emitters that exhibit spectral cross talk. Spectrally-resolved lifetime measurements can also be acquired as a function of spatial location in time-resolved hyperspectral microscopes, e.g. [15,16]. Further contrast has been demonstrated in spectrometers resolving emission with respect to excitation and emission wavelength, as well as fluorescence lifetime, e.g. [17,18]. Measurement of such excitation-emission-lifetime matrices require wavelength-tuneable pulsed excitation sources. This technique can also be extended to microscopy [19] and we recently demonstrated such an instrument applied to diamond luminescence [8]. The ability to vary the excitation wavelength enables studies of interactions between different defects. Here, we demonstrate how this can be used to offer insight into the origin of emitter quenching.
In previous studies of the quenching of charge states of nitrogen vacancies, it was shown that neutrally charged substitutional nitrogen (N s 0 ) can quench NV 0 6 , and that substitutional nitrogen (N s ) can quench NV − [4,5]. Quenching of NV − without subsequent quenching of NV 0 has not been reported, but NV 0 quenching without NV − quenching has been observed [5]. The natural lifetime of NV 0 has been reported to be ∼19 ± 2 ns [6], whilst for NV − , a double exponential decay profile has been reported with the m s = 0 lifetime being 13.7 ns and m s = ±1 being 7.3 ns, reflecting its more complex energy level structure and non-radiative transfer path [20,21]. The SiV − defect has previously been reported to present a monoexponential decay with a lifetime of 1.28 ns [22].
Here we report the adaptation and application of a previously developed instrument [8] to spectrally, and temporally analyse complex luminescence signals from nitrogen and silicon vacancies in a diamond sample, where quenching of NV − is also observed. These imaging capabilities could be useful for inspection of engineered quantum optical devices that will require detailed testing for faults and performance, noting that quantum effects will be sensitive to the specific local environment and distribution of defects rather than the bulk properties of the diamond. Automated imaging capabilities would also be useful to screen diamond samples for specific quantum electronic properties. In addition, the spectral versatility can enable the discovery and categorization of new defects, such as those synthesised through ion implantation and new growth methods [23].

Methods
The time-resolved hyperspectral microscope was previously described in detail in [8], and is schematically represented in supplementary figure 1 (stacks.iop.org/JPhysD/54/045303/ mmedia). Briefly, a supercontinuum excitation laser source (SC400-4, Fianium) generates ultrashort (picosecond) pulsed radiation at wavelengths from 400-2400 nm at a repetition rate of 40 MHz. A dichroic mirror (DMLP1000, Thorlabs) removes most of the radiation above ∼1000 nm, and spectral selection of the reflected supercontinuum beam is realised using a custom-built prism-based monochromator. The total available average excitation power is ∼800 mW in the 400-1000 nm band, where the spectral power density and spectral resolution vary as a function of wavelength, as illustrated in supplementary figure 2. Hyperspectral imaging was implemented using line illumination with the excited line of emission being imaged to the entrance slit of an imaging spectrograph (Acton SP300i, Princeton Instruments). A diffraction grating of 150 lines/mm provided ∼0.28 nm spectral resolution at an sCMOS camera (Zyla 5.5, Andor) on the imaging port, using a 10 micron entrance slit.
Time-resolved measurements were made using point illumination with a Peltier cooled photomultiplier (PMC-100-1, Becker and Hickl) and TCSPC card (SPC-830, Becker and Hickl) attached to the monochromator port of the imaging spectrograph. Inside the spectrograph, a motorized mirror could select between the two emission output paths. Excitation-emission-lifetime matrices were collected using a 50/50 non-polarising beam splitter (BSW26R, Thorlabs) in place of the microscope dichroic beamsplitter. For the emission filter, a long wavelength-pass linear adjustable filter (LWPLF, Delta Optical Thin Film) was used, which transmitted wavelengths longer than its cut-on wavelength, and which could be adjusted from 310 to 850 nm. Measurements at fixed excitation wavelengths were undertaken using a dichroic beamsplitter (525 nm cut-on, Semrock) and long pass emission filter at 514 nm (Razor Edge 514 LP, Semrock). The detected emission spectra have not been corrected for any spectral variation in quantum efficiency of the detectors, but variations in the excitation power as a function of wavelength (as shown in supplementary figure 2) have been taken into account. The data acquisition is controlled using a custom Lab-VIEW program.
For phasor analysis, the measured fluorescence decays were transformed to polar coordinates on the complex plane where the real and imaginary components are given by [12]: where G is the real component, S the imaginary component and I(t) is the intensity at time t [24]. The reference frequency ω is defined as ω = 2π Tp , where T p is the time between excitation pulses.
Frequency domain analysis of fluorescence decay profiles utilises two characteristic lifetime parameters based on the relative phase difference or the relative modulation depth of the fluorescence compared to the excitation signals, which are measured in frequency domain FLIM experiments. For monoexponential decays, the phase lifetime and the modulation lifetime are both equal to the lifetime that would be measured in the time domain. For complex decay profiles they are not equal, and each can be considered as a representative parameter of such a decay profile, analogous to an average lifetime in the time domain. In phasor analysis, an effective modulation lifetime can be calculated [11,12] from the polar co-ordinates using the equation: The modulation lifetime provides a convenient means to represent variations of complex fluorescence decay profiles that does not require fitting to an assumed model, which would be required if such variations were to be represented using average lifetimes calculated in the time domain.
The phasor analysis of the emission from each image pixel results in a point on a phasor plot, where the axes are the real and imaginary components G, S. These points should lie on, or within a 'universal circle'. Points on the universal circle correspond to image pixels presenting a mono-exponential decay profile. Image pixels presenting a double exponential emission decay profile map to points on a chord between two points on the universal circle. A more complex decay profile would be represented by points in a region within the universal circle bounded by chords linking the discrete components on the circle.
With any time-resolved measurement, it is necessary to account for the instrument response function (IRF), and this was measured using the attenuated reflection of the laser from a mirror at the sample plane onto the photomultiplier [25]. For excitation wavelength-resolved experiments, the arrival time of the excitation pulses at the sample (t 0 ) varies with wavelength due to path differences caused by the prismbased monochromator and properties of the supercontinuum source, and so a spectrally varying IRF was recorded. In the phasor calculation, the IRF of the system is taken into account using the following equation: G corrected − iS corrected = (G−iS) (GIRF−iSIRF) , where G IRF and S IRF are the real and imaginary phasor components of the IRF.
Due to practical limitations of our TCSPC system, the only way to measure a full 25 ns time window over a single pulse period is to measure over a larger window and discard some photodetection events. A window that is smaller than the pulse period can be measured without discarding photodetection events, but this results in an incomplete decay profile that impacts the phasor analysis [26]. In this work we employed a x2 frequency divider setting in the TCSPC electronics, which gives a useful measurement range of ∼46 ns. This range includes two consecutive decay profiles and the temporal delay is adjusted so that a full laser pulse period is collected over the first part of the measurement range. In order to make use of the photons detected in the second part of the measurement range, the data from the two parts of the measurement range are summed as shown in supplementary figure 3.
Iterative nonlinear fitting of the fluorescence decays was conducted using FLIMfit [27], which allows fitting to up to a 5-component exponential decay profile. The program was set to use maximum likelihood estimation (MLE), which provides less bias at low photon counts than traditional nonlinear least squares fitting based on χ 2 minimisation. The scale bar is 0.5 mm. The intensity contrast has been adjusted to enhance the visibility of signal from the CVD diamond region (see supplementary figure 4 for unsaturated image). The maximum diamond CVD emission was 3x lower than that of the HPHT substrate. The spectra in (b) are taken from the four ROIs shown in (a) and show that the CVD diamond and HPHT substrate have different emission spectra, originating from different defects/charge states. The emission spectrum from each ROI is shown using a separate vertical axis that is coloured to match the line colour for that ROI.
The model fitted to the measured decays is where β i is the population fraction of the ith decay component, τ i is its corresponding lifetime and I 0 is the intensity at time t = 0. Initial guesses for the iterative fitting were determined by fitting decays from large spatiallybinned regions or from the phasor analysis. Due to the relatively long lifetimes of diamond luminescence and high laser repetition rate, the measured decay profiles are truncated/incomplete: this is accounted for by FLIMfit for the iterative fitting and by choosing the 1st harmonic of the laser repetition frequency for the phasor analysis.
The background signal is dominated by the dark current from the photomultiplier, which is ∼160 counts per second divided across 1024 or 4096 time bins (depending on measurement settings). Due to the restriction imposed by the 25 ns laser pulse period and the long lifetimes of some components, the background cannot be accurately measured from the raw decay data and so is measured from a separate measurement with the excitation laser blocked. However, the typical photon counts measured from the sample are much higher (2-4 orders of magnitude) than the dark counts and so the background has a negligible impact on the phasor analysis and iterative fitting.
The synthetic single crystal diamond sample investigated had mass 10 mg and dimensions 1.52 × 1.92 × 0.75 mm 3 . It was grown by chemical vapour deposition on a type Ib HPHT substrate, a part of which remains attached. The sample was electron bombarded to produce vacancies and subsequently annealed. Previous studies of the luminescence from this diamond [8] indicated a predominance of nitrogen and silicon vacancies.

Hyperspectral Imaging
Hyperspectral imaging using a 4x 0.16 NA objective lens was undertaken over a ∼1.0 × 1.9 mm [2] region of the diamond encompassing the HPHT substrate (top of figure  1(a)), and the CVD grown diamond (bottom of figure 1(a)). Luminescence excitation at a wavelength of 500 nm (1.25 nm FWHM) was generated by spectral selection from a supercontinuum pulsed laser source using the prism-based monochromator and a line of illumination generated by the cylindrical lens, see supplementary figure 1. Figure 1(a) shows the calculated median wavelength of the emission (with brightness modulated by the intensity of the summed emission spectra). An emission spectrum can be extracted from any pixel or region of interest (ROI), of which examples are presented in figure 1 The spectrum from ROI 1 in the HPHT substrate shows emission from NV − , with its zero phonon line (ZPL) at 637 nm. In the CVD diamond, ROI 2 and ROI 3 show similar spectral features, with both presenting a strong NV 0 emission band with a ZPL at 575 nm. Several of the phonon side band features of NV 0 and NV − are also visible, which are discussed in detail in reference [28] for NV 0 . A contribution from SiV − with a ZPL at 737 nm is also apparent, with a ∼ 3x higher relative intensity for ROI 2 compared to ROI 3. This diamond sample was intentionally grown in bands with nitrogen dopants increasing further from the substrate, but with constant silicon dopant. This would explain this increase in SiV − emission intensity, as nitrogen acts as a donor [29] for the formation of the SiV − charge state. ROI 4 is located near the interface between the HPHT substrate and grown CVD diamond and exhibits emission from SiV − and a smaller contribution from NV − . There was no detectable emission from NV 0 , which suggests that NV 0 is not the preferred charge state at this interface. Formation of NV − in CVD diamond is important due to the quantum photonic applications of this defect, and so this phenomenon should be studied in more detail. These emission spectra suggest three main defects are present: SiV − , NV 0 and NV − . These measurements were taken at room temperature and so the phonon side bands of NV could obscure other spectral features. In all spectra, a small amount of 1st order Raman shift from the luminescence excitation was detected at 537 nm. Figure 2 presents data from emission restricted to a spectral band centred at 675 nm (2 nm bandpass), which should include emission from NV centres but exclude emission from SiV − defects. Stage scanning with spot illumination excitation at a wavelength of 500 nm was undertaken with a 40x 0.6 NA objective lens. The emission decay profiles for the FLIM map were analysed using phasor analysis and MLE fitting to a quadruple-exponential decay model using FLIMfit. Figure 2(b) presents the spatial phasor plot corresponding to the image shown in figure 2(a). A threshold of 10 000 photons was used to reject pixels where there were insufficient photon counts for accurate MLE fitting. The colour scale on the spatial phasor plot indicates the number of FLIM image pixels that map onto each location on the spatial phasor plot. Clusters of high density on the phasor plot can be analysed further, as shown in figure 4.

FLIM, lifetime fitting and phasor analysis
The modulation lifetime derived from the phasor coordinates is plotted in figure 2(c). This shows that the HPHT substrate has a significantly lower luminescence lifetime than the CVD diamond region. The mean modulation lifetime of the NV 0 defect that dominates the emission at the bottom of the CVD diamond region (blue dashed box) is 19 000 ± 393 ps, while the mean modulation lifetime of the NV − emission that dominates the HPHT substrate region (red dashed box) is 8420 ± 258 ps at this emission wavelength. The natural lifetime of NV − is shorter than NV 0 , but the mean modulation lifetime for NV − is shorter than expected.
When iterative MLE fitting was applied to the region of CVD diamond, the decay was well fitted by a monoexponential decay model with a lifetime of 19 000 ps and this is labelled τ 1 in subsequent analysis. The HPHT substrate region required a triple exponential decay model to obtain reasonably randomly distributed normalised fit residuals, and these components are labelled τ 2, τ 3, τ 4 . Figure 2(d) shows a plot of the intensity-weighted mean lifetime resulting from a maximum likelihood fit to a quadruple exponential These fixed lifetime components were determined by binning ten separate regions of interest in each of the HPHT substrate and CVD diamond and performing an MLE fit to either a single (for CVD diamond) or triple (for HPHT substrate) exponential decay model and calculating the mean values. The tolerances for these mean lifetime components were estimated to be: τ 1 = 19 000 ± 401 ps, τ 2 = 11 400 ± 188 ps, τ 3 = 4 400 ± 125 ps, and τ 4 = 615 ± 5.62 ps (mean ± standard deviation over ten regions to 3 significant figures). The phasor plot was used to inform the initial guesses for fitting the HPHT FLIM data, with the components τ 2 = 12 500 and τ 3 = 3 750 ps obtained from a line drawn through the data cloud inside the phasor circle and extrapolating to its intersection with the universal circle (dashed line on figure 2(b)). However, biexponential fits to decays from the HPHT substrate resulted in obvious structure in the fit residuals, which could be removed by changing to a three-component decay model, although this third component, τ 4, was not apparent in the phasor plot. The τ 4 component was determined to be 615 ps from iterative fits of data from ten regions. This triple exponential fit used the two initial guesses with component values of τ 2 = 12 500 and τ 3 = 3 750 ps from the phasor analysis, and the third component used a 1 000 ps initial value. These values were not fixed. For the calculation of the final intensity-weighted mean lifetime map, all four components were fixed in a quadruple-exponential decay MLE model, resulting in the intensity-weighted mean lifetime results shown in figure 2(d). Figure 3 shows the corresponding plots of the spatially varying population fractions, β i , of these components, showing that β 1 is only significant in the CVD diamond, and that the other three components are mostly confined to the HPHT substrate. The origin of these three lifetime components required further investigation.
A more detailed investigation of the spatial phasor plot is shown in figure 4. Regions on the phasor plot can be identified with corresponding pixels on the intensity and fluorescence lifetime maps. This means specific modulation lifetimes can be mapped to specific regions in the field of view. In figure  4(b), the pixels shown correspond to the red-bounded region on the phasor plot that sits on the universal circle. The mean modulation lifetime of 19 000 ± 580 ps from the phasor plot for this region agrees well with the result from FLIMfit and with the expected lifetime for NV 0 .
There is a localised cluster 'line' of points in the orangebounded region on the phasor plot that corresponds to pixels in the HPHT substrate region, as shown in figure 4(c). The modulation lifetimes derived from the phasor plot appear to vary by ∼1 ns in horizontal bands across the HPHT substrate region. A straight line drawn through these points on the phasor plot does not provide the lifetime components used in iterative fitting at the intersections with the universal circle, because the emission decay profile is more complex, as we will discuss below. These points lie within a triangle formed by chords between the three lifetime components obtained with FLIMfit; see dashed lines in (a).
The green-bounded region on the phasor plot corresponds to image pixels with complex emission decay profiles found only in a spatial region of CVD diamond grown at the interface with the HPHT substrate, as shown in the modulation lifetime map of figure 4(d). As indicated in the β maps of figure 3, this region seems to have a different composition from the rest of the CVD diamond; it is likely this originates due to NV − being the dominant charge state here, rather than NV 0 as shown in the hyperspectral image data. The initial growth of CVD diamond does not always occur under the same conditions as the growth of the rest of the sample, and this could cause a different uptake of defects.
When investigating the three lifetime components in the HPHT substrate of the diamond apparent in figure 4(c), it was noted that there are darker and lighter bands in the total emission intensity, as shown in figure 5(a). These bands were found to be correlated with the variation in modulation lifetime as shown in figure 5(b). This is confirmed in figure 5(c), which presents a plot of the total photon count and modulation lifetime across the line indicated by the black arrow in figures 5(a) and (b). A scatter plot of photon count and modulation lifetime is shown in figure 5(d), which suggests a linear relationship. The decrease in modulation lifetime with decreasing emission intensity could be explained by quenching of the NV − , which has previously been described [5]. To confirm that this correlation of lifetime and intensity was not an artefact of the TCSPC measurement, we looked for correlation between lifetime and intensity variations along the horizontal bands in the HPHT substrate and did not find a similar relationship (see supplementary figure 6). This spatial variation cannot be explained by the different lifetimes of the spin states of NV − , as the spin states are an intrinsic property of NV − , and would not be expected to exhibit any spatial variation.
FLIM data acquired for an emission wavelength of 737 nm (2 nm bandpass), i.e. centred on the ZPL of SiV − , under the same excitation wavelength and using the same microscope objective, is presented in figure 6. The phonon side band of NV 0 extends beyond 737 nm, as shown in figure 1(b) for the CVD diamond, and so this should be visible in the decay profiles. In the phasor plot in figure 6(b) there is a line from ∼1.1 ns on the universal circle to ∼19 ns on the other side of the circle, and this would correspond to a double exponential decay with these two lifetime components.  from changes in the β population fractions for the τ 1 and τ 4 decay components in the CVD diamond.
Similar to the analysis of emission detected at 675 nm, a triple exponential decay model using the same fixed lifetime components provided a reasonable fit to the emission from the HPHT substrate for luminescence detection at a wavelength of 737 nm. The phasor plot for emission at 737 nm does not directly indicate these three lifetime components, but this was supported by a detailed analysis of the fitted lifetime data and χ 2 and β maps for all five components. These maps can be found in supplementary figure 8. Thus, a total of five separate lifetime components are required to fit the luminescence decay profiles across the entire diamond sample. These components were used to calculate the mean and intensity-weighted mean lifetime maps in figure 6.
On the phasor plot in figure 6(b), there is a curved cluster of points extending from ∼1.1 ns on the universal circle to another cluster inside the universal circle (close to the phasor plot position for an 8 ns monoexponential decay), which represents the HPHT substrate emission. This arcing region on the phasor plot corresponds to NV − and SiV − emission, where the NV − concentration decreases as the distance from the substrate increases, resulting in a transition from the NV − to the SiV − dominant lifetime (see supplementary figure 9).

Spectrally resolved emission lifetime analysis
To explore the origin of the triple exponential decay in the HPHT substrate, the monochromator was used to measure a luminescence decay profile every 0.5 nm using a 40x 0.6 NA objective, and an excitation wavelength of 500 nm (1.25 nm FWHM). This can be plotted as an emission spectrumresolved phasor plot and the β values can be determined through fitting. This spectrally-resolved lifetime data was collected at two locations: ROI 1 representing the HPHT substrate, and ROI 2 the CVD diamond, as shown in figure 7(a). Figure 7(b) shows the emission spectra for each ROI. The ROI 2 spectrum from the CVD diamond is similar to the hyperspectral image spectra in figure 1(b), showing NV 0 and SiV − . ROI 1 from the HPHT substrate in figure 7(b) is similar to that of figure 1(b), except for the presence of a small component of NV 0 , in addition to that from NV − which provides the bulk of the emission. The origin of the NV 0 emission could be due to charge transfer effects caused by the higher intensity of the focussed excitation spot relative to line illumination with the same excitation laser power.
The decay data for the HPHT substrate ROI 1 shown in figure 7(c) was fitted to a quadruple decay model using MLE with lifetime components fixed to τ 1 = 19 000 ps, τ 2 = 11 400 ps, τ 3 = 4 400 ps, and τ 4 = 615 ps-as previously used for  Figure 7(c) shows that β 1 follows a spectrum which peaks at the ZPL of NV 0 . We therefore attribute β 1 to NV 0 as both the lifetime and spectrum are consistent with previously published data and this study. The remaining three β parameters (β 2 , β 3 & β 4 ) have similar spectra to one another, which suggests that these three decay components are all related to the same defect, NV − . β 1 (NV 0 ) is negligible above 650 nm, which is consistent with no NV 0 being detected in the HPHT substrate for the FLIM measurements when the emission bands were centred at 675 or 737 nm. As a four-component model was used to fit the emission over the entire sample, the NV 0 component was not found in the HPHT substrate for the FLIM maps. The variations in the different β components as a function of emission wavelength are depicted in figure 7(f), which shows the phasor modulation lifetime and intensity-weighted mean lifetimes from the MLE fits. The lifetimes obtained from the phasor and intensityweighted mean MLE fits are in reasonable agreement. The lack of variation in modulation and intensity-weighted mean lifetime in figure 7(f), also suggests only one emitting defect above 650 nm.
Luminescence emission from ROI 2 in the CVD diamond was fitted to a double exponential decay model of fixed lifetime components τ 1 = 19 000 ps and τ 2 = 1 100. This shows a shift in β for both components near 737 nm, see figure 7(c), indicating the presence of SiV − . This has the effect of reducing the modulation and intensity-weighted mean lifetimes in figure 7(f) with a minimum at 737 nm.
The CVD diamond spectral phasor plot in figure 7(e) of emission from ROI 2 shows a straight line between ∼1.1 ns and ∼19 ns. This suggests that there are only two emitting defects, NV 0 and SiV − , in the CVD diamond for an excitation wavelength of 500 nm. The HPHT substrate spectral phasor plot in figure 7(d) of emission from ROI 1, also shows a straight line for the bulk of its emission. This indicates there are two defects in the HPHT substrate diamond with strong emission for 500 nm excitation, attributed to emission from NV 0 and NV − . We note that, if the line is extrapolated to the intersection with the universal circle, the indicated short lifetime component does not lead to a good fit if included in the fitted MLE model. This requires further investigation.

Spectrally resolved emission and excitation lifetime analysis
Photoluminescence excitation-emission-lifetime (PLEEL) spectroscopy is a powerful tool to analyse complex samples. In PLEEL, the excitation wavelength is varied, and the emission spectrum analysed to determine the absorption and emission spectra of all emitters in the measured spectral range. We can also use emission lifetime contrast, e.g. to better distinguish spectrally overlapping components. Here, we implement PLEEL with the tuneable supercontinuum excitation source, and detect via the monochromator using TCSPC. An integration time of 0.5 s was used to measure each decay profile and we acquired a full excitation emission lifetime matrix over relatively wide spectral regions and with relatively high spectral and temporal resolution. The emission wavelength range covered 450 to 780 nm, and in excitation from 430 to 750 nm. A linear variable long wavelength-pass filter (LVLPF) was used as a tuneable emission filter to block the excitation radiation as the excitation source was tuned. A 50/50 beam splitter replaced the usual dichroic beamsplitter to allow broad tuning of the excitation and emission wavelengths. The spectral separation between the edge of the filter cut-on and the excitation laser line was set to 20 nm. Data was collected using a 4x 0.16 NA objective lens and a data point was collected every 2 nm in emission and 2.5 nm in excitation wavelength. Following phasor analysis and fitting in FLIMfit, the luminescence lifetime data were plotted as a function of excitation and emission wavelength.
Results from phasor analysis of this lifetime-excitationemission data are shown in figure 8 and results from MLE fitting to a quadruple-exponential decay model are shown in figure 9. The emission was measured with a spectral resolution of 2.5 nm and the spectral FWHM of the excitation varied with wavelength as shown in supplementary figure 2. Due to imperfect blocking of the excitation radiation by the filters employed, there was some bleed-through which resulted in a straight-line artefact with unity gradient in the lifetimeexcitation-emission plots of figures 8 and 9, and is also indicated by the blue region on figure 8(c). The modulation lifetime data derived from phasor analysis (figure 8) and intensity-weighted mean lifetime derived from MLE fitting to a four-component decay model ( figure 9) are consistent. Figure 8 presents data collected from the location selected on the HPHT substrate, shown in figure 8(a). Figure  8(b) shows an excitation-emission-resolved phasor plot where the colour scale bar represents the number of pixels in the excitation-emission map with that phasor coordinate. The excitation-emission phasor plot indicates, similarly to figure 7(d), a straight-line distribution originating at ∼19 ns on the universal circle (corresponding to NV 0 ) and extending to a cluster inside the universal circle, representing the complex decay profile of NV − emission.
The primary features in the excitation-emission photon count image in figure 8(c) are the laser light bleed-through (blue region) and the excitation/emission spectra of NV − (yellow region), with its ZPL excitation band indicated by the red arrow. Resonant excitation of NV − is apparent at an excitation wavelength of 637 nm (green arrow) and its corresponding emission band. In figure 8(e), the modulation lifetime excitation-emission plot shows three separate lifetime regions.  One region can be attributed to NV 0 , (indicated by dashed white rectangle in figure 8(e)), which is excited below 580 nm, and extends to ∼625-650 nm in emission, consistent with the data in figure 7(c). At 575 nm excitation (orange arrow in figures 8(e) and (f)) there is an increase in the emission lifetime that can be attributed to resonant excitation of NV 0 , which can be supported by the small peaks corresponding to phonon side band emission.
The presence in figures 8(c)-(f) of only two narrow regions (denoted by green, red and orange arrows) corresponding to resonant excitation, suggests there are only two emitting species, namely NV 0 and NV − . The more uniform region (purple pixels in figure 8(e) and blue pixels in figure 9(a)) for emission above 650 nm, cannot be attributed to NV 0 , due to the lack of significant emission above this wavelength as shown in figure 7. Therefore, this emission is attributed to NV − , and the range of lifetimes observed may be explained by quenching, perhaps due to charge effects, or interactions between both charge states of NV. This requires further investigation.
This detailed analysis suggests that the bulk of the CVD diamond has two defects; NV 0 and SiV − , and the HPHT substrate contains primarily NV − . The NV − in the HPHT substrate has a complex decay profile, and when considering the reduction in mean lifetime in regions of lower intensity, it is likely that this defect is being quenched by an unknown spatially varying defect. This defect is unlikely to be the same defect that causes quenching in NV 0 , as NV 0 in this HPHT substrate displays a mono-exponential decay with a lifetime similar to literature [6]. This quenching is not unique to the HPHT substrate as the CVD/HPHT diamond interface also exhibits quenching of NV − . Extracted excitation and emission spectra with lifetime can be seen in supplementary figure 10, where the peak of the NV 0 excitation spectra is shown coinciding with the increased lifetime peak of NV − , providing evidence that an interaction is taking place. This quenching is manifested both as spatially varying as shown in figure 5, and spectrally varying as shown in figure 8. The complexities of this reduction in lifetime could therefore suggest multiple quenching mechanisms.

Conclusions
Using the multidimensional luminescence microscope we developed, we have demonstrated that complex decay profiles of luminescence in diamond can be analysed using the phasor approach, which can inform the lifetime fitting of such data, e.g. to determine population fractions, β, and assist interpretation. Hyperspectral imaging was initially used to determine the wavelengths of interest for FLIM, and FLIM data was subsequently analysed with phasor plots and iterative fitting. Clusters of data points on the phasor plots can be used to segment the lifetime images, allowing spatial variations in lifetime components to be easily identified. This was seen with the lifetime contrast corresponding to the striations in CVD diamond attributed to SiV − and NV 0 , attributed to the known preferential incorporation of nitrogen into different growth planes.
NV − emission in the HPHT substrate was found to be complex, and modulation lifetime contrast from phasor analysis, emission lifetime analysis, and excitation emission lifetime matrices, suggest that NV − emission is being quenched, resulting in spatial variation of the NV − lifetime. Such quenching was not observed for NV 0 in this sample. The lifetime-excitation-emission data suggest that NV 0 could be involved in the quenching of NV − , noting that an increased NV − lifetime was observed where NV 0 was excited. No other emission, such as from H3, was observed in the HPHT substrate.
The flexibility of this instrument for multidimensional luminescence imaging and spectroscopy, combined with both phasor and nonlinear fitting approaches to lifetime analysis, enables complex spectroscopic signatures to be explored and mapped in detail. This can reveal spatially varying phenomena such as quenching of NV − emission, which could be important for quality assurance in the manufacture of quantum optical devices, and which will support better understanding of mechanisms for quenching of emissive defects in diamond.