Transient measurement of phononic states with covariance-based stochastic spectroscopy

We present a novel approach to transient Raman spectroscopy, which combines stochastic probe pulses and a covariance-based detection to measure stimulated Raman signals in alpha-quartz. A coherent broadband pump is used to simultaneously impulsively excite a range of different phonon modes, and the phase, amplitude, and energy of each mode are independently recovered as a function of the pump–probe delay by a noisy-probe and covariance-based analysis. Our experimental results and the associated theoretical description demonstrate the feasibility of 2D-Raman experiments based on the stochastic-probe schemes, with new capabilities not available in equivalent mean-value-based 2D-Raman techniques. This work unlocks the gate for nonlinear spectroscopies to capitalize on the information hidden within the noise and overlooked by a mean-value analysis.

with the eective polarizability operatorα and classical electric elds i ∈ {1, 2}, of envelopes E i (t), carrier frequencies ω i0 , and associated Fourier transforms The pump E 1 (t + T ) and probe E 2 (t) pulses are centered around t 10 = −T and t 20 = 0, respectively.
The ladder diagrams and level scheme of the ISRS signal are depicted in Fig. S1. For nonoverlapping pulses, the system is rst excited by the pump pulse E 1 (t + T ) from an initial equilibrium stateρ (0) to a coherent superposition of vibrational states. To rst order in perturbation theory, the interaction with the pump pulse and the resulting evolution of the system's density matrixρ(t), of elements ρ nn ′ (t), are given by with the transition energies ω nn ′ = Ω(n − n ′ ) and decay rates γ. The two terms in the above sum, ρ n ′ n ′ and ρ (0) nn , are associated with the top and bottom diagrams in Fig. S1, respectively. For nonoverlapping pulses, we can extend the above integral to t → +∞, so that Eq. (S4) can be recast as ρ nn ′ (t) = ρ 0,nn ′ + i e −iω nn ′ (t+T ) e −γ(t+T ) The pump pulse used in the experiment is a coherent pulse of Gaussian envelope E 1 (t) = E 10 e −σ 2 1 t 2 /2 e −iω10t ,  Figure S1. Ladder diagrams and level schemes for the third-order ISRS signal.
The evolution of the superposition of vibrational states generated by the pump pulse is monitored by a subsequent stimulated Raman process induced by the probe pulse E 2 (t). The frequency-dispersed heterodyne-detected signal is dened as the time-integrated rate of change of photon number in the probe pulse [2] S(ω s , T ) = − dt dN (2) whereN ω =â † ωâω denotes the number operator for a photon mode with frequency ω. The superscript (2) inN (2) ωs stresses that the modes are restricted to those occupied by the probe pulse E 2 (t). By using the Heisenberg equations of motion for the photon number operator, the signal at the detected frequency ω s reads with the time-dependent polarizability ⟨α(t)⟩ = Tr{αρ(t)} = nn ′ α n ′ n ρ nn ′ (t).
(S12) For nonoverlapping pulses, by expressing the elements of the density matrix via Eq. (S8), the time-and frequencydependent ISRS signal can be recast as contains the amplitude and phase information of the vibrational coherence generated by the pump pulse. For α n,n+1 ≈ α, this coherence term A reduces to (S15) Additional features appear in the presence of overlapping pump and probe pulses, which are responsible for the hyperbola branches in the (ω s , T ) plane observed in the experiment. These features have been previously analyzed [3] and are not the focus of this work.

S2. MODELING OF THE STOCHASTIC PROBE PULSE
The stochastic probe pulse used in the experiment is obtained from an initial coherent pulse, here modeled as a Gaussian pulseg of bandwidth σ 2 = 6.6 THz (40 fs duration). The pulse spectral phase φ(ω) is shaped by a pulse shaper with nite resolution R = 0.1 THz, leading to the shaped stochastic pulsẽ with central frequency ω 20 = 375 THz (800 nm) and peak strength E 20 . In Sec. S5, a pulse with noise present only on part of its spectral width will be considered.
The spectral phase φ(ω) is obtained by generating an array of N independent random numbers φ k for a discrete set of frequencies ω k = ω 20 + kR, k ∈ {−N/2, −N/2 + 1, · · · , N/2 − 1}. Each number is randomly distributed in the interval [−φ,φ], with vanishing expectation value ⟨φ k ⟩ = 0 and with ⟨φ k φ k ′ ⟩ = δ kk ′φ 2 /3. This is convolved with a Gaussian smoothing function with width R, resulting in The expectation value of the spectral phase is thus equal to ⟨φ(ω)⟩ = 0 (S19) and its correlation function is given by with the variance In order to interpret the spectra obtained by the covariance approach used in the main text, we introduce the twoand four-point correlation functions of e iφ(ω) , dened as respectively. To calculate these eld correlation functions, we rst introduce the correlation function ⟨φ(ω)φ(ω ′ )⟩ of the rst derivative of the spectral phaseφ(ω) = dφ/ dω. The stochastic probe pulse employed here is a wide-sense stationary process, for which only depends on the dierence ω − ω ′ . We thus introduce the functioñ and use the second-order cumulant expansion to write the two-[Eq. (S22)] and four-point [Eq. (S23)] correlation functions in terms ofp(ω) [4], leading to the following identities: The two-point correlation function F 2,sh (ω, ω ′ ) of the spectral pulseẼ 2,sh (ω) can be calculated in terms ofF 2 (ω, ω ′ ) and is equal to Eq. (S17), for values of the variance ⟨φ 2 ⟩ = (π/4) 2 and ⟨φ 2 ⟩ = (π/2) 2 , respectively. We notice in particular that, due to the method used to generate the stochastic spectral phases, these do not uniformly cover the whole phase interval [−π, π]. This is evident when ⟨φ 2 ⟩ = (π/4) 2 , and is still the case even for ⟨φ 2 ⟩ = (π/2) 2 . As a result, the amplitude of the average stochastic pulse ⟨Ẽ 2,sh (ω)⟩ does not vanish, even though its peak value reduces for increasing values of ⟨φ 2 ⟩. In time domain, this results in the appearance of a central peak surrounded by a noisy tail: this central peak survives in the envelope of the average stochastic pulse ⟨E 2,sh (t)⟩, whereas the noisy tails average out to 0. The two panels in Fig. S4 show the associated correlation functions ⟨φ The curves display the behaviour of the correlation functions for ⟨φ 2 ⟩ = (π/4) 2 and ⟨φ 2 ⟩ = (π/2) 2 . The correlation function of the rst derivative of the spectral phase shows an oscillatory behavior, with a strong positive central peak surrounded by two negative peaks. For very close frequencies ω and ω ′ lying within the width of the central peak, the associated spectral phases will have correlated rst derivativesφ(ω) andφ(ω ′ ) with the same sign.
However, for increasing frequency dierences, the spectral phases do not become directly uncorrelated: when the frequencies ω and ω ′ lie within the width of the second negative peak, this means that the associated spectral phases are still correlated, with rst derivativesφ(ω) andφ(ω ′ ) exhibiting opposite sign. As a result, the autocorrelation of the eld ⟨Ẽ * 2,sh (ω)Ẽ 2,sh (ω ′ )⟩ features one strong central peak. Out of this peak, however, the stochastic pulse still presents a nonvanishing autocorrelation function. The amplitude of this residual autocorrelation decreases with increasing values of the variance ⟨φ 2 ⟩.
In the experiment, the stochastic probe pulse was obtained by shaping its spectral phase. Nevertheless, the measured spectral intensity |Ẽ 2 (ω)| 2 features chaotic oscillations as apparent in Fig. 1b in the main text. This is due to the nite resolution Q of the optical detector. The detected eld is given by the convolution of the stochastic probe pulsẽ E 2,sh (ω) generated by the pulse shaper and the function q Q (ω) modeling the nite resolution of the optical detector   Figure S3. Numerical simulation of the stochastic probe pulses in Eq. (S17) for ⟨φ 2 ⟩ = (π/2) 2 . The dierent panels show the same properties of the pulses as in Fig. S2.
The spectral gating function q Q (ω) acts over the fast random oscillations of e iφ , aecting both amplitude and phase ofẼ 2 (ω) and leading the the intensity uctuations depicted in Fig. 1b in the main text. The two-and four-point correlation functions of the detected eld are thus given by -20 -10 0 10 20 of the stochastic probe pulses in Eq. (S17) for (blue) ⟨φ 2 ⟩ = (π/4) 2 and (yellow) ⟨φ 2 ⟩ = (π/2) 2 . and where we moved the broadband Gaussian envelopesg(ω) out of the integrals. In the following, we will model as a rectangular function of width Q = 0.1 THz and unit area, and we will use Eqs. This diers from this work, where the spectral gating function is included to reproduce the nite resolution of the optical detector, and where we set Q = R. This determines the shape of the two-and four-point eld correlation functions, which will be discussed in the following.

S3. COVARIANCE-BASED SPECTROSCOPY AND PEARSON COEFFICIENTS
In this Supplementary Section, we dene the Pearson coecients analyzed in the main text. The experiment is repeated for independent realizations of a stochastic probe, and the spectra are quantied via the Pearson coecient where I 1 (ω s1 ) and I 2 (ω s2 ) can either be the detected input intensity of the probe pulse [Ẽ 2 (ω s ) in Eq. (S29)] before interaction with the sample or the detected output intensity of the probe pulse after transmitting through the sample I out (ω s , T ) = I in (ω s ) + S(ω s , T ).
In the denition of the Pearson coecient, σ i (ω si ) denotes the standard deviation of I i (ω i ) given by In the following, we will focus on the autocorrelation of the input intensity, and the cross-correlation of input and output intensities. In all above formulas, the signalsignal correlation terms have not been included since they are weaker than the correlation between input intensity and signal.
The cross-correlation P (I in (ω s1 , T ), I out (ω s2 , T )) corresponds to three-dimensional matrix ρ c (ω IN , ω OUT , T ) displayed in Fig. 1c in the main text, where ω s1 = ω IN and ω s2 = ω OUT . As shown in Eq. (S38), this cross-correlation contains two main contributions, owing to the autocorrelation of the input pulse and the cross-correlation between input intensity and ISRS signal, respectively. This latter term, is responsible for the o-diagonal lines appearing in the bottom-right quadrant of the correlation map displayed in Fig. 1c in the main text, and is therefore the term we will focus on in the following.

S4. AUTO-AND CROSS-CORRELATIONS FOR A STOCHASTIC PROBE WITH NOISE DISTRIBUTED IN THE ENTIRE SPECTRAL BANDWIDTH
In the main text, we presented results for the partially stochastic pulses of Fig. 1b, where the pulse is shaped only at frequencies larger than 370 THz. For modeling purposes, we here start by describing the case of a pulse with noise distributed in its entire spectral bandwidth. The case of a partially shaped pulse is presented in Sec. S5.
In this Supplementary Section, we therefore consider the stochastic probe pulse of Eq. (S17), with noise distributed on the whole frequency spectrum. We set ⟨φ 2 ⟩ = (π/2) 2 and include the eect of the nite resolution of the optical detector via Eq. (S29) as a rectangular function q Q (ω) [Eq. (S32)] of width Q = 0.1 THz. Dierent choices of the function modeling the optical detector do not modify our conclusions. We consider nonoverlapping pump and probe pulses, so that the signal S(ω s , T ) contributing to the output intensity I out (ω s , T ) is given by S ISRS (ω s , T ) in Eq. (S13).
The Pearson coecients for auto-and cross-correlation of input and output intensities are exemplied for a system with a single phonon of frequency Ω = 10 THz and lifetime 1/γ = 1.6 ps.
A simulation of the autocorrelation of the detected input intensity, given by the Pearson coecient P (I in (ω s1 ), I in (ω s2 )) in Eq. (S37), is shown in Fig. S5. This is calculated in terms of the two-[Eq. (S30)] and with the standard deviation (S41) Figure S5 exhibits a strong peak along the main diagonal ω s1 = ω s2 with width determined by the pulse-shaper and optical-detector resolutions. We notice also the presence of additional, lower-intensity peaks around the central one, as visible on the right panel of Fig. S5.
When input and output intensity are cross-correlated, the resulting signal carries the additional contribution due to P I,S (ω s1 , ω s2 , T ) from Eq. (S39). With the expression for the signal in Eq. (S13), this contribution reads where we have dened the functions The correlation function P I,S (ω s1 , ω s2 , T ) depends separately on ω s1 and ω s2 due to the presence of the Gaussian envelopesg 2 (ω − ω 20 ). However, for suciently broadband pulses, the key features of P I,S (ω s1 , ω s2 , T ) are determined byG(ω s2 − ω s1 , ±Ω), which depends on ω s1 and ω s2 only via the Raman frequency This is the equivalent of δω in the main text. The functionG(ω R , ±Ω) is responsible for the o-diagonal features observed in the experimental spectra of Fig. 1c in the main text. We calculate it numerically, including the eect of the nite resolution of the optical detector, and display it in the left panel of Fig. S6. The function exhibits a structure consisting of two main peaks,G centered on ω R = 0 and ω R = ±Ω, respectively. The shape of each peak,f (ω R ), is highlighted in the right panel of Fig. S6. It contains two narrow peaks, mostly determined by the correlation length R of the stochastic pulses used. We could check numerically that the function q Q (ω), modeling the optical detector, inuences the width of the two peaks and the depth of the central local minimum appearing inf (ω). However, this does not aect the central frequency of the peaks inG(ω R , ±Ω). Scanning the Raman frequency ω R allows one to identify the energy ±Ω of the phonons in the system.
The resulting function P I,S (ω s1 , ω s2 , T ) is displayed in Fig. S7. Two peaks along ω s2 − ω s1 = ±Ω are clearly visible in the left panel, due tof (ω R ∓ Ω) inG(ω R , ±Ω). In the right panel, we highlight the oscillations displayed by both peaks as a function of time delay. The oscillations along the two peaks take place at a frequency given by Ω and have opposite phase, owing to the π shift between the initial coherences, A and A * , generated by the pump pulse. The envelopeg(ω s − ω 20 ) of the probe pulse determines the dependence of the standard deviations σ i (ω s ) upon frequency. g(ω s − ω 20 ) and, thus, σ i (ω s ) become very small far from the central frequency ω 20 . The corresponding increase in the strength of the peaks is apparent in the left panel of Fig. S7.  Figure S6. (left)G(ωR, Ω) between the spectral function of the detected probe pulse evaluated at four dierent frequencies.

S5. AUTO-AND CROSS-CORRELATIONS FOR A STOCHASTIC PROBE WITH NOISE PARTIALLY DISTRIBUTED IN THE SPECTRAL BANDWIDTH
In this Supplementary Section, we consider a stochastic probe pulse with noise only partially distributed within the spectral width of the pulse, i.e., where we have introduced the Heaviside step function θ(x) and the threshold frequency ω th = 370 THz, below which no spectral shaping is introduced. This corresponds to the experimental pulse of Fig. 1b in the main text. We assume ⟨φ 2 ⟩ = (π/2) 2 and model the nite resolution of the optical detector via the rectangular function in Eq. (S32), with Q = 0.1 THz. We consider nonoverlapping pump and probe pulses, with the signal given by the ISRS spectrum S ISRS (ω s , T ) of Eq. (S13), and study the dynamics of one single phonon of frequency Ω = 10 THz and lifetime 1/γ = 1.6 ps. For ω s1 and ω s2 both greater than the threshold frequency, the same auto-and cross-correlations are expected as in Sec. S4. However, when at least one of these two frequencies is lower than ω th , a dierent behavior is observed. This is discussed in the following for the input-intensity autocorrelation and for the cross-correlation between input pulse and signal. Figure S8 exhibits the autocorrelation function of the detected input intensity [Eq. (S37)]. For ω s1 > ω th and ω s2 > ω th , the same peak along the main diagonal ω s1 = ω s2 is found as in Fig. S5. For frequencies lower then ω th , however, where the spectral phase is unshaped, the associated standard deviation σ in (ω s ) would vanish exactly in our model, and the associated Pearson coecient P (I in (ω s1 ), I in (ω s2 )) in Eq. (S37) would diverge. To avoid this divergence and obtain simulation results in agreement with the experiment, we assume here and in the following a small, nite standard deviation ϵ at frequencies ω < ω th where the pulse does not undergo any modulation of the spectral phase. This simple assumption allows us to obtain results in agreement with the experimental analysis. Figure S9 displays the cross-correlation P I,S (ω s1 , ω s2 , T ) between the input intensity of the probe pulse and the ISRS signal. In order to understand the behavior of the cross-correlation function, we rst observe that for ω s1 < ω th , where the input intensity does not feature any added noise, no correlation can exist between I in (ω s1 ) and the signal Figure S8. Pearson coecient P (Iin(ωs1), Iin(ωs2)) of the autocorrelation of the detected input intensity for a partially modulated stochastic probe pulse. Figure S9. Contribution of the correlation function between the detected input intensity of the probe pulse and the ISRS signal to the Pearson coecient P (Iin(ωs1), Iout(ωs2, T )) for a partially modulated stochastic probe pulse. The cross-correlation function is displayed (left) as a function of the two frequencies ωs1 and ωs2 for T = 500 fs, and (right) as a function of the Raman frequency ωR and time delay T for a mean frequency (ωs1 + ωs2)/2 = 374 THz. The colorbar in the left panel was chosen to highlight the weaker peaks appearing in the top-right quadrant.
S ISRS (ω s2 , T ), independent of the value of ω s2 . On the other hand, the signal is determined by the action of the stochastic probe pulse at frequencies ω s2 and ω s2 ± Ω. For cross-correlations to appear, it is necessary that also these frequencies lie above the threshold frequency ω th . In Fig. S9, three peaks can be distinguished. In contrast to Fig. S7, however, these peaks only appear when both contributing frequencies lie above the threshold frequency: (i) the central peak at ω s1 = ω s2 is only present when both frequencies are larger than ω th ; (ii) the o-diagonal peak along ω s1 = ω s2 − Ω only appears if ω s1 > ω th and ω s2 − Ω > ω th ; and (iii) the o-diagonal peak along ω s1 = ω s2 + Ω only appears for ω s1 > ω th and ω s2 +Ω > ω th . This creates a clear asymmetry in the cross-correlation function with respect to the Raman frequency ω R , as visible in Fig. S9. Also in this case, we have assumed a small nonvanishing variance of the input and output intensities also for frequencies below ω th , in order to avoid a vanishing denominator and divergences. This agrees with the experimental results, explaining the large strength of the peaks in the bottom-right quadrant of the left panel of Fig. S9 compared to the weaker peaks in the top-right quadrant.

S6. COMPARISON WITH TRADITIONAL PUMP-PROBE SPECTROSCOPY
We have performed mean-value based pump-probe measurements on the same alpha-quartz sample as a comparison to the covariance-based pump-probe spectroscopy technique. Table S1 summarizes the phonon dephasing times extracted for the three most prominent phonon modes within the pulse bandwidth used in these measurements. The mean-value and covariance-based measurements were performed one after the other, so experimental parameters could be kept consistent. The dephasing times are consistent within error margins between the two measurement techniques.  Table S1. Phonon dephasing times measured using covariance based and mean-value based pump-probe techniques. Figure S10. Standard deviation of a 15x15 pixel section of the covariance map for ∆T = −3 ps in the quadrant where signal normally appears. This is a measure of the background due to spectrally uncorrelated noise in the detection system.

S7. UNCORRELATED BACKGROUND NOISE IN 2D CORRELATION MAPS
Spectrally uncorrelated noise introduced by the detection system (e.g. detector noise and electrical readout noise) produces a background contribution in the eventual 2D covariance maps that can obscure weak signals from the sample.
However, this undesirable contribution to the measured spectra is relatively small compared to the real uctuations in the probe pulse which contain the sample-induced correlations. In the correlation maps, the uncorrelated detection noise presents as a random signal nearly independent of ω OU T and ω IN , with positive and negative values distributed normally around zero. We can thus characterize the background noise level with a standard deviation of the values in areas of the covariance map where there is not signal present. This standard deviation of values decreases for increasing number of unique spectra acquired.
To quantify this eect in our measurement, we have taken the standard deviation (σ ρ ) of all points in a 15x15 pixel box within the lower right quadrant of the cross-correlation map (as indicated by the black dashed box in g. 1c in the main text) for a single pump-probe delay. A pump-probe delay of −3 ps was chosen to ensure no signal was present. In Fig. S10 σ ρ is plotted as a function of the number of unique spectra included in the calculation (N ), showing a rapid decrease in σ ρ below N = 800, followed by a slower decrease thereafter, tangentially approaching σ ρ = 0.01. A thorough analysis of this noise is beyond the scope of this work, but it demonstrates that the signal (which has peak values of ρ = ±0.1 − 0.2) can be clearly resolved even with relatively low N . A related question is how ρ can be connected to physically relevant quantities measured in traditional mean-value based experiments, such as ∆T /T . Unfortunately there are intrinsic dierences between the measurement techniques that hinder a direct connection. For example, ρ depends on two frequencies, whereas ∆T /T is measured at only a single frequency. Furthermore, the degree of correlation induced by a given signal depends on experimental parameters such as the amplitude and spectral width of the uctuations introduced in the probe, and over what fraction of the probe bandwidth they are applied.
With that said, we have attempted to make a comparison between the results reported here with mean-value ∆T /T measurements performed on the same sample in similar conditions. In this case, we nd that a correlation value of ρ = 0.1 is comparable to ∆T /T ≈ 0.01, but stress that this is not a universal relationship, and future work should be undertaken to understand in what circumstances ρ can be linked to ∆T /T or ∆R/R in a quantitative way.