LIFETIME-RESOLVED PHOTON-CORRELATION FOURIER SPECTROSCOPY

The excited state population of single solid-state emitters is subjected to energy fluctuations around the equilibrium driven by the bath and relaxation through the emission of phonons or photons. Simultaneous measurement of the associated spectral dynamics requires a technique with a high spectral and temporal resolution with an additionally high temporal dynamic range. We propose a pulsed excitation-laser analog of Photon-Correlation Fourier Spectroscopy (PCFS), which extracts the lineshape and spectral diffusion dynamics along the emission lifetime trajectory of the emitter, effectively discriminating spectral dynamics from relaxation and bath fluctuations. This lifetimeresolved PCFS correlates photon-pairs at the output arm of a Michelson interferometer in both their time-delay between laser-excitation and photon-detection T and the time-delay between two photons τ . We propose the utility of the technique for systems with changing relative contributions to the emission from multiple states, for example, quantum emitters exhibiting phonon-mediated exchange between different fine-structure states.


Introduction
The spectral dynamics of single emitters can broadly be categorized as fluctuations and relaxation. Fluctuations are temporally stochastic variations around the equilibrium configuration of any chemical system interacting with its environment at non-zero temperature. Relaxation refers to the system's return to the ground state equilibrium after preparation of a non-equilibrium state, for example, through a laser-driven creation of an excited state population. For single optical emitters, fluctuations can manifest as spectral diffusion, which is spectral jumping occurring from nanoseconds to seconds that is reflective of the microscopic interaction of the bath with the excited state. [1,2] Relaxation occurs via irreversible phonon-or spin-mediated dissipation or spontaneous emission of photons, processes typically observed from picoseconds to microseconds. [3] Simultaneous measurement of the full relaxation and fluctuation dynamics for single emitters requires a technique with high spectral and temporal resolution with an additionally high temporal dynamic range from picoseconds to seconds. No single such technique exists and different approaches present their own strengths and weaknesses. Streak-cameras can resolve the spectral evolution along the photoluminescence lifetime trajectory with picosecond resolution, but fail to resolve any dynamics beyond a few nanoseconds and typically lack single-molecule sensitivity. [4] CCD-based single-molecule emission spectroscopy can have millisecond temporal resolution, especially if auto-correlation of individual spectra recorded at high sweep-rates is performed, but fails to resolve even faster fluctuations and does not discriminate spectral dynamics along the lifetime trajectory. [5,6] Hong-Ou-Mandel (HOM) spectroscopy has been used to measure the coherence of single photons through two-photon interference. HOM can resolve energy fluctuations through a decrease in photon-coalescence efficiency on fast (picosecond to nanosecond) timescales and lifetime-resolved photon-presorting is straightforward. [7,8] However, HOM spectroscopy is exclusively suitable for single-photon emitters at low temperatures as it requires photon-coherences near the transform limit. Moreover, the dynamic range of HOM is practically limited to a few nanoseconds delay-time between the two interfering photons, rendering HOM of limited utility for measuring fluctuations occurring over many orders of magnitude in time. Cross-arXiv:2102.03706v1 [quant-ph] 7 Feb 2021 correlating spectrally-filtered photons provides higher temporal dynamic range, but is limited by the finite bandwidth of optical filters restricting the technique to broad lineshapes and spectral diffusion with large energetic spread. [9] Fourier spectroscopy is not bound by limitations in spectral bandwidth or two-photon delay-times as each photon selfinterferes. As a result, Fourier spectroscopy can readily be married with photon correlation to provide spectral readout with high temporal dynamic range and with arbitrarily high temporal resolution only limited by photon shot-noise. [10] This technique, Photon Correlation Fourier Spectroscopy (PCFS), has now been established as a powerful tool for the study of optical dephasing and spectral fluctuations of single emitters at low and room temperatures. [11,12,13]Despite the success in characterizing spectral fluctuations, so far PCFS was not able to resolve any spectral changes associated with relaxation of the system back to the ground state, for example phonon-mediated relaxation or energy transfer between different states. Here, we propose a pulsed excitation-laser analog of PCFS that readily extracts relaxation and fluctuation dynamics from single emitters. The proposed technique is shown in Figure 1. In conventional PCFS, all photons emitted after the continuous-wave laser excitation are used to compile the spectral correlation (a), the auto-correlation of the spectrum compiled from photon pairs with a temporal separation of τ along the macrotime axis t of the experiment and an energy difference of ζ. In lifetime-resolved PCFS, photon-pairs are additionally correlated in the microtime T after pulsed laser excitation. Specifically, the photons are binned according to their microtime, sometimes referred to as the TCSPC channel, and spectral correlations are calculated using these microtime-separated photons (b). The technique can be readily implemented using picosecond single-photon counting equipment as shown in Fig. 1c and requires high-throughput post-processing photon-correlation analysis. We show through numerical simulation that this lifetime-resolved PCFS technique can separate the lineshapes and spectral diffusion dynamics of systems with more than one emissive state as long as the relative weights of the emission from different states changes over the course of the photoluminescence lifetime. More broadly, lifetime-resolved PCFS can in principle extract spectral fluctuations and relaxation with temporal resolutions only limited by the IRF response time of single-photon counting modules.

Lifetime-resolved PCFS
In PCFS, the interferometer path-length difference is adjusted to discrete positions inside the coherence length of the emission and periodically dithered on second timescales and over a multiple of the emitter center wavelength. [10] The dither introduces anti-correlations in the intensity cross-correlation functions of the output arms that encode the degree of spectral coherence at a given center position. The lineshape dynamics is thus encoded in the intensity correlations as a function of the time-separation between photons τ . We show in the Supplementary Information that the PCFS equations can straightforwardly be expanded to include spectral dynamics along the microtime T . The central observable in lifetime-resolved PCFS for a spectrum s(ω, t, T ) dependent on the microtime t and macrotime T is given by the spectral correlation p(ζ, τ, T ) as where ... represents the time average. Equation 1 describes the central observable in lifetime-resolved PCFS and can be understood intuitively as a histogram of photon-pairs with a shared microtime T , a macrotime separation of τ , and an energy separation ζ. The form and interpretation of p(ζ, τ, T ) depend on the dynamics of the emissive system and will be discussed in the following sections. We first discuss the general form of the spectral correlation for a system undergoing spectral diffusion in section 2.2. We then discuss two universal systems that map onto many specific real-world scenarios. First, we consider a system of two uncoupled and lifetime-distinct radiating dipoles undergoing uncorrelated spectral diffusion in section 2.3. Second, we consider a system of two coupled radiating dipoles subject to population exchange and correlated spectral diffusion in section 2.4.

The effect of spectral fluctuations on the spectral correlation
We consider equation 14 for spectral fluctuations δω(t, T ) that present along the macrotime axis of the experiment around the center frequency ω 0 of a spectrum. We can write for the spectrum s(ω, t, T ) = s(ω, T ) ⊗ δ(ω − δω(t, T )), where ⊗ is the convolution, s(ω, T ) the undiffused spectrum, and δω(t, T ) the time-dependent shift from the center wavelength. Spectral fluctuations can be characterized by the correlation function C(τ ) = δω(t, T )δω(t + τ, T ) . The canonical form of any spectral correlation can then be recast as reflecting the transition from the undiffused spectral correlation (absent any fluctuations p(ζ, τ → 0, T )) to the diffused spectral correlation p(ζ, τ → ∞, T ) with the evolution of C(τ ). Note that for τ → 0, δω(t 1 , T ) = δω(t 2 , T ) and the spectral correlation thus reduces to the homogeneous spectral correlation p(ζ, τ → 0, T ) = ∞ −∞ s(ω, T )s(ω − ζ, T )dω .

A lifetime-distinct doublet undergoing Gaussian spectral fluctuations
We discuss a system of two uncoupled and lifetime-distinct radiating dipoles undergoing uncorrelated spectral diffusion and involving states |A and |B . The system's energy diagram is shown in figure 2a (inset). The microscopic interpretation involves a system with two emissive states coupled to different bath fluctuations. We show how lifetimeresolved PCFS can separate the homogeneous lineshape and spectral diffusion parameters of the the two transitions. The different emission lifetimes result in microtime-dependent relative weights of emission intensity originating from states |A and |B after equal populations have been prepared through laser excitation. We decompose the overall dynamic spectrum of the system s(ω, t, T ) into microtime-dependent components as s(ω, t, T ) = a(T )s A (ω, t) + b(T )s B (ω, t), where a(T ) and b(T ) are the relative probabilities of a given photon originating from either state |A or |B , and show that the spectral correlation expands as (see Supplementary Information). The terms quadratic in a(T ) and b(T ) represent the spectral auto-correlations of the individual states p AA and p BB , while the cross-terms involving p AB represent the cross-correlation of the spectra s A (ω, t, T ) and s B (ω, t, T ). The form of the spectral correlation can be understood intuitively because the spectral correlation is compiled from pairs of photons with origins drawn from the four possible combinations of |A and |B . Importantly, the left-and right-sided correlations p AB and p BA are not identical unless s A (ω) and s B (ω) share the same center frequency ω 0 and are symmetric in ω. Spectral diffusion is a ubiquitous process observed for many single emitters. Common descriptions of single-emitter spectral diffusion are the non-Markovian and discrete Poissonian Wiener process [14] or the mean-reverting Ornstein-Uhlenbeck process. [15] These processes describe spectral diffusion phenomenologically and for simplicity we consider a simple non-Markovian Poissonian Gaussian jumping model (GJM). [12] The GJM process is characterized by a time-invariant probability density for discrete spectral jump occurrence to a new spectral position drawn from a Gaussian probability distribution function over ω. For the two states |A and |B as denoted in the subscripts, we A,B for the probability of a given spectral shift at a point in time. Here, we have introduced the spectral fluctuation term δω A,B from earlier. The microscopic interpretation of this process is the time-stochastic variation of the bath assuming discrete conformations coupling to the system. The corresponding fluctuation correlation function can be written as C(τ ) = e −τ /τc and is described by an exponential decay with a characteristic spectral jump time of τ c . When the two states diffuse independently of each other, no correlation is present and C AB (τ ) = 0. In this case, δω A (t)δω B (t + τ ) = δω a (t) δω B (t + τ ) = 0 because independently diffusing emissive states will not be correlated and the cross-terms p AB and p BA in equation 3 only reflect the crosscorrelations of the inhomogeneous components p AB/BA (ζ, τ → ∞). Absent any memory of spectral fluctuations even at early τ , the time average over the spectral-correlations of all random configurations is the cross-correlation of the inhomogeneously broadened (diffused) spectra where σ A and σ B are the widths of the Gaussian probability envelopes of the diffused distributions of states |A and |B . We numerically simulate the system of independently-diffusing optical transitions with parameters commensurate with typical experimental cryogenic single-molecule spectroscopy (see Supplementary Information). The time-domain results of the simulation are discussed in Figure 2. The configuration of the system is shown in Figure 2(a). The corresponding lifetime exhibits biexponential decay behavior as expected. In (b) and (c), we compare the cross-correlation functions for two different slices with microtime ranges of T = 0 − 100ps and T = 2000 − 7000ps, where |A and |B are the dominant emissive states, respectively. Unlike for the static doublet discussed in the Supplementary Information, the cross-correlations g (2) X (τ ) indicate spectral dynamics evident from the loss of anti-correlation at longer τ . As we specify different jumping rates for the two states, the decay of the spectral coherence evident in (b) and (c) occurs at different τ . The PCFS interferogram derived from the cross-correlations (see Supplementary Information for the derivation) for photons emitted with a time constant of < 100ps is shown in (d) and informs on the loss of photon-coherence between 1µs and 1ms owing to the energy fluctuations of the photons emitted < 100 ps after laser excitation. In Figure 3 we discuss the same simulation results in the spectral domain. In (a), we show the full-width-at-halfmaximum (FWHM) of the spectral correlation for both T and τ , a representation that makes immediately obvious the differences in the homogeneous linewidths at early τ and the differences in spectrally-diffused linewidths at late τ . For completeness we also show p(ζ, T ) (d) and p(ζ, τ ) (b) for fixed τ and T , respectively. These two representations inform on the spectral evolution owing to spectral diffusion and changing relative emission contributions from different states, respectively. (c) displays the evolution of the spectral correlation from the narrow homogeneous spectrum with a Lorentzian lineshape to the diffused Gaussian lineshape. One capability of lifetime-resolved PCFS is the ability to extract the homogeneous linewidths of different lifetimedistinct states in the presence of fast spectral diffusion. We demonstrate this ability through a global fit to the T-dependent spectral correlation. We define a model for the fit as a linear combination of two Lorentzians and a Gaussian with floating linewidths parameters. The relative amplitudes p AA ,p BB , and p AB,BA are calculated according to 3 taking the weights a(T ) and b(T ) from fits to the emission lifetime into account. p AA ,p BB , and p AB,BA are also displayed in (d).
We apply a global fit to the slices of the spectral correlation p(ζ, τ = 60µs, T ) along T as shown in (e),(f) and (g). The cross-correlation p AB,BA present as a broad Gaussian background superimposed with the homogeneous Lorentzian spectral correlations p AA,BB as introduced in equation 3. The width of this Gaussian component is σ AB ≈ σ 2 A + σ 2 B . The homogeneous lineshape parameters parsed into the numerical model are extracted by the fit within photon shotnoise, thus validating the approach adapted herein. We note that in PCFS, the high temporal resolution achieved through photon-correlation comes at the cost of the loss of the absolute phase of the spectral information. In other words, both the asymmetry of the lineshape and the center frequency of s(ω) is lost in the spectral correlation p(ζ). The unambiguous reconstruction of s(ω) from p(ζ) is therefore impossible and the spectral correlation is typically fit with a model parametrizing a suitable form for the underlying emission spectrum, as we adapted herein. [13,12]

A dynamic doublet with population transfer and spectral fluctuations
We now turn to a system of two coupled radiating dipoles undergoing population exchange and subject to correlated spectral diffusion. A specific example would be solid-state quantum emitters undergoing incoherent and phononmediated population transfer after non-resonant excitation. [16] In quantum emitters, disentangling the relaxation rate and coherence times of the different fine-structure states in the presence of spectral diffusion is important for a detailed understanding of the dephasing process as phonon-mediated population exchange constitutes an important dephasing process in the solid-state. [3] We depict the system's energy diagram in Figure 4(a), which exhibits two excited states with equal oscillator strengths and an irreversible relaxation rate k from the higher to the lower-lying state. In this system, photon emission from the higher-lying state |A will start immediately after population of the state. Emission of the lower-lying state |B requires further relaxation and is often phonon-mediated. [3] The relative population of states |A and |B will thus change during the emission lifetime of the overall system as long as the relaxation rate k is faster than the radiative rate 1/T 1 of both |A an |B . The population dynamics of the system can be described by the following set of coupled equations: with the solutions: We show the effect of the changing relative cross-correlation probabilities between states |A and |B (p AA ,p BB ) in Figure 4 (b). Despite the T -invariant exponential population decay constant leading to a monoexponential photoluminescence lifetime of the overall system, the relative weights of p AA and p BB are changing with T .We show the spectral correlation of the lifetime-resolved PCFS experiment with indiscriminate T in (c). On timescales shorter than the spectral diffusion time τ , the fine-structure states are well-separated. At late τ , the broad diffused lineshape obfuscates the fine-structure splitting. We demonstrate that lifetime-resolved PCFS can recover the lineshape parameters of the homogeneous doublet by applying a least-squares fit of a suitable model to the T-dependent spectral correlation as shown in (d),(e),(f). The model consists of two Lorentzians. with the floating linewidths Γ 1 , Γ 2 , an energy offset Ω and a relaxation rate k, which determines the temporal change of the relative emission contributions of A and B . We recover all model parameters within photon shot-noise thus validating the utility of lifetime-resolved PCFS to extract the coherences and relaxation rates of different emissive fine-structure states. We note that the observation of early-τ multiplets in the spectral correlation compared to the broad Gaussian background in section 2.3 is the signature of correlated spectral diffusion dynamics between the two states. Our simulations suggest that measuring the photon-coherences of quantum emitters exhibiting spectral fluctuations and different emissive fine-structure states will provide an avenue to study quantum emitter optical dephasing through both fluctuations and population exchange between different electronic states.

Conclusions
We propose a new photon-correlation spectroscopic technique that extracts spectral fluctuations along the lifetimetrajectory of single emitters. The technique works through time-correlation of photons detected at the output arms of a variable path-length difference interferometer in both the microtime and macrotime domain and can be implemented using standard picosecond photon-counting electronics. We show that lineshape and fluctuation parameters can be extracted from the fits to the lifetime-resolved spectral correlations. Our technique opens up multiple frontiers in singleemitter spectroscopy. We emphasize that our technique is general, but point to its special utility in quantum emitter research enabled by the high spectral resolution required to resolve photon-coherences at low temperatures. Experimental efforts will be directed towards probing the fluctuation dynamics of non-stationary systems and investigation of the decoherence processes in quantum emitters. Specific materials are readily available such as emissive defects in diamond and emerging 2D materials as well as semiconductor nanostructures.

Acknowledgements
The lead author of this study (H.U., study conception, derivation, modeling and interpretation) was initially funded by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering (award no. DE-FG02-07ER46454) and funded by Samsung Inc. (SAIT) during the completion of the study. We thank Weiwei Sun, David Berkinsky, Alex Kaplan, Andrew Proppe, and Matthias Ginterseder for critically reading the manuscript and their feedback.
X (τ ) for different optical path-length differences δ 0 and microtimes T , where |A and |B are the dominant emissive states, respectively. The loss of coherence with increasing τ is evident from the reduction in anti-correlation. This coherence loss occurs at earlier τ for early-T photons (emission predominantly from |A , (b)) compared to late-T photons (emission predominantly from |B , (c)). The PCFS interferogram G (2) (δ, τ ) for early-T photons is shown in (d) and reflects the evolution from the exponential homogeneous dephasing at early τ to the spectrally-diffused Gaussian dephasing at late τ . We show the evolution of the weights of auto-and cross-correlations between states along T in (d).The weights are derived from the relative amplitude of the two exponential components of the photoluminescence decay in (Fig.2a). Taking p AA ,p BB , and p AB,BA into account, we apply a global fit to the spectral correlation along T to recover the lineshape parameters of the undiffused system as shown in (e),(f) and (g). The broad underlying Gaussian component in (f) reflects the cross-correlation of the diffused distributions of |A and |B and has a width of σ ≈ σ 2 A + σ 2 B . Figure 4: Lifetime-resolved PCFS simulation of two coupled dipoles undergoing population transfer and interacting with the same bath resulting in collective spectral diffusion of the doublet (a). We introduce a phonon-mediated relaxation rate between the upper and lower state of k relax = 1/80ps −1 . As the radiative rates of the two states are chosen to be equal, the emission lifetime follows a monoexponential decay behavior despite changing relative populations of |A and |B with the microtime (b). The spectral correlation irrespective for all photons irrespective of their microtime is shown in (c) and demonstrates the transition from a triplet at early τ to the spectrally-diffused distribution at late τ . The fine-structure splitting Ω, the linewidths Γ A,B , and the relaxation rate k can be recovered through lifetime-resolved PCFS and a global fit of the slices along T with a fixed macrotime correlation of τ = 8µs (d),(e) and (f).

Derivation of Lifetime-Resolved PCFS
We derive PCFS classically, that is using the field rather than the photon description of light. In practical PCFS experiments, the finite integration time renders the photon anti-bunching timescales outside of reasonable signal-tonoise regimes, justifying our classical approach adapted herein. [SI1,SI2] We denote the microtime as T , the macrotime of the absolute experimental clock as t, and the time-separation between photons as τ . We show how the experiment can access the spectral correlation p(ζ, τ, T ), a histogram of photon pairs with a temporal separation τ along the macrotime axis, a separation from the laser excitation pulse T , and an energy separation of ζ.
The intensity distribution at the two interferometer outputs a and b for an arbitrary spectrum s(ω, t, T ) evolving over the microtime T after laser excitation and the macrotime t of the absolute experimental clock at a given center path-length difference δ 0 is given by where δ(t) represents a time-dependent path-length difference defined by the chosen scanning trajectory of the interferometer and exhibiting the form δ(t) = δ 0 + ∆δ(t), where δ 0 is the center path-length difference and ∆δ(t) is a time-dependent term introduced experimentally through periodic dithering of one of the interferometer arms. F T {s(ω, t, T )} δ(t) is the Fourier transform of the time-dependent single-emitter spectrum. Substitution of 9 into the canonical form for the intensity cross-correlation function g where ...
The exact functional form of ∆δ(t) can be chosen arbitrarily as long as it is experimentally ensured that the pathlength difference for each recorded correlation function is scanned symmetrically around the center position δ 0 : Tac 0 ∆δ(t)dt = 0, where T ac is the data acquisition time at each δ 0 . Under this condition, terms linear in the Fourier transforms of the form F T {s(ω, t, T )} δ(t) or F T {s(ω, t + τ, T )} δ(t+τ ) vanish and the expansion of equation 11 can be recast as g (2) × (δ 0 , τ, T ) = The first term is just the intensity auto-correlation function g || (τ, T ). We then introduce the PCFS interferogram G (2) (δ 0 , τ, T ) defined as: where we have used the spectral correlation p(ζ, τ, T ), the auto-correlation of the spectrum compiled from photon pairs with a separation τ along the macrotime axis and a separation T from the laser excitation pulse. We note that equation 13 is only valid for δ(t) ≈ δ(t + τ ), i.e. for inter-photon lag-times much shorter than the dither period. Explicitly, for a microtime-and macrotime-variant spectrum s(ω, t, T ), the spectral correlation p(ζ, τ, T ) can be written as Equation 14 describes the central observable in lifetime-resolved PCFS and can be understood intuitively as a histogram of photon-pairs with a microtime T , a macrotime separation of τ , and an energy separation ζ. The form and interpretation of p(ζ, τ, T ) are dependent on the dynamics of the emissive system and are discussed in the manuscript.
The lifetime-distinct spectral correlation One instructive system conducive to investigation with lifetime-resolved PCFS is a system with two transition dipoles radiating with distinct oscillator strengths for the states |A and |B . The different radiative lifetimes result in microtimedependent relative weights of emission intensity originating from states |A and |B under δ-like non-resonant excitation. Substituting s(ω, t, T ) as s(ω, t, where a(T ) and b(T ) are the relative probabilities of a given photon originating from either state |A or |B , and separating into the sum of time-averaged terms .... yields for the spectral correlation. We can write more compactly: The terms quadratic in a(T ) and b(T ) represent the spectral auto-correlations of the individual states p AA and p BB , while the cross-terms involving p AB represent the cross-correlation of the spectra s A (ω, t, T ) and s B (ω, t, T ). The form of the spectral correlation can be understood intuitively because the spectral correlation is compiled from pairs of photons with origins drawn from the four possible combinations of |A and |B . Importantly, the left-and right-sided correlations p AB and p BA are not identical unless s A (ω) and s B (ω) share the same center frequency ω 0 and are symmetric in ω.

A lifetime-distinct static doublet
We now consider the expected form of the spectral correlation for a static doublet, where s A,B (ω, t 1 ) = s A,B (ω, t 2 ) for all t, in other words time-invariant spectral lines for each state |A and |B . Per definition, C(τ ) = 0 as δω(t) = 0 for all t. The corresponding terms of the spectral correlation are invariant in τ and equation 3 reduces to reflecting the absence of any dynamic spectral changes.
We numerically simulate the experiment (see Methods section) to showcase the technique and resulting spectral correlation for the static doublet as shown in Figure 1. We introduce two states with 100 and 1000 ps radiative lifetimes T 1 and different coherence times T 2 of 30 and 60ps with exponential decoherence (Lorentzian lineshape).
The corresponding fluorescence lifetime as compiled from the numerical photon-arrival time data is shown in (a) and reflects the biexponential decay corresponding to the emission from the two states. (b) shows the cross-correlation function g (2) X (τ ) for all photons irrespective of the microtime T . The individual correlation functions are τ -invariant as expected in the absence of spectral fluctuations. The degree of anti-correlation increases towards a maximum for δ = 0. The spectral correlation p(ζ, τ ) is shown in the contourplot in (c). The higher noise-level at early τ in (b) and (c) can be attributed to photon-shot noise. The lower panel in (c) shows the lifetime-resolved PCFS data for τ = 60µs. The evolution from the narrow (long coherence, fast lifetime) to broad (shorter coherence, slower lifetime) is consistent with the evolution of the emission probabilities from state |A and |B . We also plot the relative cross-correlation probability of photons originating from |A and |B . At early T , a significant fraction of photon-correlation counts stem from the cross-correlation between |A and |B as indicated by the side-peaks in the spectral correlation. These side-peaks are maximal around T ≈ 300ps, where the relative weights of emission from states |A and |B are equal. In PCFS, the high temporal resolution achieved through photon-correlation comes at the cost of the loss of the absolute phase of the spectral information. In other words, both the asymmetry of the lineshape (not not apply for the systems considered in this work) and the center frequency of s(ω) is lost in the spectral correlation p(ζ). The unambiguous reconstruction of s(ω) from p(ζ) is therefore impossible and the spectral correlation is typically fit with a model parametrizing a suitable form for the underlying emission spectrum. [SI3,SI2] Using a similar approach for lifetimeresolved PCFS, the lineshape parameters for the two different states can be extracted by applying a global fit to the spectral correlation of a spectral doublet with relative weights of the emission lines of a(T ) and b(T ). For the lifetime-distinct doublet, a(T ) and b(T ) can be extracted as the relative weights of the exponential fits describing the biexponential decay dynamics as shown in Figure 1 (c). We validate this approach by applying a global fit to the slices of the spectral correlation p(ζ, τ = 60µs, T ) along T using equation 11 and assigning a Lorentzian lineshape to |A and |B as shown in Figure 2e,f and g.

Numerical Simulation Methods
We use a custom MATLAB library to numerically simulate the relaxation-resolved PCFS experiments. All code is made publicly available on github. [SI4] The numerical simulations have three broad components. i) the simulation of the photon-emission of a given system, ii) the transcription of these photon-streams into time-tagged single-photon data and iii) the PCFS analysis of the time-tagged single-photon data. We first model the photon-streams with the PhotonSimulator.m class by writing the photon-emission stream into arrays containing the macrotime information t, microtime information T , center wavelength λ and a lineshape function g(t). For any given time increment of the experiment 1/f rep , where f rep is the laser repetition frequency, the probability of photon detection is drawn from a Poissonian distribution for reasonable total average detection count rates of around 10 5 photons/second and laser repetition rates of f rep = 10 − 20M Hz: p(λ, k) = λ k e −λ k! with k = 1 and λ = 1x10 5 cps/f rep . By omitting higher-order (multi-photon) detection events after the same laser excitation pulse with k ≥ 2, the photon-stream reflects the quantum behavior of the single emitter. The total time of photon emission simulated at each stage position was chosen to be 30-40 s, commensurate with a total integration time for the experiment of around 1-2 h. For the lifetime-distinct doublet, each photon is randomly assigned a state |A or |B and the corresponding center wavelength λ A,B and g(t) are chosen respectively. For the coupled dipole system, the population is drawn from the microtime-dependent population distribution that changes owing to the relaxation from state |A to |B . The microtime T is drawn from an exponential distribution with the respective emission lifetime T 1 A,B as p(T ) ∝ e −T /T 1 A,B . A single-mode photon exhibits a Lorentzian lineshape, and the corresponding lineshape function g(t) will present as exponential dephasing in the time domain. Our simulations are limited to Lorentzian lineshapes. The adaptation to non-pure single-photons i.e. photons in superposition states of multiple spectral modes, is straightforward through the choice of a Gaussian form for g(t). We introduce Gaussian spectral diffusion by drawing the probability for jump occurrence from the boolean set {0, 1} with a time-invariant probability for each laser pulse. For each jump, an energy offset δω is drawn from a Gaussian distribution and added to the color λ of the photon of the respective state. The action of the PCFS experiment is to assign photons of a given color, coherence time, and lineshape a detector output flag. This assignment is performed with the function PhotonAnalyer.m, which draws the exit flag from the boolean set ∈ {0, 1} with the respective probabilities p 0,1 given by: p 0,1 (δ(t), λ) g(x) = 1/2 ± g(δ(t)/c)cos( 2πδ(t) λ ).
The resulting photon arrival-time data is presorted according to the microtime interval of choice and the different microtime-sorted photon-streams are cross-correlated using a MATLAB implementation of the multi-tau algorithm to perform the PCFS analysis. [SI5] The details of the PCFS analysis are described in Utzat et al. [SI2].
X (τ ) at different stage positions transcribing the degree of spectral coherence into anti-correlations are shown in (b). The invariance of the correlation functions indicates the absence of any spectral fluctuations. The spectral correlation p(ζ, τ ) is shown in (c) for all T . A slice along T for a given τ is shown in (d). (e) shows the evolution of the weights p AA ,p BB and 2p AB/BA of the different contributions to the spectral correlations, highlighting the possibility of extracting the lineshape parameters from lifetime-resolved PCFS data. The corresponding interferogram is shown in (f) for three different T , where p AA , p AB,BA , and p BB are dominant, respectively. Taking the respective weights into account, the evolution of p(ζ, T ) with T can be fit with two Lorentzian peaks with widths Γ 1 , Γ 2 and their energy separation Ω, (g-i).