Heralded spectroscopy reveals exciton-exciton correlations in single colloidal quantum dots

Multiply-excited states in semiconductor quantum dots feature intriguing physics and play a crucial role in nanocrystal-based technologies. While photoluminescence provides a natural probe to investigate these states, room temperature single-particle spectroscopy of their emission has so far proved elusive due to the temporal and spectral overlap with emission from the singly-excited and charged states. Here we introduce biexciton heralded spectroscopy, enabled by a single-photon avalanche diode array based spectrometer. This allows us to directly observe biexciton-exciton emission cascades and measure the biexciton binding energy of single quantum dots at room temperature, even though it is well below the scale of thermal broadening and spectral diffusion. Furthermore, we uncover correlations hitherto masked in ensembles, of the biexciton binding energy with both charge-carrier confinement and fluctuations of the local electrostatic potential. Heralded spectroscopy has the potential of greatly extending our understanding of charge-carrier dynamics in multielectron systems and of parallelization of quantum optics protocols.


Introduction
Over the past three decades, numerous types of semiconductor nanocrystals with varying compositions, shapes, sizes and structures have been fabricated and studied, 1-3 with some even making their way into mass produced consumer products. 4 Since the energy of a charge carrier in a nano-confined solid is quantized, nanocrystals are often referred to as 'artificial atoms'. In a further analogy to atomic physics, photo-exciting such a quantum dot (QD) generates a Hydrogen-like electron-hole state, an exciton, which is typically bound even at room temperature, due to the increased Coulomb interaction. However, unlike atoms and molecules, semiconductor nanocrystals include another readily-excited manifold of states, multiexcitons; multiple electron-hole pair states. 5 In the lowest energy multiexcitonic state, the biexciton (BX), a strong exciton-exciton interaction is imposed by the confining potential of the nanoparticle. In the resulting energy ladder of ground, single exciton (1X) and BX states, shown in Figure 1a, the BX state energy is somewhat offset from twice the energy of the 1X state, by the BX binding energy (ε b ). 5 This binding energy is considered positive (attractive interaction) when E BX < E 1X , where E k is the energy difference between state k and the state right beneath it in the ladder.
Cascaded relaxation from the top to the bottom of this ladder can yield a pair of photons; the first around E BX and the second around E 1X . Such an emission from self-assembled Indium Arsenide (InAs) QDs, for example, is a leading candidate for efficiently generating on-demand entangled photon pairs. 6 On the other hand, avoiding the excitation of the BX state is key for using the same type of QDs as high-purity singlephoton sources. 6 More conventional light-based applications that stand to vastly benefit from the incorporation of colloidally synthesized QDs, such as light emitting diodes (LEDs), 7,8 lasers, 9,10 displays 4 and photovoltaics, 11 also require characterization and control of the energy and dynamics of the BX state. For instance, non-radiative Auger recombination often dominates the BX relaxation dynamics. 12 Therefore, to achieve low threshold lasing from QDs, sophisticated heterostructures have been designed to either increase the BX energy via Coulomb repulsion (avoiding its occupation), 13 or to reduce the Auger recombination rate. 14 Furthermore, the Auger induced low BX quantum yield sets a saturation boundary for the achievable light fluence of nanocrystal-based LEDs. 15 The above-mentioned interest in the BX-1X ladder inspired substantial spectroscopic efforts to character- Figure 1: Obstacles in measuring the 1X-BX energy ladder. a) The energy diagram of the ground, single exciton (1X) and biexciton (BX) states in a nanoparticle. The energy difference between the BX and 1X states (E BX ) is smaller than the difference between the 1X and ground states (E 1X ) by the biexciton binding energy (ε b ). b) Schematic of the thermal broadening of spectral lines (intensity normalized for clarity). At room temperature the 1X-ground (blue) and BX-1X (red) transitions emission lines are broadened to approximately ∼k B T ≈ 26 mev > ε b . As a result, the two emission spectra substantially overlap. c) Schematic dependence of the 1X energy on nanoparticle size. An ensemble measurement (with a ±5% size variance) at room temperature roughly includes a mixture of the depicted spectra. d) Scheme of spectral drift in the emission lines of a single nanocrystal throughout the measurement time. Apart from the stochastic random walk of the spectral lines (diffusion), discrete spectral jumps typically accompany blinking events.
ize the BX binding energy in multiple material systems such as III-V, [16][17][18] II-VI 19,20 and lead halide perovskite 21,22 semiconductor nanocrystals as well as atomically thin films of transition metal dichalcogenides (TMDC). 23 However, conditions under which the BX-1X ladder can be directly probed are restrictive. While at cryogenic temperatures the very stable and narrow 1X and BX emission lines of a single self-assembled InAs or colloidal cadmium chalcogenide QDs can be discerned, 17,18,24 in most other cases this was not achieved due to several fundamental limitations. First, at room temperature, both emission lines are thermally broadened well-beyond typical BX binding energies (Figure 1b). Spectral features are further broadened in ensemble measurements due to nanoparticle inhomogeneity ( Figure 1c). Additional inhomogeneous broadening is caused by spectral fluctuations (Figure 1d). In many nanocrystals, this includes not only the spectral diffusion due to spurious electric fields, but also spectral jumps due to other states contributing to fluorescence intermittency. 25,26 More indirect methods to probe the BX state rely on power-dependent measurement of photoluminescence (PL) (either in a time-resolved or quasi-continuous-wave manner) 20,22,23,27 or of transient absorption. 21,28 While careful modeling and analysis of these measurements provided important spectroscopic information, it has often led to large variance in BX binding energies measured in different studies. For example, while some works found very high (∼100 meV) BX binding energies in CsPbX 3 (X = Cl, Br, I) nanocrystals, 21 recently a substantially more stringent bound on its magnitude (|ε b | < 20 meV) has been reported. 22,29 Such discrepancies are mainly due to the difficulty in modeling the different mechanisms that affect PL and absorption at high excitation powers, such as charging, oxidation, blinking and photo-induced damage. 22 Furthermore, these methods demand disentangling the spectral contributions of the BX and 1X states, which becomes increasingly difficult to perform when these features strongly overlap for small BX binding energies.
While isolating the BX state in the spectral domain alone is a convoluted task, isolating it in the time domain is conceptually simple. A detection of a photon pair emitted from a single nanocrystal (following a short excitation pulse) pinpoints a relaxation cascade: 17,18 first from the BX to the 1X state and then from the 1X to the ground state ( Figure 1a). Separating the two consecutive emissions according to their detection times, on a nanosecond scale, and measuring their respective spectrum, facilitates a direct measurement of E BX and E 1X for a single nanoparticle. However, currently available instrumentation does not enable practical implementation of this simple scheme. Namely, a standard spectrometer cannot provide the required temporal resolution since it relies on cameras with a maximal frame rate of ∼10 3 fps. This can be addressed by replacing the camera with a monolithic single-photon avalanche diode (SPAD) array; a technology that has achieved a considerable performance boost over the past decade. 30,31 In such a novel spectrometer, termed here spectroSPAD, the spectral information of single photons can be measured and correlated with sub-nanosecond temporal resolution.
Here we present the spectroSPAD system and use it to measure spectral correlations in BX-1X emission of single CdSe/CdS/ZnS core/shell/shell QDs with meV precision. Thanks to the temporal resolution and high sensitivity of the spectroSPAD, our measurement scheme overcomes all of the aforementioned obstacles (thermal broadening, spectral diffusion, blinking and low BX yield) and easily separates the BX emission from the misleadingly similar trion ('grey') charged state emission. While for the particular sample under study the average value of the binding energy is ∼6 meV, we disentangle the inhomogeneous size effect and show that its value in individual QDs correlates with the 1X band edge transition energy. Furthermore, we follow the temporal fluctuations of the BX binding energy for a single nanocrystal and find that those correlate with the spectral diffusion of the 1X transition.  Figure 2: Sketch of the spectroSPAD experimental apparatus. A microscope objective is used to focus a pulsed laser beam on a single QD and collect the emitted fluorescence. At the microscope output, light is passed through a spectrometer setup with a linear SPAD array detector. Inset (i) is an optical image of part of the detector array. Each blue circle represents a single pixel. Scale bar is 30 µm. Insets (ii) and (iii) show two possible analyses of a single nanocrystal spectroSPAD measurement. In (ii), detections are binned according to detection time in the millisecond scale (horizontal) and photon energy (vertical); whereas in (iii) the horizontal axis represents the detection temporal delay from the preceding excitation laser pulse in nanosecond scale. Color-scale corresponds to the number of counts in each temporal-energy bin. Optical elements: objective (L1), tube (L2), collimating (L3) and imaging (L4) lenses; blazed grating (BG).

Results and Discussion
Apparatus. In recent years, considerable efforts were invested in the design of time-resolved light spectrometers with high sensitivity. 32-36 As a replacement for the standard CCD camera, different research groups adopted photo-multiplier tube (PMT) arrays, 33 superconducting nanowire single photon detectors (SNSPDs) 35-37 or SPAD arrays. 32,34,38,39 While these implementations harbor great potential for applications such as Raman spectroscopy and on-chip quantum communications, none is able to provide the combination of high overall detection efficiency, low dark counts, and parallel time and spectrum detection at single-photon level. The spectroSPAD spectrometer ( Figure 2) achieves precisely that by employing a high-performance linear SPAD array as a detector in a Czerny-Turner spectrometer. While a detailed description of the experimental setup is given in Supporting Section S1, we provide here a brief account. A microscope with a high numerical aperture objective is used to focus pulsed laser illumination on a single QD, and to collect epi-detected fluorescence. This signal is spectrally filtered from the excitation laser with a dichroic mirror and a dielectric filter (not shown), and imaged by a second lens. This image serves as the input for a spectrometer setup -a 4f system with a blazed grating at the Fourier plane. At the output image plane of the spectrometer, a monolithic linear SPAD array is placed, such that each pixel is aligned with the image of a different wavelength range. The detector's photon detection efficiency (PDE) is ∼8% at 620 nm and ∼11% at 530 nm (this can be improved, see Discussion) and the median dark count rate (DCR) is ∼33 counts per second per pixel. We analyze the signal of only 40 out of the 512 pixels available in the array, thereby spanning approximately 80 nm around a center wavelength of 620 nm. This results in a spectral resolution of ∼2 nm (6−7 meV). Single-photon detections are time-tagged by an array of 64 time-to-digital converters (TDCs) programmed in the firmware of a field-programmable gatearray (FPGA). The TDCs are synchronized with the excitation laser. Finally, time and wavelength tagged data is analyzed with a dedicated MATLAB script.
Insets (ii) and (iii) show possible visualizations of fluorescence data collected by the system from a single QD, as 2D detection histograms (see Supporting Section S2 for QD details). The spectrum over time is seen in (ii), where the time of detection spans the horizontal axis (10 ms time-bins) and energy the vertical axis (6−7 meV energy-bins). The effects of thermal broadening, spectral diffusion and blinking, discussed above, can be clearly observed through the width of the spectral peak at each time-bin (35−50 meV FWHM), the temporal jitter of the spectral peak position and the variation of emitted intensity, respectively. To achieve spectrally dependent fluorescence decay curves, shown in (iii), the same dataset is analyzed by binning the detections according to their delay from the preceding excitation pulse. Note that a full horizontal binning (FHB) of either histogram is equivalent to a standard spectrometer measurement, and a full vertical binning   (FVB) to a standard analysis from a single-SPAD measurement.
1X spectral dynamics. Prior to analyzing the BX state, it is important to first study the 1X state. Indeed, simultaneous acquisition of both temporal and spectral data enables a more in-depth analysis of the fluorescence spectral dynamics. Figure 3 employs such an analysis to identify and quantify the spectral broadening effects shown in Figure 1b,d for a single QD measurement. Figure 3a shows the total fluorescence intensity collected over all array pixels, at 1 ms time-bins. The intensity trace features a characteristic blinking behavior -stochastic switching between a bright ('on'), dim ('grey') and dark ('off') fluorescent states. 40,41 The presence of these three states is clearly evident in Figure 3b, a histogram of fluorescence intensities over the entire 5-minute measurement. Some time-bins are classified as 'mixed', interpreted as time-bins where the QD spent comparable time in different states. Analyzing the spectrum of each intensity state (Figure 3c), reveals a clear red-shift of the 'grey' state spectra, as well as an order-of-magnitude shorter fluorescence lifetime (see Supporting Section S3). This supports an identification of the 'grey' state as emission from a charged exciton state (identified as a negative trion in past work 41 ). According to ensemble measurements, the BX state is also expected to present a shorter lifetime and shifted emis-sion compared with the 'on' state, making it difficult to differentiate the two contributions. In addition to blinking, the emission also exhibits spectral diffusion. Figure 3d shows the spectral evolution of the neutral 1X emission ('on' state) over time. Each dot represents the momentary mean photon energy over a 1 ms time-bin ( E 1X 1 ms ), colored according to the local density of data-points for clarity. The red line represents a Gaussian-weighted (σ = 10 s) moving average of these values ( E 1X 10 s ). This smoothed trend shows a gradual wavelength change, predominantly towards shorter wavelengths, possibly due to oxidation. 42 The fast spectral diffusion dynamics are evident in the distribution of E 1X 1ms around this moving average, ∆E 1X E 1X 1 ms − E 1X 10 s . These faster dynamics are typically attributed to rapid fluctuations in the local electrostatic potential, leading to a shift in emission energy according to the quantum confined Stark effect. 43 The above results emphasize the difficulty in isolating the BX state emission. Namely, its rather weak contribution is overshadowed by spectral broadening, and especially by the 'grey' state emission, which overlaps with it in both spectral and temporal domains. As a result, even a comprehensive analysis of the 2D lifetimespectrum data (see Supporting Section S4) was unable to resolve the BX state spectrum. 1X spectral peak (eV)  Figure 4: Heralded spectroscopy. a) 2D histogram of photon pairs following the same excitation pulse, according to the energy of the first (E BX ) and second (E 1X ) photons (vertical and horizontal axes, respectively), over a 5-minute measurement. Dashed green line serves as a guide to the eye, marking same energy for both photons (undetectable by the system) b) Spectrum of the BX (red dots), 1X (blue rings), and all 'on' state detections (grey area, normalized). Solid red and dashed blue lines are fits of the BX and 1X spectra, respectively, to a Cauchy-Lorentz distribution. Binding energy, estimated as the difference between BX and 1X spectral peaks, is ε b = 9.3 ± 1.0 meV. c) Binding energy as a function of 1X spectral peak for 30 QDs. Error bars depict 68% confidence intervals.
Heralded spectroscopy. To directly probe the BX emission, pairs of photon detections following the same excitation pulse are post-selected. Such paired events are the result of an excitation to the BX state, and two subsequent radiative relaxations ( Figure 1a). We note that, due to the low quantum yield of the BX (∼9%, see Supporting Section S5 and reference 44 ), this is not the most probable route for relaxation from the BX state. Yet, its occurrence provides sufficient signal for our analysis. Applying this post-selection (see details in Supporting Section S6) to the 5-minute single-QD acquisition shown in Figure 3, yields ∼1.4 × 10 3 pairs over 1.5 × 10 9 excitation pulses (in agreement with theory, see Supporting Section S5). The 2D spectrum of photon pairs, showing the distribution of the energy of the first emitted photon as a function of that of the second, is shown in Figure 4a. The distribution is clearly centered below the diagonal, indicating BX binding. Note that events where both photons of a single cascade impinge on the same pixel are not detected by the system due to pixel dead time (∼25 ns). Figure 4b highlights the first insight that can be derived by such an approach -the BX spectrum (red dots, FHB of panel a) is red-shifted with respect to the 1X spectrum (blue rings, FVB of panel a). The agreement of the 1X spectrum with the overall spectrum of the 'on' state (grey area), corroborates this distinction. The BX binding energy for this particular QD, estimated as the difference between the BX and 1X spectra peaks (extracted by fitting a Cauchy-Lorentz distribution, shown as lines in Figure 4b), is ε b = 9.3 ± 1.0 meV (68% confidence interval). We note that the identification of the 1X and BX spectral peaks is done here without ambiguity. While previous studies required a power dependence series to correctly assign the 1X and BX states, 16,20 heralded spectroscopy obviates this requirement. More importantly, this approach super-resolves the few-meV separated 1X and BX spectral peaks despite their ∼50 meV FWHM, and clearly distinguishes between the overlapping BX and 'grey' state emission. In fact, owing to the unprecedented sensitivity of this method, measuring QDs featuring lower ε b than almost all previous measurements of II-VI semiconductor QDs did not incur any additional challenge (see Discussion section below).
Thanks to the single-nanocrystal nature of this method, it is not limited to measuring ensemble averaged properties, but can also observe their distribution within the ensemble. Figure 4c shows that the BX binding energy increases with the 1X spectral peak position for 30 QDs taken from the same sample. This can be explained as a result of variation in the physical size of the synthesized QDs. For the QDs investigated in this work, a higher energy 1X spectral peak is likely associated with a thinner CdS shell. A thinner shell also corresponds to further confinement of the electrons in the core and an increased Coulomb interaction between charge carriers, leading to a higher BX binding energy. This trend is in agreement with ensemble measurements for CdSe/CdS seeded nanorods. 27 2 meV ∆E 1X window. The red line represents a linear fit of these medians, emphasizing the positive crosscorrelation of ε b and ∆E 1X . The slope of this line is 0.59±0.08 (68% confidence interval); i.e. for each 2 meV red-shift (blue-shift) of the 1X emission spectral peak, the binding energy is lower (higher) by roughly 1.2 meV. Figure 5b, a histogram of the ε b median slope values for 30 QDs, shows that this positive correlation is evident for all QDs measured. This result suggests that the BX binding energy, much like 1X emission, is subject to the quantum confined Stark effect. A higher local field associated with a red-shift of the 1X emission, is correlated with lower BX binding energy (weaker attraction). This can be attributed to the spatial separation of holes and electrons induced by the external field. Notably, at high enough fields, the sign of the BX binding energy flips, indicating repulsive interaction of the excitons. This observation agrees with past results on the effect of charge separation in type-II QD heterostructures on the BX binding energy. 45 Furthermore, it strengthens the assertion that spectral diffusion indeed originates from fluctuations in the local electrostatic potential.
Discussion. Due to the practical and fundamental importance of the BX state in II-VI semiconductor nanoparticles, extensive experimental work has been dedicated to the estimation of its energy, 5,19,20,27,[45][46][47][48][49] with ε b values spanning -100's to +10's meV (see Supporting Section S7). In core/shell CdSe/CdS nanocrystals, ε b has been shown to transition continuously from positive (attractive) to negative (repulsive) values when tuning the core or shell diameter. 27 This is indicative of a transition from type-I to type-II or quasi type-II architecture, where electrons and holes are separated to different layers of a heterostructure. 45 The few-meV ε b values measured in this work are in good agreement with particles near this transition. It is worth noting, that for such low ε b , ensemble measurements become increasingly challenging, as they demand resolving two highly overlapping spectral peaks. Ensemble techniques are, therefore, more readily applied to measure particles exhibiting tens of meV BX binding energies, which are indeed reported more often for II-VI nanocrystals (see Supporting Section S7). More importantly, ensemble methods require the delicate analysis of a power-series measurement, and therefore prone to systematic biases due to power dependent charging and absorption cross-section heterogeneity. In fact, such heterogeneity in the sample was suggested as the source of recent controversy regarding the BX binding energy in perovskite nanocrystals. 21,22 The background-free (isolated BX spectrum) and singleparticle nature of heralded spectroscopy afford an unprecedented, sub-meV, sensitivity in BX binding energy measurement at room temperature, while at the same time circumventing the above mentioned biases.
While enabling previously inaccessible measurements, at the present there are two limiting factors that should be taken into account when implementing onchip SPAD arrays in correlative-spectroscopic setups. First, the detection efficiency (∼10% PDE) is still low compared with common camera alternatives. However, this is in part due to the circuit board design implementation used here, dictating signal collection from only every other pixel and a relatively low gain voltage on the diodes. The next system iteration (currently in development) specifically tailored for the purpose of spectroSPAD, is expected to achieve >30% photon detection efficiency, similar to SPAD arrays implemented with the same technology. 44 This will give rise to an order of magnitude increase in the correlation signal level due to its quadratic dependence on the detection probability. Moreover, detection efficiency is complemented by the very low noise level of state-of-the-art SPAD arrays. Unlike commonly used cameras, these arrays avoid readout noise and feature median dark count rates well below 100 counts per second per pixel. 50 A second factor to consider, common in detector arrays and specific to photon correlation analysis, is inter-pixel crosstalk. Due to the close-packing of pixels, a detection in one pixel has a small probability to lead to a false detection in a neighbouring pixel, and hence a false photon pair. Any bias due to this effect is mitigated here by a combination of the chip design and temporal gating (see Supporting Section S6), bringing the crosstalk probability down to ∼10 −5 , and of a statistical correction as described in reference 44 .

Conclusions
Heralded spectroscopy of BX emission cascades enabled us to perform direct observation and unambiguous identification of emission from multiply-excited states of single QDs at room temperature. In addition to avoiding the pitfalls of indirect and ensemble approaches, by separating the BX and 1X emission in the time domain we greatly extend the range of nanomaterials in which we can observe the BX state, and allow new insights into exciton-exciton interaction within the single nanoparticles. Our study reveals a positive correlation between exciton-exciton attraction and tighter chargecarrier confinement of the single QDs. We also unveil a fluctuation of this attraction strength, correlated with the fast fluctuations of the local electrostatic potential, and significant enough to lead to exciton-exciton repulsion. These capabilities represent a new tool to probe QD physics, and can lead to better design of QD-based technologies where multiexcitonic states typically play a major role.
All this is enabled by constructing the spectroSPAD -a SPAD-based correlative spectrometer extending the temporal resolution limit of standard spectrometers by several orders of magnitude. This tool is not only useful for probing the physics of charge-carrier dynamics, but can also address current challenges in quantum optics and communication.
Supporting Information Details of the spectro-SPAD system and the linear SPAD array; details of the quantum dots used in this work including synthesis, sample preparation and excitation saturation estimation; fluorescence decay lifetime by intensity state analysis; 2D lifetime-spectrum analysis; BX quantum yield estimation; heralded spectroscopy analysis and correction details; published values of BX binding energies in II-VI nanocrystals.

S1 spectroSPAD system details
This section presents a detailed overview of the experimental apparatus, and further technical details of the linear SPAD array detector at its core.

S1.1 spectroSPAD setup overview
The system is built around a commercial inverted microscope (Eclipse Ti -U, Nikon). A pulsed diode laser (470 nm, 5 MHz, LDH-P-C-470B, PicoQuant) is focused through an oil immersion objective (x100, 1.3 NA, Nikon) on a single QD. Illumination power density at the sample plane is ∼140 W/cm 2 leading to ∼66% probability to excite at least one exciton (per pulse, see subsection S2.2). The same objective is used to collect the emitted fluorescence, while back-scattered laser light is filtered by a dichroic mirror (505 LP, Chroma) and a long-pass dielectric filter (488 LP, Semrock). At the output of the microscope, the spectrometer consists of a collimating lens, a blazed grating (235 g/mm, 5.06°b laze, 53-*-790R, Richardson) and an imaging lens, resulting in 3.9 × 10 −5 reciprocal linear dispersion and ∼6 A spectral resolution (FWHM). At the spectrometer output image plane, a 512 pixel linear SPAD array (an upgraded version of the sensor described in reference 1 , see details in subsection S1.2), is placed such that the active pixel pitch is ∼2 nm in wavelength (every second pixel is active). An FPGA with an implemented TDC array (synchronized with the laser) assigns timestamps and pixel addresses to single detections in 40 pixels of the array. The trace of detections is then analyzed by a dedicated MATLAB script, implementing temporal and intensity corrections (see section S6) and the analysis schemes. S1.2 SPAD array technical details a b Figure S1: The linear SPAD array. a) Electrical circuit of a SPAD array pixel. b) Optical image of the detector array, mounted with microlenses. Each blue square represents a single pixel. Scale bar is 100 µm.
Since the creation of the first single-photon avalanche diode (SPAD) in complementary metal-oxide semiconductor (CMOS) in 2003, research in the field of SPADs and SPAD image sensors has led to the creation of the first integrated array in 2004, followed by a wide range of SPAD based image sensors with advanced functionality and continuously increasing speed.
The sensor used in this paper comprises an array of 512 SPAD pixels based on the design proposed in reference 2 . Each pixel comprises a SPAD quenched and recharged passively through a poly resistor. The SPAD is interfaced to the exterior of the chip through a circuit described in Figure S1a, including capacitive decoupling, a clamp to V dd , and a low-threshold buffer. The purpose of this circuitry is to ensure low threshold of detection of the avalanche, thus optimizing jitter while controlling noise. Quenching resistor R q is designed to present a sufficiently high impedance to the anode of the SPAD, while minimizing the avalanche current, so as to control the overall power consumption of the chip. V op is set to V BD + V EX , where V BD ≈ 24.4 V and V EX ≈ 1.6 V are the breakdown and excess bias voltages, respectively.
The chip was mounted directly on a board with the SPAD outputs wire-bonded and connected to a fieldprogrammable gate array (FPGA), which hosts the TDCs that enable the time characterization of the response. The TDC array is an improved version of the earlier implementation detailed in reference 1 . By multiplexing over 64 TDC channels, up to 256 pixels of the array can be temporally correlated. Extending the system to two FPGA boards, enables simultaneous read out from all 512 pixels. However, to reduce DCR, only 40 pixels were used in this work. The DCR reduction is achieved both by collecting DCR from fewer pixels, and by avoiding the few 'hot pixels' in the array that feature exceptionally high DCR (typically 'hot pixels' are defined as those that feature DCR at least two orders of magnitude higher than the median; here the 'noisiest' pixel DCR was just a factor of 3 above the median (104 cps)). The choice of 40 pixels resulted in overall DCR that is about an order of magnitude lower than the detected 'on' state emission in intensity measurements. After temporal gating (see section S6) the number of DCR induced pairs was also about an order of magnitude lower than the total number of detected pairs. Pixels exhibit an average jitter of 105 ps (FWHM) and a median DCR of 33 cps, their native fill factor (without microlenses) and pitch are 25.1% and 26.2µm, respectively. Figure S1b shows a detail of the SPAD array, including microlenses. Microlenses were deposited on the chip to enhance effective fill factor and thus overall photon detection efficiency (PDE).
The recorded temporal response after the TDC (Figure S2) features two peaks: A narrow main peak accounting for ∼90% of the counts, and a secondary broader peak delayed by ∼3.5 ns, accounting for the remaining ∼10%. The presence of the secondary peak is specific to the particular implementation used in this work, and is not generally evident in similar detector arrays. Preliminary results suggest that this irregular temporal response can be eliminated by better firmware design in the next system iteration (currently in development). The QDs investigated in this work featured 1X decay lifetimes significantly longer than this artifact (see section S4). Hence, any ambiguity in order of arrival can be negated by temporal gating (as detailed in section S6) with minimal loss of signal.

S2 Quantum dots used in this work
This section describes the synthesis of the QDs used in this work, the sample preparation scheme and an estimation of excitation saturation.

S2.1 Quantum dot synthesis and sample preparation
Colloidal CdSe/CdS/ZnS core/shell/shell QDs were synthesized by the following protocol: A cadmium oxide (CdO), n-tetradecylphosphonic acid (TDPA), and 1-octadecene (ODE) mixture was heated to 280°C in a three-neck flask under argon environment. Next, a stock solution of trioctylphosphine selenium (TOPSe) was rapidly injected. The temperature was then reduced to 250°C until the particles reached the desired diameter. A layer-by-layer growth technique in a one-pot synthesis method 3 was used for shell growth of cadmium sulphide (CdS) and zinc sulphide (ZnS). This resulted in spherical QDs with an outer diameter of 5.3 ± 0.6 nm (see Figure S3). Some of the QDs (< 10%) are slightly elongated up to an aspect ratio of 1:1.5. The quantum yield (QY) was measured with an absolute photoluminescence QY spectrometer (Quantaurus-QY, Hamamatsu), and is 90%.
Samples were prepared by spin coating a glass coverslip with a solution of these QDs dispersed in a 3wt% solution of poly(methylmetacrylate) (PMMA) in toluene.

S2.2 Quantum dot excitation saturation
To estimate QD excitation saturation, single QDs were illuminated at varying intensities. Figure S4a shows such an experiment, where over time the excitation intensity was increased at a steady rate from 28 W/cm 2 to 280 W/cm 2 and then back down. Higher illumination intensities typically result in more time spent in the 'grey' and 'off' states. Hence, to assess single-excitation saturation, it is essential to identify and estimate the peak occurring intensity of the 'on' state alone. This was achieved by creating an intensity histogram at each value of the excitation power, smoothing the histogram with a Gaussian filter and finding the smoothed histogram peak. The peak occurring values are shown in Figure S4b, together with a fit to a saturation model: 4 where P is the 'on' state peak, I excitation power, I sat saturation power and A asymptotic 'on' state peak (last two are the fit parameters). This simple model assumes negligible contribution to the intensity from multiexcitonic recombination, an assumption justified by the measured value of g (2) (0) ∼ 0.1 (see Figure S8). The data agrees well with the fit and I sat = 129 ± 12 W/cm 2 (68% confidence interval). The probability to excite n excitons following a single pulse can be estimated from the Poissonian distribution: using the distribution parameter λ = I Isat . At the excitation power used in this work (∼140 W/cm 2 , dashed purple line in Figure S4), the probability of exciting at least a single exciton following a laser pulse is ∼66% and of exciting at least twice to form a biexciton ∼30%.  Figure S4: Quantum dot excitation saturation. a) Saturation measurement. A single QD illuminated at increasing intensities from 28 W/cm 2 to 280 W/cm 2 and then back down (10 s, ∼25 W/cm 2 steps). Each point represents the detected intensity at a 5 ms time-bin, colored according to the local density of data-points for clarity. b) Peak occuring 'on' state intensity for each illumination power (blue circles) and a fit to a saturation curve (red solid line). The excitation power used in this work is marked by a purple dashed line. Figure S5 shows fluorescence decay curves for the different intensity states identified in Figure 3 of the main text. Each trace is a sum over all detector pixels for a specific intensity state. The 'on' and 'grey' state decay traces are both dominated by a single exponential term. However, the 'grey' state decay is an order of magnitude faster (τ on ≈ 20 ns and τ grey ≈ 1.35 ns, see fit details in next section). An irregular instrument response function (IRF, see subsection S1.2) results in an artifact visible at the first ∼4 ns of each curve.

S4 2D lifetime-spectrum analysis
Apart from heralded spectroscopy, presented in the main text, the single-particle, spectro-temporal information provided by the spectroSPAD can reveal connections between the spectral and the dynamical characteristics of nanocrystal fluorescence. One example of such an observation is presented in the 2D histogram shown  Figure S5: Fluorescence decay by intensity state.
Histogram of detection delays from the excitation pulse for the intensity states identified in Figure 2 of the main text. Note the different decay lifetime scales of the 'on' and the 'grey' state.
in the top panel of Figure S6: Photon detections (here only from the 'on' blinking state), are binned according to both their arrival time (with respect to excitation pulse) and their energy. In such a spectrum-lifetime dataset, one can differentiate the spectra of different lifetime components. While this type of data is commonly measured for an ensemble of particles (with a scanning monochromator), this is the first demonstration of such a measurement for a single QD.
In the following we attempt to distill the spectrum and lifetime of the biexciton (BX) state from the multicomponent data. We show that even when observing the 'on' state data alone, it is extremely difficult to separate the contributions of the charged and the BX state. As a result, we claim that determining the BX binding energy from such an analysis is a challenging task, and is prone to ambiguities.
To analyze the spectrum-lifetime data, we perform a global fit with the same two exponential decaying components for each of the energy-bins (detector pixels). That is, we fit every row in the matrix presented in the top panel of Figure S6, with a sum of two decaying exponential functions. As a result, for n pix pixels, we have 2 + 2·n pix fit parameters, one lifetime parameter and n pix amplitudes (denoting the spectrum) for each fit component. Note, that in order to reduce the effect of the irregular IRF (see subsection S1.2), we omit a 4.3 ns portion of the data in all fits (the apparent time gap in Figure S6).
Fitting the 'grey' state spectrum-lifetime data (not shown here), we observe that most of the contribution (> 80%) arises from a short lifetime component (1.35 ns) which we associate with the well-known short lifetime of the charged state. The respective spectrum of this state is shown in Figure S7 (purple dots) together with a fit to a Cauchy-Lorentz distribution centered at 1.99 eV (purple line). Next, we turn our focus back to the 'on' state spectrum-lifetime data shown in the top panel of Figure S6. Fitting the data with two decaying exponentials, we obtain a short (∼0.5 ns) and a long (∼20 ns) lifetime components. The long lifetime matches the common lifetime of the single exciton state (1X) in CdSe QDs and indeed its associated spectrum ( Figure S7, black dots) resembles that of the 'on' state spectrum presented in figures 3 and 4 of the main text. The peak of a fitted Cauchy-Lorentz distribution (black line) is at 2.007 eV, very close to 2.004 eV obtained from the fit of the heralded 'on' state spectrum.
Unlike the straightforward analysis of the long lifetime contribution, interpreting the short lifetime component is quite a challenging task. First, its integrated intensity is very small, ∼1/400 of that of the 1X term. In addition, at short time delays the residuals of the fit ( Figure S6, bottom panel), are substantial and systematic (switch from positive to negative sign around 0.4 ns), indicating that our model does not fully account for all features. Most importantly, the spectrum of the short lifetime component ( Figure S7, green dots) coincides with that of the charged state.
Considering the intensity proportion of the short lifetime component, the expected contribution to the data from BX recombination is given by where g (2) (0) is the second-order correlation function at zero delay time (antibunching factor) and p(N ≥ k) is the probability to populate k or more excitons after a single laser pulse, given by a Poisson distribution. In the case of the present measurement, the probability of exciting at least a single exciton was 66% and that of at least two excitons was 30% (see subsection S2.2). The expected calculated intensity proportion is more than an order of magnitude larger than that of the short lifetime fit component. It is therefore evident that the short lifetime fit component in the 'on' state doesn't match the expected contribution of the BX in both the spectral position and in amplitude. We could not pinpoint the reason for the lack of a fast decaying signal for the BX recombination in the spectrum-lifetime fit. However, we suspect that both the short lifetime of the BX decay and the systematic errors in the fit result in a lack of accuracy for such an analysis. For example, even slight differences in IRF between different pixels can result in a systematic bias that can be overcome only with a complex characterization and modeling procedure.
While a more sophisticated numerical analysis protocol (e.g. taking the IRF into account) may lead to better results, the above analysis shows the challenges of interpreting such a data set. Namely, the contribution of the BX state is extremely weak and its characteristics overlap other states. For example, a fit that includes three decaying exponentials for the 'on' state resulted in over-fitting of the results and the emergence of unnatural spectral features. In contrast, heralded spectroscopy, presented in Figure 4 of the main text, unambiguously isolates the BX component from the 1X and charged 'grey' state, enabling a reliable and robust estimation of the BX spectra.

S5 Biexciton quantum yield
The BX QY to 1X QY ratio can be estimated from the second order correlation of photon arrival times seen in Figure S8. The figure was generated and corrected by applying the method described in reference 5 to the spectroSPAD measurements. Briefly, the measurement described in the main text can be viewed as a multiple arm Hanbury Brown and Twiss correlation setup, where the combination of spectral thermal broadening and the spectrometer are used instead of beamsplitters to multiplex the detection. Second order correlation of arrival times is calculated by histograming delays between detections (arrival times) in different pixels, and temporal and intensity corrections are calculated from the same data used in section S6. The result are peaks separated by the pulse repetition period (200 ns), and broadened by the fluorescence decay lifetime (τ ∼20 ns). The zero delay peak is significantly lower indicating antibunching -the lower probability of two photon detections following the same pulse. The ratio of BX QY and 1X QY, is equal to the normalized second order correlation at zero time delay (g (2) (0)), calculated as the ratio between the zero delay peak area and the mean area of all the nonzero delay peaks. The value for this specific QD estimated from Figure S8 is g (2) (0) = 0.09 ± 0.02 (standard error). The artifact seen near zero time delay originates from shot noise on the significant crosstalk correction. A thorough study of identical QDs from the same synthesis with a more accurate setup showed that this value is within the expected range of g (2) (0) = 0.10±0.02 (standard deviation). 5 The 1X QY is ∼90% (see section S2), and hence the BX QY is ∼9%.
The ratio between heralded photon pair detections and all detections (α BX ), under the excitation conditions used in this work (see subsection S2.2), can be estimated as where p det is the probability to detect a photon per laser pulse from a single QD. The factor of 0.45 is the result of the ratio between of the probability of exciting at least a single exciton (66%) and that of at least two excitons (30%) (see subsection S2.2). During the 'on' state of the QDs measured in this work, values of p det ∼ 10 −2 and g (2) (0) ∼ 10 −1 are measured, resulting in α BX ∼ 2 · 10 −4 . This value is in excellent agreement with the α BX value extracted by the heralded approach ((2.0 ± 0.7) · 10 −4 for the 30 QDs shown in Figures 4c and 5b of the main text), further corroborating the reliability of these results.

S6 Analysis details
This section details the heralded spectroscopy parameters and the temporal and intensity corrections appended to the acquired data.

S6.1 Heralded spectroscopy
Photon pairs were time gated to support correct identification of BX and 1X emission and reduce the contribution of background due to dark counts. BX emission was gated to the first 5 ns following the excitation pulse. Due to the short lifetime of the BX state, this gating leads to a negligible loss of signal accompanied with a significant reduction of noise. 1X detections were gated to 5 − 60 ns delay from the BX detection. The upper bound serves to reduce noise while accommodating the longer lifetime of the 1X state emission (see section S4). The lower limit filters out possible misidentification of BX and 1X due to the IRF (see subsection S1.2). The upper limits for BX and 1X also assert that only photon pairs following the same pulse are taken into account. Emission spectral peaks were then estimated by a fit to a Cauchy-Lorentz distribution. In Figure 4 of the main text, E BX and E 1X correspond to the energies of the first and second photons, respectively, of the pairs identified by the time gating described above. In Figure 5a of the main text, the momentary mean 1X energy ( E 1X 1 ms ) and averaged 1X energy ( E 1X 10 s ) are estimated at the time of each such pair event (as described in Figure 3d of the main text). ∆E 1X E 1X 1 ms − E 1X 10 s is used as an estimator of the momentary spectral shift (horizontal axis). The momentary BX binding energy is estimated as the difference between E 1X 1 ms and E BX (vertical axis). Using these estimators (i.e. E 1X 1 ms instead of E 1X ) is beneficial when estimating momentary spectral shift and binding energy for each pair event separately, as the 1 ms averaging averages over the thermal broadening distribution of 1X emission, reducing its effect on the final results.

S6.2 Temporal Corrections
The time-to-digital converter (TDC) architecture assigns timestamps with a mean interval of ∼18 ps (the detector jitter is larger, see subsection S1.2). However, as detailed elsewhere, 1 the timestamps are not uniformly spaced but rather each span a 0 − 92 ps range of arrival times (most spans are within 18 ± 12 ps). This non-uniformity was characterized by illuminating the detector with temporally featureless halogen light, and recording the occurrence of each timestamp as a measure of the relative time duration it spans. The correction was then implemented statistically by assigning to each recorded raw-timestamp a corrected-timestamp chosen at random from the respective time span. In addition, timestamps recorded for each detector pixel are differently delayed from the TDC trigger. This was characterized by illuminating the detector directly with the <160 ps FWHM excitation laser pulse, and adding a per-pixel timestamp delay so that the recorded pulse peaks in all detectors temporally overlap.

S6.3 Intensity corrections
Two sources of false detections and detection pairs had to be considered in the analysis of heralded spectroscopy. The first, dark count rate (DCR), was recorded and subtracted from the intensity trace (per pixel). The expected number of DCR-photon detection pairs was estimated and subtracted from the photon pair signal (DCR-DCR pair occurrence is negligible). The second source of false photon pairs, detector crosstalk, was characterized and corrected statistically, by the protocol detailed in reference. 5 In Figure 5a of the main text, each data-point represents a single photon pair event. Hence, only in this figure, intensity corrections could not be implemented with this statistical approach. The number of DCR and crosstalk induced pairs in this figure can be estimated to be ∼10% of the overall data-points. Hence, the use of a median estimator (red crosses in Figure 5a of the main text) mitigates the effect of noise induced outliers. Furthermore, to avoid biases where noise might be more significant than signal, the median values shown and considered for the fit were only for ∆E 1X energy-bins including at least 1% of the total signal.
All corrections were verified to be stable over time.

S7 Published values of BX binding energy
In order to compare the measured value of BX binding energy in this work to previously reported values in II-VI nanoparticles we have summarized a survey of the literature in Table S1. We note that different works use opposite conventions with respect to the definition of sign in biexciton binding energy. Here, we adopt the convention in which a positive binding energy (ε b > 0) refers to a BX spectral peak at a lower energy with respect to the 1X spectral peak due to attractive interaction between the excitons. In the case of single particle techniques, such as the current work, we note only the ensemble average BX binding energy in the table. This work * uPL -PL upconversion, TRPL -Time resolved PL, TA -Transient absorption, CS -Single particle PL spectroscopy at cryogenic temperatures, TGPL -Transient grating PL, 2DFS -Two dimensional fluorescence spectroscopy, HS -Heralded spectroscopy.  Figure S7: The spectra of different lifetime components. Following the fit procedure described in section S4, we obtain a spectrum for each lifetime component separately. The black and green dots show the long and short lifetime spectral contributions to the 'on' intensity state data, respectively. The purple dots show the spectrum of the charged state obtained through the fit of the 'grey' intensity state spectrum-lifetime data. The matching color lines present Cauchy-Lorentz fits for each spectrum. For comparison with the heralded spectroscopy method, described in the main text, we present the fit of the 1X (BX) spectrum for the same QD with a blue (red) dashed line (as shown in Figure  4b of the main text