Direct measurement of key exciton properties: energy, dynamics and spatial distribution of the wave function

Excitons, Coulomb-bound electron-hole pairs, are the fundamental excitations governing the optoelectronic properties of semiconductors. While optical signatures of excitons have been studied extensively, experimental access to the excitonic wave function itself has been elusive. Using multidimensional photoemission spectroscopy, we present a momentum-, energy- and time-resolved perspective on excitons in the layered semiconductor WSe$_2$. By tuning the excitation wavelength, we determine the energy-momentum signature of bright exciton formation and its difference from conventional single-particle excited states. The multidimensional data allows to retrieve fundamental exciton properties like the binding energy and the exciton-lattice coupling and to reconstruct the real-space excitonic distribution function via Fourier transform. All quantities are in excellent agreement with microscopic calculations. Our approach provides a full characterization of the exciton properties and is applicable to bright and dark excitons in semiconducting materials, heterostructures and devices.


Abstract
Excitons, Coulomb-bound electron-hole pairs, are the fundamental excitations governing the optoelectronic properties of semiconductors. While optical signatures of excitons have been studied extensively, experimental access to the excitonic wave function itself has been elusive.
Using multidimensional photoemission spectroscopy, we present a momentum-, energy-and time-resolved perspective on excitons in the layered semiconductor WSe 2 . By tuning the excitation wavelength, we determine the energy-momentum signature of bright exciton formation and its difference from conventional single-particle excited states. The multidimensional data allows to retrieve fundamental exciton properties like the binding energy and the exciton-lattice coupling and to reconstruct the real-space excitonic distribution function via Fourier transform.
All quantities are in excellent agreement with microscopic calculations. Our approach provides a full characterization of the exciton properties and is applicable to bright and dark excitons in semiconducting materials, heterostructures and devices.
Excitons, bound electron-hole quasi-particles carrying energy and momentum but no net charge, are fundamental excitations of semiconductors and insulators arising from light-matter interaction. 1 An initial excitonic polarization induced by a light field (often referred to as coherent excitons and e.g. detected by optical absorption spectroscopy) rapidly loses coherence with the driving field and dephases into a population of incoherent excitonic states. 2,3 The generated excitons propagate in solid-state materials through diffusion 4,5 and eventually release their energy e.g. in the form of luminescence (photon), lattice excitation (phonon) or dissociation into single charged quasi-particles. [6][7][8][9] Understanding exciton physics is of capital importance for advanced photonic and optoelectronic applications including photovoltaics. Layered transition metal dichalcogenide (TMDC) semiconductors exhibit rich exciton physics even at room temperature due to strong Coulomb interaction. 10 Excitons in TMDCs feature large oscilla-tor strength 11 and their inter-and intra-band dynamics have been extensively investigated. [12][13][14] Moreover, strong spin-orbit coupling and broken inversion symmetry in each crystalline trilayer lead to a locking between spin, valley and layer degrees of freedom, which started a surge of valley physics studies. [15][16][17] A large portion of the research on excitonic phenomena in TMDCs adopts optical spectroscopic techniques, 10,12,[14][15][16][18][19][20][21] which only access bright excitonic transitions with near-zero momentum transfer. While techniques such as time-resolved THz spectroscopy also allow probing optically dark excitons via internal quantum transitions, 13 finite-momentum excitons which lie outside the radiative light cone remain inaccessible to such methods. This limitation is overcome by time-and angle-resolved photoemission spectroscopy (trARPES), a spectroscopic tool accessing excited states, including excitons, in energy-momentum space and on ultrafast timescales. 17,[22][23][24] Here, we reveal the characteristics of the excitonic wave function in the photoemission signal of the prototypical layered TMDC semiconductor 2H-WSe 2 and establish that all fundamental exciton properties are encoded in the trARPES signal's energy, time, and momentum dimensions: the exciton binding energy, its self-energy as a measure of the excitonlattice coupling, as well as the real-space distribution of the excitonic wave function. given energy-momentum-plane with high counting statistics using the HA, and alternatively resolve both in-plane momentum directions yielding a 4D photoemission signal I(E kin , k x , k y , t) of the entire valence band with the MM. 22,25 Fig.1 (b-d) and (e-g) show snapshots of the 3D and 4D data with 1.55 eV excitation, respectively, at three selected time delays: i) prior to optical excitation, showing the ground-state band structure of WSe 2 from the Brillouin zone (BZ) cen- With the HA, the equilibrium band structure in a finite momentum window is found at negative pump-probe delay, t = −100 fs, showing the spin-orbit split valence bands near the K-points. Upon 1.55 eV excitation, the excited state dynamics, i.e., the intervalley scattering from K to Σ valley, are representatively shown in (c) t = 0 fs and (d) t = 100 fs. (e) Four-dimensional (4D) band structure mapping, I(E, k x , k y , t), with the MM showing the band dispersion within the whole Brillouin zone from its center Γ to the K valleys at its corners. The same evolution of the excited state is shown for (f) t = 0 fs and (g) t = 100 fs, respectively. All the excited states signal are scaled for clarity.
ter Γ (only shown in the MM data) to the BZ boundary K points (b,e); ii) upon optical excitation resonant with the A exciton absorption (the first excitonic state), featuring excited-state signal at the K and Σ valleys (c,f); and iii) at t = 100 fs after optical excitation, with excited-state signal mostly at the Σ valleys (d,g). In the following, we identify the excitonic features in the excitedstate photoemission signal and quantify the exciton properties retrieved from the energy, time and momentum dimensions of the 4D trARPES signal.
Photoemission signature of excitons. During the photoemission of an electron bound in an exciton, the electron-hole interaction diminishes, i.e. the exciton breaks up, as a single-particle photoelectron is detected while a single-particle hole is left behind in the material. 23 To identify the signature of the excitonic electron-hole interaction in photoemission spectroscopy, we compare the signal of excitons with that of single-particle excited states. For generating excitons, we excite with 1.55 eV photons (1/e 2 bandwidth = 43 meV), in resonance with the low-energy side of the A-excitonic absorption of bulk WSe 2 . 12 Fig.2(a) shows the excited-state signal integrated in the first 25 fs after pump-probe overlap. The data reveals a vertical transition at the K point through an excited-state signal localized in energy and momentum. In contrast, the above-band-gap excitation with 3.1 eV photon energy generates a population higher in the conduction band, which rapidly redistributes to all bands and valleys of the lower conduction band (the equivalent scenario applies to the holes in the valence band). Fig.2(b) shows this excitedstate signal in the K valley 100 fs after excitation, where carriers have redistributed in energy and momentum. This signal resembles the dispersion of a single-particle band with an effective mass of m * = 0.55 m e , in good agreement with electronic band structure calculations. 26 The energetic positions of the excitonic and single-particle states at the K point are determined as the center of mass of the energy distribution curves (EDCs), see Fig.2(b). The excited-state signal upon resonant excitation of the A exciton is centered 100 ± 3 meV below the center of the single-particle band. Such a signal below the single-particle band has been predicted as photoemission signature of excitons and the energy difference can be associated with the exciton binding energy E b . 2,3,9,27 A calculation of the A-exciton binding energy in bilayer WSe 2 based on the screened Keldysh-like potential (see SI for details) yields E b = 91.3 meV, in very good agreement with the experimental value. It is important to note that we retrieve the exciton binding energy directly from measuring the absolute energies of many-body and single-particle states with a single photoemission experiment, in contrast to combining different experimental methods 18 or by comparing photoemission signals with electronic structure calculations. 24 The observation that the excitonic binding energy is measurable as energy loss of the photoelectron confirms that the hole final states indeed are identical to single-particle holes, which further implies that the localized hole of the exciton transforms to Bloch-like single-particle states during the photoemission process.
To set the stage for discussing the exciton dynamics, we emphasize a signal appearing as replicas of the upper (VB1) and lower (VB2) spin-orbit split valence bands in Fig.2(a) shifted by the photon energyhω pump . This signal only appears during temporal pump-probe overlap and we attribute it to a photon-dressed electronic state due to coherent coupling to the optical driving field. Since the employed s-polarized pump light (polarization parallel to the sample surface) suppresses laser-assisted photoemission, the experimental configuration selectively probes the coherent excitonic polarization induced by the pump field. 2,3,28 Formation and decay dynamics of bright excitons. The bright A excitons at the K point are not the lowest-energy excitons in WSe 2 but can relax their energy further by scattering in momentum space. We extract the quasiparticle dynamics within three regions of interest (ROIs) from the trARPES data in Fig.2(a), representing the coherent excitonic replica of VB1, the excitonic state at the K valley, and the Σ valley population. The respective time traces in Fig.3(a) reflect three types of quasi-particle dynamics: the dephasing of the coherent excitonic polarization (black), the buildup and relaxation of a bright exciton population at K (blue) and the carrier accumulation of dark states at Σ (yellow).
The observed carrier dynamics imply the following microscopic processes as sketched in Fig.3(b). First, the interaction of the initial valence band state |i with the near-resonant optical light field creates a coherent excitonic polarization (dashed line), which quickly dephases into an optically bright exciton population |n , offset by the pump detuning ∆ a . The decoherence Figure 2: Signatures of excitons versus quasi-free carriers in trARPES. (a) Upon the arrival of the pump photons athω = 1.55 eV, the excited states at K and Σ are populated. During the pump-probe overlap, sideband replica of the two topmost valence bands, VB1 and VB2, are visible (highlighted by the green dashed lines). Inset: schematic of the bright exciton at the K valley and dark exciton at the Σ valley. (b) Using the above-band-gap pump at 3.1 eV, the parabolic dispersion of the conduction band at the K valley is clearly observed. As shown in the EDCs at K (right panel), the single-particle (3.1 eV pump photon energy; blue) is ∆E ≈ 100 meV higher than the excitonic signal at 1.66 eV (1.55 eV pump photon energy; red). The excited state signals are scaled for clarity. process occurs with the pure dephasing time T * 2 . These bright excitons undergo rapid scattering into the optically dark Σ-point state |s on the timescale T 1 . We model these processes and the photoemission signals from these states into the continuum final states |f and |g using a five-level extension to the optical Bloch equations 29,30 (OBE, see SI). Based on a multivariate least-squares fitting procedure, we can describe the dynamics of coherent and incoherent exciton contributions, obtaining a coherent exciton dephasing time of T * 2 = 17 ± 9 fs and a population lifetime for the bright A-exciton population of T 1 = 18 ± 4 fs. The extracted dephasing time corresponds well to microscopic calculations. 14,31 To evaluate the mechanism governing the bright exciton scattering, we performed ab initio calculations of the single-particle self-energy of WSe 2 . At low excitation densities, the electron self-energy is dominated by electron-phonon interaction which is computed using density functional perturbation theory (DFPT), taking into account the electronic screening of the lattice motion (see SI). The imaginary part of the momentum-resolved self-energy is shown in Fig.3(c) encoded by the color scale. From the calculation, we obtain Im(Σ el−ph ) = 13.1 meV at the conduction band minimum and Im(Σ h−ph ) = 2.6 meV at the valence band maximum of the K valleys. While a rigorous description of exciton-phonon coupling requires treatment on the basis of excitonic eigenstates, in the weak coupling limit, i.e., small self-energy renormalization due to the electron-hole interaction, the exciton-phonon self-energy is dominated by its incoherent contribution. 32 In this case, the exciton-phonon interaction can be approximated as sum of the single-particle-phonon interactions. Our calculated value Im(Σ el−ph )+Im(Σ h−ph ) ≈ 16 meV agrees with the experimental exciton self-energy Im(Σ ex ) = 18 ± 4.8 meV determined according to Im(Σ ex ) =h/(2 T 1 ). This agreement with theory shows that the exciton lifetime provides a quantitative measure of the strength of its interaction with the lattice and supports the assumption of a dominating incoherent self-energy contribution.
Momentum-and real-space distribution of A excitons. Our 4D trARPES data not only provides the energy-momentum dynamics of excitons but also contains direct amplitude information about exciton wave functions. In Fig.4(a), we display the early-time excited-state momentum distribution I(k x , k y , t = 0 fs) of the K valleys, by integrating in energy over the CB. Signals from other valleys are filtered out in order to focus on the A excitons (see SI). The total photoemission intensity is proportional to the squared transition dipole matrix element, which connects the initial state wave function ψ i to the photoemission final state ψ f , via the polarization operator A · p. Here, A is the vector potential of the light field and p is the momentum operator. Within the plane wave approximation (PWA) for the final state, the matrix element takes the form where k is the wave vector of the photoionized electron. According to Eq.1, the matrix element is proportional to the amplitude of the Fourier transform (FT) of the initial state wave function.
Therefore, the momentum distribution of the photoemission signal I(k x , k y ) can be used to retrieve the real-space probability density of the electron contribution to the two-particle excitonic wave function, i.e., the modulus-squared wave function, I(r x , r y ), with a suitable assumption for the missing phase information.
A similar reconstruction of electronic wave functions from ARPES spectra has been previously demonstrated for occupied molecular orbitals in the ground state of crystalline organic films and chemisorbed molecular monolayers. 33,34 Here, we extend this technique into the time domain and apply it to reconstruct the excitonic wave function in WSe 2 . Assuming a constant phase profile across the BZ as a lower-limit wave function extension (see SI), we retrieve the exciton probability density via 2D FT as shown in Fig.4(a-b). The reconstruction exhibits a broad isotropic real-space exciton distribution carrying high-frequency oscillations, corresponding to the hexagonal periodic lattice structure of WSe 2 . To resolve the isotropic exciton wave function envelope more clearly, the 1D real-space carrier distribution without the oscillatory pattern is shown in Fig.4(d), obtained by FT of only one of the six K valleys, yielding a value of r exp WSe 2 = 1.74 ± 0.2 nm for the excitonic Bohr radius.
To verify the method of reconstructing excitonic wave functions, we performed microscopic calculations of trARPES spectra. The momentum-resolved description of the exciton is based on a many-particle treatment of the Coulomb interaction between electron-hole pairs and the exciton-phonon scattering dynamics 3 (see SI). The momentum distributions of the bright Kexcitons calculated within the PWA for the final state is shown in Fig.4(c). We find a very good agreement to the experimental momentum distribution curve (MDC) taken along the dashed line in Fig.4(a)), supporting our assumption that the trARPES spectrum contains the fingerprints of the excitonic wave function and justifying the use of the PWA. Furthermore, the calculated In the SI, we reconstruct the real-space exciton density distribution with non-constant intervalley and intravalley phase profiles, where we find a broadened exciton distribution in the case of an intravalley varying phase. Therefore, we note here that the real-space reconstruction of the exciton density with a constant phase is suitable for topologically trivial solid-state wave functions. However, the winding of the phase in topologically non-trivial materials leads to an additional expansion of the carrier density distribution, requiring explicit momentum-dependent phase information. In general, the phase of the excitonic wave function might additionally be reconstructed through iterative phase retrieval algorithms. 35 We envision that future developments will allow retrieving the phase as well as orbital information of excitonic wave functions by utilizing dichroic observables 36-38 in trARPES.
In this work, we provide a comprehensive experimental characterization of an excitonic state with trARPES. The interactions governing the formation of this prototypical many-body state are observable as energy renormalization in comparison to single-particle states, while its interaction strength with other quasi-particles is reflected in the excited state's lifetime. These quantities are intimately connected to the real and imaginary parts of the many-body state's selfenergy and our approach establishes experimental access to these elusive quantities. Moreover, we retrieve real-space information of the excitons by Fourier transform of its momentum distribution, establishing the measurement of wave function properties of transient many-body states with 4D photoemission spectroscopy. Our approach is applicable to all exciton species occurring in a wide range of inorganic and organic semiconductors, van der Waals heterostructures and devices. Its extension to other many-body quasi-particles in solids appears straightforward.
Data availability We provide the full experimental dataset as well as the details of the data analysis on the data repository Zenodo. Also, we provide the source code of our data analytics on GitHub. H.H. and A.R. calculated the single-particle self-energy. R.P.X. developed the 4D data processing code. All authors contributed to the manuscript. Competing Interests: The authors declare that they have no conflict of interest.

Supplementary Notes
This section includes the experimental methods, procedure of reconstructing the real-space excitonic wavefunction and energy calibration of trARPES spectrum.

Experimental Methods
Time-and angle-resolved photoemission spectroscopy. The whole setup and the computational workflow for data processing have been described in detail elsewhere. 1-3 Our laser system is a home-built optical parametric chirped-pulse amplifier (OPCPA) delivering 15 W (λ pump = 800 nm) at 500 kHz repetition rate. 4 The major part (80%) of the OPCPA output is used to drive high-order harmonic generation (HHG) by tightly focusing the second harmonic of the laser pulses (400 nm) onto a dense Argon gas jet. Out of the generated XUV frequency comb, a single harmonic (7th order, 21.7 eV) is isolated by a combination of a mulilayer mirror and propagation through a 400 nm thick Sn metallic filter. The remaining part of the OPCPA output serves as the optical pump beam, with a transform-limited pulse duration of 35 fs. Another pump wavelength used in the experiment, λ pump = 400 nm, is the frequency-doubled fundamental pump light generated using a barium borate crystal. In the measurement, the pump fluence of 800 nm is F 800 = 1.3 mJ/cm 2 and that of 400 nm is F 400 = 85 µJ/cm 2 . All measurements are performed at room temperature. The optical pump and probe beams are focused at the sample position in the ultra-high vacuum chamber which is equipped with two type of photoelectron analyzers, a conventional hemispherical analyser (HA) from SPECS GmbH and a novel time-of-flight momentum microscope (MM) from Surface Concept and SPECS GmbH. This combination of complimentary electron analyzers enables high quality and efficient data collection within the full Brillouin zone (with the MM) and regions of interests (with the HA), and the exploration of ultrafast dynamics. The MM collects photoelectrons within a wide emission angle using an extraction lens, simultaneously recording the photoemission spectrum at both in-plane momentum directions k x and k y using a delay line detector. In a pump-probe scheme, the MM thus directly provides the 4D-dimensional photoemission intensity I(E kin , k x , k y , t) in an efficient way as shown in the Fig.1. In contrast, the HA yields a single energy-momentum cut in a fixed experimental configuration, effectively allowing for a much higher electron detection rate within a particular momentum range, which allows us to investigate the delicate quasiparticle scattering dynamics with enhanced signal-to-noise ratio.
In our experiment, the energy axis of trARPES spectra are aligned with the ground state of valence band maximum (VBM) at the K valley.
Sample preparation. Bulk WSe 2 is a purchased crystal from HQ Graphene, which is firstly glued on top of a copper sample holder and then cleaved at room temperature and a base pressure of 5x10 −11 mbar. The sample is further handled by a 6-axis manipulator (SPECS GmbH) for trARPES measurements.

Determination of the excitonic distribution function
Matrix element effects in the ARPES spectrum Within the sudden approximation, the intensity of the ARPES spectrum can be written as the product of a transition matrix element M k f,i , the one-electron removal spectral function A(k, ω) and the electronic distribution function f (k, ω): As discussed in the main text, the matrix element is defined as a transition between the initial and final state wave functions|M k f,i | 2 = | ψ f |A · p|ψ i | 2 , mediated by the vector potential of the exciting electromagnetic field A and the momentum operator p. Assuming the final state as a plane wave, this component can be simplified as |M k f,i | 2 ∝ |A · k| 2 | e ik·r |ψ i | 2 , containing two elements: (1) the polarization factor |A · k| between the potential of the optical field A and the photoelectron wave vector k, and (2) the Fourier transformation of the initial state wave function ψ i (k). The first term is largely determined by the experimental geometry, leading to a slight modulation of the spectral weight based on the projection of polarization to the final state wave vector. 5 Therefore, it would unavoidable introduce a intensity modulation in the practical use of hemispherical electron analyser, because of the sample rotations during the momentum scan. In our setup, the momentum microscope eliminates the experiment-induced intensity variation by collecting the photo-ionized electrons of a wide angle range, such that the momentum map can be achieved in a fixed experimental geometry. The matrix element component is naturally encoded with the information of initial state as discussed in the main text.
The photoemission intensity modulation due to the orbital texture has been demonstrated in the recent study of time-reversal dichroism in ARPES. 6 Under the plane wave assumption (PWA) of final states, the real-space carrier reconstruction using its momentum distribution achieves the experimentally long-pursued goal: mapping the fundamental shape and size of electronic wavefunction. Nevertheless, the validity of the Fourier imaging approach, more specifically speaking, the PWA of the final states has been continuously debated, which is considered as applicable under the favorable conditions of simple orbital components and high photon energy. 7 The good agreement of the microscopy calculation of trARPES spectrum using PWA of final states with the experimental results, shown in the main text Fig.4, supports the viability of real-space reconstruction of exciton in our system. The data preparation procedure The momentum maps as shown in Fig.4 (a)  and some of the Σ valleys of the second BZ ( Fig. S.1(a)). By applying a circular momentum mask, we isolate the six K valleys as shown in Fig. S.1(b). Next, the population intensity at each K valley is normalized by the local maximum to remove the geometric polarization factor ( Fig. S.1(c)), which comes from the coupling between light field and photoionized electron determined by the experimental geometry. In this measurement, we used linearly polarized light at 65 • incidence angle with respect to the sample surface. In step (ii), the experimental momentum resolution effect is deconvolved from the measured data to obtain the intrinsic momentum distribution. Our momentum resolution of METIS is, δk = 0.063Å  S.1(g)). Noted, we present the intrinsic momentum distribution of six K valleys with the averaged lineshape after removing the momentum resolution effect. Within the 2D FT process, we transfer the square root of the intensity of the 2D momentum map I(k x , k y ) to the real-space wavefunction ψ(r x , r y ) and then square the wave function to present the probability |ψ(r x , r y )| 2 . The neglected phase information will be discussed below.
In analogy to the electron Bohr radius in hydrogen atom, 8 we estimate the exciton Bohr radius based on the average radius of exciton distribution r 2 = |r · ψ(r)| 2 . The 1D realspace exciton distribution |ψ(r)| 2 is shown in the main context Fig.4(d) and the exciton Bohr radius r WSe 2 is defined as the distance of the peak intensity in r 2 . By calculating the exciton distribution at six K valleys, we obtain the exciton Bohr radius r WSe 2 = 1.74 ± 0.2 nm, which the standard deviation is a measure of variability between valleys.
Momentum dependent phase profile of exciton wave function In the FT procedure described in the main text, we assumed a constant phase profile of the excitonic wave function, which yields a lower limit for the excitonic distribution function. This relation will be rational- While more advanced approaches such as phase-retrieval schemes have been employed to reconstruct the amplitude and phase of ARPES spectra in iterative algorithms, 9 here we discuss the influence of inter-and intravalley-dependent phase variations based on symmetry considerations and simulations. With linearly polarized excitation, bulk WSe 2 preserves both time reversal symmetry and spatial inversion symmetry. The time reversal operatorT introduces a phase reversal between K and K' point, which suggests the Bloch wave function can be written as ψ K (r) = e ik·r φ K (r)e iθ and ψ K (r) = e ik·r φ K (r)e −iθ , respectively.
To fulfill these conditions, we construct a periodical varying phase mask of e iθ and e −iθ  with the momentum-independent phase profile ( Fig. S.2(a)) which was applied for the evaluation in the main text, we find no influence on the width of the exciton distribution ( Fig. S.2(h)).
Note that the high frequency oscillations under the envelope function vary with the intervalley phase reversal, which can be regarded as analogous to the interference pattern changing with the relative phase differences between two coherent emitters (K and K' valleys).

Optical Bloch equation fitting
Optical Bloch equation (OBE) modelling has been employed to describe two-photon photoemission from metallic surface states. 10,11 In general, the use of OBE in solid state systems and in particular in photoemission hinges on the important result that the coherent excitation of a quasi-continuum of states may be treated in the same way as the incoherent limit of a pure two level system, significantly simplifying the full harmonic description of the quantum process. 12 It is therefore possible to apply OBE to photoemission processes within semiconductors when the conditions allow schematizing the system as a set of atomic levels. In particular, it is well suited for pump excitation wavelengths very close to the direct bandgap of the semiconductor. In this conditions, hot carriers have very little excess energy and intra-valley relaxation dynamics can be neglected.
OBE offer a tool to disentangle the formation and decay of a coherent polarized state of the solid from the intrinsic lifetime of electronic states. This allows for a much more precise assessment of intrinsic lifetimes than models based on rate equations that ignore the decoherence processes, and tend to grossly overestimate the real values, especially when the lifetime is close to, or even shorter than, the temporal resolution, i.e. pump-probe pulse cross-correlation. Here we explain how the three level OBE model is extended in order to determine the dephasing time of two-photon photoemission and the intrinsic lifetime of bright excitons at the WSe 2 K-valley.
Model construction For the derivation of the equations for the 3 level system, we use the formula described by Martin Weinelt's. 10 To simplify the formalism, the energy variables are expressed in terms of "detunings", i.e.: In the application of OBE to photoemission, ∆ a and ∆ b are used with a different conceptual approach. As the pump creates transitions between states in the solid, ∆ a = 0 is intended to arise either from wavelength detuning, or from energy gap detuning between Bloch states in the solid (for example, a change in the bandgap due to temperature variation). ∆ b refers, instead, to vacuum states whose energy can vary continuously. Because of this, ∆ b = 0 is used, in the following, as a way to map the measured ARPES spectra in the final detected state: ω b and E f are considered fixed, and ∆ b = 0 allows to produce transitions from states energetically above or below E n to the observed final state. The energy-resolved spectrum can therefore be obtained by continuously varying ∆ b .
The derivation of OBE is based on the Liouville-Von Neumann equation: The equation describes the evolution of the density operator under a perturbationV with unperturbed HamiltonianĤ 0 . Note that has been set to 1, as will be in the following. The perturbing fields are described in the dipole approximation by the quantities: where the temporal envelope functions E a,b (t) of both pump and probe pulses are assumed to be Gaussian distributions and D is the dipole operator. In order to account for intrinsic population decay and decoherence term in eq.2, we add the damping operatorΓ, that can be represented in matrix form as:Γ The terms Γ j and Γ * j describe the decay and the dephasing rates of state j, respectively. In eq.5, the intrinsic lifetime of excitons can be characterised by the diagonal termΓ n , describing the population decay rate. The initial and the final state are assumed to have an infinite lifetime, i.e., Γ i = Γ f = 0. On the other hand, the off-diagonal terms of intermediate state contains half of the decay rate and the complex part Γ * n + Γ * i , which is the so-called "pure-dephasing rate", describing the decay of quantum coherence between the levels n and i.
By substituting eq.3-5 in eq.2, and applying the rotating wave approximation (neglecting high frequency oscillating terms), we obtain the following optical Bloch equations: where we have set: This is the set of equations needed to describe a three level system. To adapt the formalism to include K-Σ incoherent scattering, we need to create an additional two level system coupled to the main three level system via the K decay rate Γ n . We have to build an equation for the Σ state (subscript s), one for the photoemitted final state (subscript g), and one for the corresponding coherence (subscript sg). To meet the experimental reality, we also insert a backscattering parameter B, allowing to account for all mechanisms transferring population from Σ to K (eq. 7 has to be modified, see eq. 13). For a schematic of the 5-level system, see Fig. S.4(a).
nf − Γ n ρ nn +Γ n Bρ ss (13) Here, we define ρ (2) sg = e −iω b t ρ sg . These equations constitute the model used for the numerical fit of the data.
Assumptions In order to simplify the fitting procedure, we assume that the initial state dephasing time is equal to the excited states 11 , i.e., Γ * i = Γ * n = Γ * s , and that the final states dephasing time is infinite (Γ * f = Γ * g = 0). In this way, the only fitting parameters of the Gamma matrix are the inverse lifetime Γ n and the pure dephasing rate Γ * n of the bright excitonic state at K, together with the backscattering fraction B. However, a large number of parameters relative to the optical excitation require further optimization, in particular the determination of the time of coincidence of the pump and probe pulses. This quantity must be defined with a precision much higher than the measurable cross-correlation FWHM in order to achieve good confidence in the fit and determine lifetimes shorter than the duration of the pulses.
Given the aforementioned approximations, reducing the band structure to a system of non-dispersing levels, a ploy may be used to achieve this. As mentioned in the main text, the density matrix treatment of the photoemission process upon resonant pumping allows correctly treating the mixing of the two types of quantum processes to reach the final state: the first one consists of a coherent two-photon process with an intermediate virtual state, and the second one involves a two-step process, where the A exciton is created at the K valley and then photoemitted to the final state by the probe photon. The former can be separated from the latter by wavelength detuning of the pump to values smaller than the band gap (∆ a < 0). If the detuning becomes very large (∆ a << 0) the intermediate states of the two processes are very far apart in energy and they can be observed separately by solving the dynamics at two different ∆ b . At probe detuning ∆ b = 0 the dynamics is only determined by the two-step process involving the formation of a real population in the K valley; while at probe detuning ∆ b = −∆ a the final state population arises from the coherent two-photon process with a detuned virtual state. The evolution in time of the latter allows to precisely assess the coincidence between the pump and probe pulses.
To achieve a good fit, it is therefore advantageous to evaluate ρ f f once with a very large detuning (which we may define as ∆ a ) and extract its time evolution at ∆ b = −∆ a to isolate the coherent process dynamics and accurately determine the pump-probe coincidence. We can further evaluate ρ f f a second time at a pump detuning equal to the experimental pump detuning (which we define as ∆ a ), and extract its time evolution at ∆ b = 0 to determine the population dynamics of the K-exciton state. A third curve, arising from ρ gg evaluated at probe detuning ∆ b = 0 and pump detuning ∆ a , completes the set of three curves that can be extracted from the model and fitted to the experimental data, with the result plotted in Fig.3(a) of the main text.
The effective coincidence time is the fourth (and last) fitting parameter.   The latter procedure also allows observing the cross correlation between parameters and results of the fit, as displayed in Fig. S

Self-energy calculation
For the density functional electron-phonon calculation we used the EPW code 13 as part of the quantum espresso package 14 using norm-conserving pseudopotentials and the local density approximation. The hexagonal structure of WSe 2 was relaxed with lattice constants of a=6.2020 a.u, c=3.9548 a.u. and the density was computed with the pw.x program with a planewave cutoff of 160 Ry and a Monkhorst-Pack-grid of 42x42x16 k-points. The phonons were evaluated using the ph.x program on a q-mesh of 6x6x2. Based on a wannierization of the electronic structure using the wannier90.x program, 15 the EPW code interpolates the phonon and electron-phonon coupling quantities to a 20x20x10 q-mesh to evaluate the self energy, shown in the main text.  The electron-phonon lifetime is evaluated at each band n and k-point k as, where the sum runs over all phonon modes (q, ν) and electronic bands . g nmν are the electronphonon coupling matrix elements and f and n are the Fermi and Bose distributions, respectively.

Microscopic calculation of exciton signatures in trARPES
The microscopic description of the time-and angle-resolved photoemission spectroscopy including the Coulomb interaction between electrons and holes and the exciton-phonon scattering dynamics is based on a many-particle Hamiltonian and the Heisenberg equation of motion formalism. We treat the quasi-particle band structure within a parabolic approximation, with effective masses stemming from ab initio calculation established in literature, 16 and include the electron-light interaction for the VIS pump and the XUV probe pulse semi-classically in length gauge. In order to take into account the electron-hole Coulomb interaction we exploit the unit operator method exploiting the completeness relation of the Fock space. 17,18 Hence, the conduction band electron operators are expressed uniquely by electron-hole pair excitations, which are transformed to the exciton picture by introducing exciton relative and center-of-mass momentum. The eigenenergies and wave functions of the excitons are computed by numerically solving the Wannier equation.
with the index ξ v/c standing for the valley of the involved valence or conduction band electron, the reduced mass m ξvξc and the exciton state µ. The screening of the Coulomb potential V q , due to the dielectric environment, is treated beyond the Rytova-Keldysh framework. [19][20][21] To model the bulk we assume that the XUV pulse irradiates only the first two layers of the crystal.
According to this, we take a bilayer band structure from ab initio calculations 22 restricted to the 1s-A-exciton. The first term accounts for photoemission of valence band electrons P vf k ,k = v †ξv k f k with the corresponding Rabi-frequency Ω vf k . The second and third terms describe the exciton-assisted photoemission yielding the excitonic signals in trARPES, with the corresponding Rabi-frequencies Ω cf ξc− k = Ω cf ξc k ϕ * ξcξv k exp(iω vis t) and Ω cf ξc k,Q = Ω cf ξc k ϕ * ξcξv k −α ξc ξv Q . The exciton wavefunction is denoted by ϕ * ξcξv k −α ξc ξv Q with center-of-mass momentum Q and mass ratio α ξc ξv = m e /(m e + m h ). The second term couples the excitonic coherence P * ξcξv 0 , excited by the VIS pump pulse, with the photoemission transition of valence band electrons. In contrast, the third term is driven by the excitonic population caused by the phonon-induced decay of the excitonic coherence. [24][25][26] We find a T * 2 time of 18.1 fs. The exciton dynamics describe the phonon-induced relaxation and thermalization throughout the Brillouin zone expressed by a Boltzmann-like scattering equation. 27 For the exciton-phonon interaction the underlying electron-phonon deformation potential matrix elements for acoustic and optical phonons are taken from ab initio calculations. 28 In a first approximation, we use for the bilayer the same electron-phonon matrix elements as for the monolayer, which is supported by similar temperature dependent shifts in absorption experiments when going from mono-to bilayer. 26 The trARPES signal for small pump-probe delays is determined by the second term. The excitonic coherence oscillates with the excitation energy and decays with increasing time forming incoherent excitons, lying at the exciton energy, which explains an observable time dependent shift of the excitonic trARPES signal for detuned excitations. At large pump-probe delays the third term is the origin of the excitonic trARPES signal.
The slight deviation of the MDC between theory and experiment might be explainable by the spectrally broad pump pulse used in the experiment, which generates a hot exciton distribution of KK excitons (electron and hole situated at the K-point), which broadens the momentum distribution curve even at τ = 0.