Probing the fluctuations of the optical properties in time-resolved spectroscopy

TL;DR

This paper introduces a statistical method in time-resolved spectroscopy to measure the dynamic noise of electronic currents in nonequilibrium materials, revealing insights into thermalization and fluctuation-dissipation relations.

Contribution

It proposes a novel framework for analyzing optical pulse fluctuations to extract dynamical information about electronic and other degrees of freedom in nonequilibrium states.

Findings

Allows measurement of current noise dynamics in ultrafast experiments

Reveals how electronic systems obey fluctuation-dissipation theorem during thermalization

Provides a general approach extendable to magnetic and vibrational excitations

Abstract

We show that, in optical pump-probe experiments on bulk samples, the statistical distribution of the intensity of ultrashort light pulses after the interaction with a nonequilibrium complex material can be used to measure the time-dependent noise of the current in the system. We illustrate the general arguments for a photo-excited Peierls material. The transient noise spectroscopy allows to measure to what extent electronic degrees of freedom dynamically obey the fluctuation-dissipation theorem, and how well they thermalize during the coherent lattice vibrations. The proposed statistical measurement developed here provides a new general framework to retrieve dynamical information on the excited distributions in nonequilibrium experiments which could be extended to other degrees of freedom of magnetic or vibrational origin.

Equations

$${\sigma }_{\text{bulk}}={ϵ}^2\int dtd{t}^{\prime }{E}_{+}(t)E_{-}(t^{\prime })⟨\delta j(t,\mathbf{R})\delta j(t^{\prime },\mathbf{R}){⟩}_{\text{ret}},$$

$$I_n=$$

$$\times ⟨T_{\bar{\tau }}[j(1)\cdots j(n)]T_{\tau }[j(1^{\prime })\cdots j(n^{\prime })]{⟩}_{\text{ret}}.$$

$${\sigma }_{\text{bulk}}\propto {ϵ}^2\int dtd{t}^{\prime }⟨\delta J(t_{\text{ret}})\delta J(t_{\text{ret}}^{\prime })⟩{\dot{E}}_{+}(t){\dot{E}}_{-}(t^{\prime }),$$

$$C(t_p,\omega )=2b(\omega ,\beta )\text{Im}{\chi }^R(t_p,\omega ),$$

$$H=-t_0\sum _{⟨ij⟩\sigma }{c}_{i\sigma }^{†}{c}_{j\sigma }+{\omega }_0\sum _j{b}_j^{†}{b}_j+g\sqrt{2}\sum _j{X}_j(n_j-1),$$

$$S=$$

$$-i\sum _{\sigma }{\int }_{C}d{t}_{1}d{t}_2{c}_{\sigma }^{†}(t_1)\Delta (t_1,t_2)c_{\sigma }(t_2)$$

$$F(t)=\sqrt{2}g\sum_{\sigma}\big{(}\langle c_{\sigma}^{\dagger}(t)c_{\sigma}(t)\rangle-0.5\big{)}.$$

$$h_{loc}(t)=-\sqrt{2}g⟨X(t)⟩.$$

$$\Sigma (t,t^{\prime })=i{g}^{2}G(t,t^{\prime })D(t,t^{\prime }),$$

$$D(t,t^{\prime })=-2i⟨T_C\tilde{X}(t)\tilde{X}(t^{\prime })⟩.$$

$$\begin{split}\big{(}i\partial_{t}&+\mu-h_{loc}\big{)}G(t,t^{\prime})-\\ &\big{(}\Delta(t,t^{\prime})+\Sigma(t,t^{\prime})\big{)}*G(t,t^{\prime})=\delta_{C}(t,t^{\prime}).\end{split}$$

$$P(t,t^{\prime })=-2i{g}^{2}G(t,t^{\prime })G(t^{\prime },t),$$

$$\big{(}1-D_{0}(t,t^{\prime})*P(t,t^{\prime})\big{)}*D(t,t^{\prime})=D_{0}(t,t^{\prime}).$$

$$D_0(t,t^{\prime })=$$

$$\displaystyle+\theta_{C}(t,t^{\prime})e^{-i\omega_{0}(t-t^{\prime})}\big{]},$$

$${\Delta }_a(t,t^{\prime })=v(t)G_b(t,t^{\prime })v(t^{\prime }),$$

$${\Delta }_b(t,t^{\prime })=v(t)G_a(t,t^{\prime })v(t^{\prime }),$$

$$\chi (t,t^{\prime })=i⟨T_{C}j(t)j(t^{\prime })⟩.$$

$$\delta ⟨j(t)⟩={\int }_{C}d\bar{t}\chi (t,\bar{t})\delta A(\bar{t}),$$

$${\hat{\psi }}_k=\left(\begin{array}{c}{c}_{k,a}\\ c_{k,b}\end{array}\right),$$

$$\displaystyle\langle j(t)\rangle=-2i\sum_{k\in RBZ}\text{tr}\big{[}v_{k-A}\hat{\sigma}_{1}\hat{G}_{k}(t_{+},t_{-})\big{]},$$

$$\displaystyle\chi(t,t^{\prime})=2i\sum_{k\in RBZ}v_{k}^{2}\,\text{tr}\big{[}\hat{\sigma}_{1}\hat{G}_{k}(t_{+},t^{\prime})\hat{\sigma}_{1}G_{k}(t^{\prime},t_{-})\big{]},$$

$$\sum _{k\in RBZ}f({ϵ}_k)={\int }_0^{\infty }\rho (ϵ)f(ϵ),$$

$$\sum _{k\in RBZ}{v}_k^{2}f({ϵ}_k)={\int }_0^{\infty }D(ϵ)f(ϵ),$$

$$\hat{h}(t)$$

$$\hat{\Sigma }(t,t^{\prime })$$

$$I_1=$$

$$I_2=$$

$$\times g(2,\bar{2})g(2^{\prime },{\bar{2}}^{\prime })⟨T_{\bar{\tau }}[j(\bar{1})j(\bar{2})]T_{\tau }[j({\bar{1}}^{\prime })j({\bar{2}}^{\prime })]⟩.$$

$${\mathcal{D}}_{1,1^{\prime }}=ϵ\delta ({\mathbf{r}}_1-\mathbf{R})\delta ({\mathbf{r}}_1^{\prime }-\mathbf{R}){\int }_0^{\infty }d\omega e^{-i\omega (t_1-t_1^{\prime })}$$

$$E_{ind}(1)=\int d\bar{1}g(1,\bar{1})j(\bar{1}).$$

$$⟨T_{\bar{\tau }}[\delta j(\bar{1})\delta j(\bar{2})]⟩⟨j({\bar{1}}^{\prime })⟩⟨j({\bar{2}}^{\prime })⟩,⟨\delta j(\bar{1})\delta j({\bar{2}}^{\prime })⟩⟨j({\bar{1}}^{\prime })⟩⟨j(\bar{2})⟩,$$

Full text

Pump-probe experiments are the prime way to study condensed matter out of its equilibrium state on timescales of femto- and picoseconds. In optical pump-probe experiments, ultrashort pulses are used in pairs. The pump triggers the dynamical response and the probe is used to detect changes in the optical properties of the sample. By and large, the experiments performed to date measure the intensity variation of probe pulses after their interaction with the sample for each pump-probe delay, following the protocol of averaging over many stroboscopically repeated experiments Giannetti2016 . Little attention has been given to the accurate measurement of the fluctuation of the intensity of reflected (or transmitted) probe pulses.

In the basic implementation of a “statistical” pump-probe set-up, the intensity of every single probe pulse is separately acquired with low-electronic-noise detectors, for every pump-probe delay. This allows to measure both the average of the intensity, which gives the usual pump-probe signal (e.g. the relative variation of the reflectivity $\frac{\Delta R}{R}(t_{p})$), and its statistical distribution to all orders. Recent technological advances enabled experiments delivering such full statistical information on intensity fluctuations Esposito2015 ; Tesi ; Misochko2001 ; Riek2017 , thereby providing an experimental benchmark to address the following interesting, yet largely unexplored, question: what is the spectroscopic information carried by the fluctuations of the intensity of ultrashort light pulses reflected (or transmitted) by complex materials out of equilibrium?

It is known that the photon counting noise in the radiation emitted by nanoscopic emitters contains valuable information on their steady state transport properties Beenakker2001 ; Lebedev2010 , but for pump-probe measurements in bulk samples it is still unclear how to connect the fluctuations of a measurement which is spatially averaged over the optical beam section to microscopic properties of the material. Proof of principle experiments in this direction have shown that, for transparent materials, i.e. materials where a large optical gap freezes the electronic degrees of freedom and an effective photon-phonon Raman coupling is the leading interaction, the statistical fluctuations of the probe intensity can be used to map the temporal evolution of vibrational observables Benatti2017 ; KelvinThesis ; Esposito2015 ; Tesi . Nevertheless, in the more general setting of pump-probe experiments on complex absorbing materials, where the photo-excitation can trigger a dynamical response in the electronic system, a formalism capable of linking the intensity fluctuations of the probe pulses to the microscopic properties of the material is still unavailable. In this paper we provide such formalism and show how the knowledge of the fluctuations of the optical properties of a system can give access to the spectrum of current fluctuations, which is not accessible from the average intensity in a pump-probe experiment. The general theoretical result will be illustrated with numerical simulations and experimental data for photo-induced coherent lattice vibrations in a charge-density wave system.

Intensity fluctuations — A pump-probe experiment measures the number $m$ of photo-counts in a detector, during a detection time-window which is very large compared to the probe-pulse duration (in general, larger by at least three orders of magnitude). In the statistical pump-probe experiment one obtains in addition the variance $\Delta m^{2}=\sigma[I]+\langle m\rangle$, which can be written as the sum of a shot-noise contribution $\langle m\rangle$ proportional to the intensity itself, and a sample-dependent contribution $\sigma[I]$. Because the induced fields at the detector are linearly related to their sources, i.e., the current density $j$ in the sample, the probability distribution $\mathcal{P}(m)$ can be written in terms of the time-ordered current correlation functions Fleischhauer1998 . In the following we show that the nonequilibrium current fluctuations generated by the pump give rise to a contribution $\sigma_{\text{bulk}}$ to $\sigma[I]$ which can be identified by the statistical pump-probe experiment on a bulk sample, because it is of leading order in the sample volume. This contribution can be understood as the interference term $\langle|\int dtE_{+}(t)\varepsilon(t)|^{2}\rangle_{\text{noise}}$ of the (conventionally averaged) reflected probe field $E(t)$ at the detector ($E_{+}$ is the positive frequency component $\sim e^{-i\omega t}$), a field $\varepsilon(t)$ which is generated at the detector by a classical noise current $\eta$ with a Gaussian distribution and variance given by the current fluctuations in the sample, $\langle\eta(1)\eta(1^{\prime})\rangle_{\text{noise}}=\langle\delta j(1)\delta j(1^{\prime})\rangle$. (Here and in the following, $1\equiv(\bm{r}_{1},t_{1})$ denote space-time points). Since Maxwell’s equations relate the current density $j(1)$ at point $1$ to the induced fields $E_{\text{ind}}(2)=\int d1g(2,1)j(1)$ at a point $2$ via a linear kernel $g(2,1)$, the above statement is equivalent to

[TABLE]

where $\langle j(1)...\rangle_{\text{ret}}\equiv\int d\bar{1}g(1,\bar{1})\langle j(\bar{1})...\rangle$ is the expectation value of the currents propagated to space-time point $1$, $\bm{R}$ is the location of the detector, and $\epsilon$ is proportional to the detector volume and efficiency. To lowest order in the probe amplitude, the correlation function in Eq. (1) can be evaluated without the effect of the probe. Hence the probe pulse can be chosen to project out the transient current noise spectrum of the nonequilibrium state generated by the pump at suitable times $t,t^{\prime}$ (or as a function of time and frequency).

To prove Eq. (1) we start from Ref. Fleischhauer1998 , which gives the factorial moments $I_{1}=\langle m\rangle$, $I_{2}=\langle m(m-1)\rangle$ of the photo-count distribution $\mathcal{P}$ in terms of the sources,

[TABLE]

Here $T_{\tau}$ ($T_{\bar{\tau}}$) is the (anti) time-ordering operator acting on the sources, and $\mathcal{D}_{1,1^{\prime}}$ is the detector response function. For a small and wide-band detector at site $\bm{R}$, which absorbs at all positive frequencies, the latter is simply $\mathcal{D}_{1,1^{\prime}}=\epsilon\delta(\bm{r}_{1}-\bm{R})\delta(\bm{r}_{1}^{\prime}-\bm{R})\int_{0}^{\infty}\!d\omega\,e^{-i\omega(t_{1}-t_{1}^{\prime})}$. Differently from the standard quantum theory of photo-detection Mandel1958 ; Kelley1964 , Eq. (2) is not obtained within the rotating wave approximation, and is thus applicable down to ultra-short (possibly single-cycle) pulses Fleischhauer1998 . The variance of the photo-count is $\Delta m^{2}=\langle m^{2}\rangle-\langle m\rangle^{2}=I_{2}-I_{1}^{2}+\langle m\rangle$, so that $\sigma[I]=I_{2}-I_{1}^{2}$. To evaluate $\sigma[I]$, we shift the operators $j(t)=\langle j(t)\rangle+\delta j(t)$, and expand $\sigma[I]$ in the fluctuations. In a bulk sample, terms like $\int d^{3}{\bm{r}_{1}}d^{3}{\bm{r}_{2}}\cdots\langle\delta j(\bm{r}_{1})\delta j(\bm{r}_{2})\cdots\rangle$ are proportional to the volume $V$, since they are obtained as the derivative of the extensive free energy with respect to an external vector potential. Hence, the dominant contribution ($\sim V^{3}$) to $\sigma[I]$ is given by the second order terms like $\langle j(1)\rangle\langle\delta j(2)\delta j(3)\rangle\langle j(4)\rangle$ (first and zeroth order terms vanish by construction). Since $\langle j(\bm{R},t)\rangle_{\text{ret}}\equiv E(t)$ is just the classical reflected field at the detector, and the contraction with the detector function projects out it’s positive and negative frequency parts, $\int d1^{\prime}\mathcal{D}_{1,1^{\prime}}\langle j(1^{\prime})\rangle_{\text{ret}}\equiv E_{+}(t_{1})$ and $\int d1^{\prime}\mathcal{D}_{1^{\prime},1}\langle j(1^{\prime})\rangle_{\text{ret}}\equiv E_{-}(t_{1})$, one arrives at Eq. (1) footnote1 (see appendix).

While Eq. (1) holds for an arbitrary geometry, to obtain specific results we adopt the standard setup and neglect propagation effects in the sample, so that $E_{\text{ind}}(t,\bm{R})=-\dot{A}(t)$ with $A(t)\propto\int_{\text{sample}}d^{3}\bm{r}\,j(\bm{r},t_{\text{ret}})$, where $t_{\text{ret}}=t-\Delta t$ is shifted by the propagation time $\Delta t$ from the sample to the detector. This gives

[TABLE]

where $J$ is the current integrated over the sample volume.

When the system is in a thermal equilibrium state, the current fluctuations $C(t,t^{\prime})=\langle J(t^{\prime})J(t)\rangle$ are related to the response function $\chi^{R}(t,t^{\prime})=i\theta(t-t^{\prime})\langle[J(t),J(t^{\prime})]\rangle$ of time-resolved optical spectroscopy Eckstein2008 by a fluctuation-dissipation theorem (FDT) Kubo1957 ,

[TABLE]

where $f(t_{p},\omega)=\frac{1}{2\pi}\int ds\,\,e^{i\omega s}f(t_{p}+s/2,t_{p}-s/2)$ is the Fourier transform with respect to difference time, and $b(\omega,\beta)=(e^{\beta\omega}-1)^{-1}$ is the Bose function with inverse temperature $\beta$. In an equilibrium state, we take $\dot{E}_{\pm}(t)\sim e^{\mp i\omega t}$, so that $\sigma_{\text{bulk}}\propto C(\omega)$, which is therefore negligible at room temperature for probe photon energies in the visible and infrared spectral range. Only the shot noise is then relevant, as it is usually assumed in quantum optics for the reflection of light by a mirror. Note that in our language, contributions to $\sigma[I]$ on top of $\sigma_{\text{bulk}}$ which are subleading in the volume $(\sim V^{2})$, can be nonzero in equilibrium, scale like the intensity $E^{2}$ and add to the shot noise of the light coming from an “imperfect reflector”. When the sample is out of equilibrium, instead, $\langle J(t)J(t^{\prime})\rangle$ can significantly contribute to the total fluctuation of the intensity, and its measurement can therefore give information on the nonequilibrium state.

To illustrate such approach both from the theoretical and experimental point of view, we concentrate on a Peierls charge-density wave (CDW) system. A prototypical realization is bismuth Peierls1991 , where an impulsive excitation generates coherent vibrations of the $A_{1g}$ phonon, which modulates the out-of-equilibrium reflectivity Misochko2001 . A minimal theoretical model for the Peierls system is the Holstein model, which describes free electrons coupled to a dispersion-less vibrational mode (Einstein phonon). At half filling, the Hamiltonian is

[TABLE]

where $t_{0}$ is the nearest neighbor hopping, $c^{\dagger}_{j\sigma}$ creates an electron at lattice site $j$ (spin $\sigma$), $b_{j}$ and $b^{\dagger}_{j}$ are the bosonic operators of a phonon mode at frequency $\omega_{0}$, and $g$ is the coupling constant between the electronic density $n_{j}=n_{j\uparrow}+n_{j\downarrow}$ and the phonon coordinate $X_{j}=(b^{\dagger}_{j}+b_{j})$. Below the critical temperature, the system is in an insulating and symmetry broken phase, with a staggered charge disproportionation on neighboring sites and a gap at the Fermi level. We have numerically studied the out-of-equilibrium dynamics of the Holstein model using the non-equilibrium dynamical mean-field theory NonEqDMFT within the self-consistent Migdal approximation Murakami2015 , on a bipartite lattice (see appendix). We choose $g=0.34$, the inverse temperature $\beta=25$, and $\omega_{0}=0.2$. The free bandwidth $W=4$ sets the energy and time-scale ($\hbar=1$). The system is brought out of its equilibrium state by a short modulation of the hopping $t_{0}\to t_{0}(1+\delta h\,e^{-(t-0.16)^{2}/0.32})$, see dotted line in Fig 1a, which impulsively creates a conduction band electron population, analogous to a pump pulse impinging on a sample. (The precise excitation mechanism is not important, since we focus on the subsequent relaxation dynamics, which depends mainly on the excitation density controlled by $\delta h$.)

Before discussing the nonequilibrium current fluctuations, we verify that this model captures the main features which are known from experiments on CDW systems Misochko2001 , i.e., a coherent oscillation which is increasingly damped and softened for larger excitation density. A good observable to study this dynamics is the displacement $\langle X\rangle$ of the atoms from the translationally invariant positions of the high-temperature phase ($\langle X\rangle$ takes opposite values on the two sub-lattices of the bipartite lattice). As shown in Fig. 1a, for $\delta h$ below a critical value $\delta h_{c}\simeq 1.7$, $\langle X(t)\rangle$ coherently oscillates around a new non-zero position. ($\delta h_{c}\simeq 1.7$ corresponds to an pump-induced conduction band population $n_{ex}=0.05$.) As the excitation density increases, the oscillations are damped more quickly, and the frequency $\omega^{*}$ becomes lower, i.e. the vibrational mode is softened (see Fig. 1b). For $\delta h>\delta h_{c}$, the melting of the CDW phase is induced, with oscillations around the undistorted atomic position $\langle X\rangle=0$ (Fig. 1a). While the melting of the CDW phase has been observed in experiment Teitelbaum2016 , we will report measurements below the threshold and therefore focus on the regime $\delta h<\delta h_{c}$ for the following theoretical analysis.

Besides the distortion, we calculate the optical susceptibility $\chi^{R}(t_{p},\omega)$ and the fluctuations $C(t_{p},\omega)$, as defined above Eq. (4). At equilibrium, the optical susceptibility $\chi^{R}(\omega)$ of the CDW phase displays a gapped optical band (Fig. 2a). After the excitation, the out-of-equilibrium $\text{Im}\chi(t_{p},\omega)$ oscillates between a gapped spectrum, which is partially filled by intraband transitions due to the excited conduction band electron population (Fig. 2a, $t_{p}=76$), and an almost gapless spectrum when the atoms are closest to their translationally invariant positions $\langle X\rangle=0$ (Fig. 2a, $t_{p}=52$). The spectrum $C(t_{p},\omega)$ shows that current fluctuations exist only at low frequencies (black line) in equilibrium, but extend to high frequencies when a non-thermal population is created.

A common assumption is that the electrons in a solid, after an intense photo-excitation, can be considered in an effectively thermal state. In systems with strongly coupled electrons and phonons, this assumption need not hold, in particular when phonons itself are far from equilibrium. Theoretically, a check of the FDT (4) provides a natural way to access the thermal nature of a nonequilibrium state Kogoj2016 ; NonEqDMFT . The dashed lines in Fig. 2 show the right-hand side of Eq. (4) for two representative pump-probe delays ($t_{p}=52,76$), with a best-fit inverse temperature $\beta_{f}=7.1$ and $6.3$, respectively. The most pronounced violation of the fluctuation dissipation relation is visible after one full oscillation (yellow curve), i.e. when $\langle X(t)\rangle$ is maximum (see also the blue line in the inset of Fig. 2b, which shows a fit error between $C(t_{p},\omega)$ and the right-hand side of Eq. (4)). One can understand this behavior as follows: Initially, a non-thermal population is created at the top of the lower band and at the bottom of the upper band of the system. As the atoms move closest to $\langle X\rangle=0$ and the gap is the smallest, electrons thermalize and the occupation of the bands closely resembles the Fermi-Dirac distribution, with a high electronic temperature. But when the gap revives after one period, a nonthermal distribution far from the Fermi energy is restored, because reopening of the gap itself prevents electronic thermalization along with the evolution of $\langle X\rangle$.

The measurement of the fluctuations of the intensity of the probe pulses now allows to obtain the fluctuations of the current in the sample experimentally. In the present case, the nonthermal nature of the state after one oscillation period can be detected by probing the fluctuations at a specific high frequency. We evaluate Eq. (3) assuming a generic pulse $\dot{E}(t)\propto\cos(\omega_{p}(t-t_{p}))e^{-(t-t_{p})^{2}/\delta t^{2}}$ with frequency $\omega_{p}=1.3$ larger than the gap, probe time $t_{p}$, and probe pulse duration $\delta t=2.4$ footnote02 The result, shown by the black curve in the inset of Fig. 2b, closely resemble the error in the fit of the FDT. These numerical predictions can be compared to experimental results obtained in experiments on Bi single crystals with the statistical pump-probe set-up Esposito2015 . The system is excited and probed at a photon energy of 1.5 eV, with a probe duration of 80 fs. Figures 3a and b show the measured variance $\Delta m$ and $\operatorname{\frac{\Delta R}{R}}$ for two different fluences. The pump generates a coherent amplitude mode oscillation, which is seen in $\operatorname{\frac{\Delta R}{R}}\propto\langle X\rangle$. Also $\Delta m$ mainly follows this oscillation. This is expected because in addition to the contribution $\sigma_{\text{bulk}}$ from the nonequilibrium current fluctuations, $\Delta m$ contains shot noise proportional to the intensity itself. Moreover, there is a slight variation with $\operatorname{\frac{\Delta R}{R}}$, because the measurement of $\Delta m$ is calibrated to the equilibrium reflectivity of the sample. These contributions to $\Delta m$ are quantitatively unknown, but they are proportional to $\operatorname{\frac{\Delta R}{R}}$, and one can clearly see that a rescaled $\operatorname{\frac{\Delta R}{R}}$ cannot describe the time-dependent fluctuations $\Delta m$. In particular, we find an enhancement of the fluctuations at early times, and, at higher fluence, after one oscillation period, i.e., around the time when the atoms revive to their most distorted positions. This is well in agreement with the theoretical prediction for the fluctuation contribution $\sigma_{\text{bulk}}$ which comes from the current noise. The experimental data therefore show the that electrons have not thermalized at each time during the early CDW oscillations of the system.

In conclusion, in this letter we have derived an expression for the fluctuations of the intensity of light reflected by a sample in an ultra-fast pump-probe experiment. We have shown that fluctuations of the reflectivity in addition to shot noise contributions can measure the nonequilibrium current noise spectrum as a function of time and frequency on the femtosecond timescale. In a most straightforward way, current fluctuations can reveal the nonequilibrium state of the sample, such as the effective temperature of the relevant degrees of freedom, or whether the system can be considered in a quasi-thermal state at all. We illustrated our general results for the case of coherent vibrations in solids, both via numerical calculations and statistical pump-probe experiments on bismuth single crystals. We find current fluctuations at high energies, which reveals that nonthermalized electronic distributions emerge as the CDW gap oscillates. An intriguing future direction would be to also characterize the quantum state of the low energy mode. While the theoretical simulations predicts that is does not only coherently oscillate but also becomes squeezed, this is only indirectly reflected in the probed response at high energies. In general, the results of this manuscript show that measuring the fluctuations in the reflectivity opens up an independent way to characterize the nonequilibrium dynamics in bulk solids, with intriguing future applications in the characterization of transient superconducting or nonthermal symmetry broken states.

Appendix A Nonequilibrium DMFT for the Holstein model

In this section we present the nonequilibrium DMFT setup for the Holstein model. Apart from the symmetry breaking, the formalism is analogous to what has been explained in Ref. Murakami2015 . We therefore only state the equations, and do not provide a detailed derivation.

In DMFT, the Holstein model [see Eq. (5) of the main text] is mapped to a set of Anderson-Holstein impurity models (one for each inequivalent lattice site), with action

[TABLE]

on the Keldysh time-contour $C$. (For an introduction to nonequilibrium DMFT and to the Keldysh formalism, see Ref. NonEqDMFT ). In Eq. (6), the first term is the local part of the lattice Hamiltonian, which involves the coupling of the electrons at the impurity site to the coordinate $X=(b^{\dagger}+b)/\sqrt{2}$ of the local oscillator, and $\Delta(t_{1},t_{2})$ is the hybridization function, which is determined self-consistently below.

The impurity model is solved using the self-consistent Migdal approximation Murakami2015 , where also the vibrational mode evolves as a consequence of the interaction with the electrons. In the symmetry broken phase, the coordinate $X$ acquires a nonzero expectation value. The expectation value is determined by the exact equation of motion $\frac{d^{2}}{dt^{2}}\langle X(t)\rangle=-\omega_{0}^{2}\langle X(t)\rangle+F(t)$, with the time-dependent force $F(t)$

[TABLE]

In turn, there is a time-local (Hartree) contribution to the electronic self-energy, i.e., a self-consistent on-site potential,

[TABLE]

Furthermore, we include the leading order self-consistent diagrammatic corrections in the expansion in terms of the fluctuations $\tilde{X}=X-\langle X(t)\rangle$. The second-order electronic self-energy is

[TABLE]

where

[TABLE]

(We consider the spin-symmetric phase and omit spin indices $\Sigma_{\sigma}$ and $G_{\sigma}$.) With this, the Dyson equation for the electronic Green’s function reads

[TABLE]

To include the back-action of the electrons on the phonons on the same diagrammatic level, we include the phonon self-energy (polarization operator)

[TABLE]

and solve the phonon Dyson equation in the form

[TABLE]

Here $D_{0}(t,t^{\prime})$ is the non-interacting phonon propagator,

[TABLE]

where $b_{\beta}=1/(e^{\beta\omega_{0}}-1)$ is the Bose function.

In the present case of a two-sublattice symmetry broken phase, we have two inequivalent impurity models (6), which represent sites on the $a$ and $b$ sublattice, i.e., all quantities, $G$, $\Delta$, $h_{loc}$, $\langle X(t)\rangle,$ $\Sigma$, $P$, will additionally depend on the sublattice $a,b$. For the particle-hole symmetric case, we have $\langle X(t)\rangle_{a}=-\langle X(t)\rangle_{b}$. We use a bipartitle lattice with a semielliptic density of states, in which the DMFT self-consistency is given by NonEqDMFT

[TABLE]

where $v(t)$ is the time-dependent profile of the hopping amplitude. This closes DMFT equations.

Optical conductivity: The main quantity of interest in this work is the current correlation function $C(t,t^{\prime})=\langle j(t)j(t^{\prime})\rangle$, and the optical susceptibility $\delta\langle j(t)\rangle/A(t^{\prime})$, which characterizes the long-wavelength current ($q\rightarrow 0$) in response to an applied time-dependent vector potential $A(t)$. Since the current operator is given by $j=-\delta H/\delta A$, the latter response function is given by $\chi^{R}(t,t^{\prime})=i\theta(t,t^{\prime})\langle[j(t),j(t^{\prime})]\rangle$. (For simplicity of notation, we are omitting cartesian components $x,y,z$.)

The current-current correlation function is obtained from the lattice Green’s function by direct generalization of the expressions presented in Ref. Eckstein2008 . Both response and correlation function are obtained from the contour-ordered current-current correlation function,

[TABLE]

which is defined as the response of the current to an arbitrary variation of the vector potential along the Keldysh contour,

[TABLE]

omitting a diamagnetic contribution which is time-local and thus irrelevant for the discussion of the dynamic properties discussed in this paper. Current fluctuations are given by the greater and lesser component, $\chi^{>}(t,t^{\prime})\equiv\chi(t_{-},t^{\prime}_{+})=i\langle j(t)j(t^{\prime})\rangle$ and $\chi^{<}(t,t^{\prime})\equiv\chi(t_{+},t^{\prime}_{-})=i\langle j(t^{\prime})j(t)\rangle$ ($t_{\pm}$ is on the upper/lower branch of the Keldysh contour), and $\chi^{R}(t,t^{\prime})=\theta(t,t^{\prime})(\chi^{>}(t,t^{\prime})-\chi^{<}(t,t^{\prime})$.

In the symmetry broken phase, the lattice has a unit cell with two sites $a,b$, and a reduced Brilluoin zone (RBZ). We introduce the spinor,

[TABLE]

and the momentum-dependent Green’s function $G_{k}$ then becomes a $2\times 2$ matrix, $\hat{G}_{k}(t,t^{\prime})=-i\langle T_{C}\hat{\psi}_{k}(t)\hat{\psi}_{k}^{\dagger}(t^{\prime})\rangle$. The hopping term takes the form $H_{hop}=\sum_{k}\hat{\psi}_{k}^{\dagger}\epsilon_{k-A}\hat{\sigma}_{1}\hat{\psi}_{k}$, with the Pauli matrix $\hat{\sigma}_{1}$; $\sum_{k}$ is a sum over the reduced Brillouin zone (RBZ). The vector potential is added by the Peierls substitution $\epsilon_{k}\to\epsilon_{k-A}$, so that the current $j=-\delta H/\delta A$ is given by

[TABLE]

where $v_{k}=\partial_{k}\epsilon_{k}$ is the band velocity, and the factor $2$ is for spin. Following Ref. Eckstein2008 , we take the variation $\delta A(t^{\prime})$, using that vertex corrections to the current correlation function vanish in DMFT. This gives the susceptibility (evaluated at $A=0$),

[TABLE]

which is the usual bubble diagram of the Green’s functions.

In DMFT (and when we consider only a modulation of the hopping amplitude, as in the manuscript), $\hat{G}_{k}$ depends on $k$ only via the dispersion $\epsilon_{k}$, i.e., $\hat{G}_{k}(t,t^{\prime})\equiv\hat{G}_{\epsilon_{k}}(t,t^{\prime})$. In the particle-hole symmetric case on a bipartite lattice, the RBZ corresponds to positive values of $\epsilon_{k}$. The momentum sum can then be represented by integrals

[TABLE]

where $\rho$ and $D$ depend on the lattice. We assume a semi-elliptic density of states $\rho(\epsilon)=\sqrt{4-\epsilon^{2}}$, and the corresponding form for $D(\epsilon)$ as defined in Ref. Eckstein2008 . The lattice Green’s function is evaluated on a grid of momentum points, solving the Dyson equation $\hat{G}_{\epsilon_{k}}(t,t^{\prime})=(i\partial_{t}+\mu-\hat{h}(t)-\hat{\Sigma}(t,t^{\prime}))^{-1}$, where $\hat{h}$ and $\hat{\Sigma}$ in the {a,b} basis are

[TABLE]

with the sublattice-dependent $\Sigma_{a,b}$ and $\langle X\rangle_{a,b}$.

Appendix B Derivation the reflectivity fluctuations

In this section we present explicit steps of the fluctuation expansion leading from the general expression for the moments $I_{n}$ of the photon-count [main text, Eq. (2)] to the variance of the intensity [main text, Eq. (1)]. We start from the expression for $I_{n}$,which was been derived by Fleischhauer Fleischhauer1998 ; footnote ,

[TABLE]

Here

[TABLE]

is the detector response function, and $g$ is the linear kernel which relates the induced field $E_{ind}$ and the current $j$ by a solution of Maxwell equations,

[TABLE]

In Eq. (26) and (27), we insert the expansion $j(1)=\langle j(1)\rangle+\delta j(1)$. Terms which are first order in $\delta j$ vanish by construction, because $\langle\delta j\rangle=0$. Zeroth order terms are identical in $I_{2}$ and $I_{1}^{2}$, and thus vanish in the variance $I_{2}-I_{1}^{2}$. Third and fourth order terms, such as Eq. (27) where the correlation function in the integrand is replaced by $\langle T_{\bar{\tau}}[\delta j(\bar{1})\delta j(\bar{2})]\delta j(\bar{1}^{\prime})\rangle\langle j(\bar{2}^{\prime})\rangle$, are not considered here as explained in the main text, because they scale differently with the sample volume. To second order in $\delta j$, Eq. (27) has six terms, where the current correlation function in the integral takes one of the following combinations

[TABLE]

Here the time-ordering operator can be dropped whenever it acts on c-numbers $\langle j\rangle$, such as for the second term, $\langle T_{\bar{\tau}}[\delta j(\bar{1})\langle j(\bar{2})\rangle]T_{\tau}[\delta j(\bar{2}^{\prime})\langle j(\bar{1}^{\prime})\rangle]\rangle=\langle\delta j(\bar{1}^{\prime})\delta j(\bar{2})\rangle\langle j(\bar{1})\rangle\langle j(\bar{2}^{\prime})\rangle$. Of the six terms in Eq. (30), the last two are cancelled by corresponding terms in the expansion of $I_{1}^{2}$. For the remaining four, one can evaluate integrals in (27) which correspond to a contraction of the current expectation values $\langle j\rangle$ with $\mathcal{D}$,

[TABLE]

using Eqs. (28) and (29) in the first step. Similarly,

[TABLE]

Inserting Eqs. (30), (31), and (32) into Eqs. (27) and (26) we get

[TABLE]

[TABLE]

In the main text the first term is not discussed because it would vanish by averaging over a carrier envelope phase $\varphi$ ($E_{\pm}\sim e^{\mp i\varphi}$). The remaining term is Eq. (1) of the main text.