Sample preparation<\/p>\n
To fabricate the WSe2\/PTCDA heterostructure (Extended Data Fig. 1a<\/a>), hBN was first exfoliated onto a 0.1% niobium-doped SrTiO3 substrate and an approximately 50-nm-thick flake was identified by optical microscopy. In a parallel procedure, WSe2 monolayers were directly exfoliated onto a silicone gel film (DGL Film, Gel-Pak) and identified through optical microscopy and Raman spectroscopy. Afterwards, a monolayer WSe2 flake was transferred from the silicone gel film onto the hBN flake on the SrTiO3 substrate (Extended Data Fig. 1b<\/a>). After introduction into ultrahigh vacuum (<5\u2009\u00d7\u200910\u22129\u2009mbar), the sample was annealed at 670\u2009K for 2\u2009h to ensure a clean sample surface. The bare monolayer WSe2 was analysed with the momentum microscope in real space (Extended Data Fig. 1c<\/a>) and reciprocal space (Extended Data Fig. 2e<\/a>), showing the expected characteristic features of monolayer WSe2, that is, the spin-split valence bands at the K valley and a single parabolic band at the \u0393 valley below the global VBM at the K valley (Fig. 1b<\/a> and refs. 20<\/a>,23<\/a>). Moreover, the clear separation between the top and the bottom valence band (Fig. 1c<\/a>) indicates the high-quality of the WSe2 with only contributions of inhomogeneous broadening.<\/p>\n Subsequently, approximately a monolayer of PTCDA was thermally evaporated onto the sample, which was maintained at room temperature (base pressure <1\u2009\u00d7\u200910\u22129\u2009mbar). The deposition rate was monitored with a with a quartz crystal microbalance and calibrated using the known deposition of PTCDA onto a Ag(110) crystal surface. On the Ag(110) surface, the first monolayer of PTCDA is adsorbed in a brickwall structure whereas additional layers grow in a Herringbone structure with a different superstructure that can be analysed by low energy electron diffraction (LEED)51<\/a>. By step-wise evaporation of PTCDA onto Ag(110) and recording of the LEED pattern, the evaporation rate was determined and was then used to deposit a monolayer PTCDA onto monolayer WSe2. The successful deposition of a monolayer PTCDA onto monolayer WSe2 was confirmed by the observation of additional spectral weight in the static ARPES data, which can be attributed to the HOMO of PTCDA (Fig. 1c,d<\/a>) and backfolded WSe2 bands (Extended Data Fig. 1e<\/a>), which are caused by the adsorbed molecular PTCDA superstructure (Fig. 1b<\/a>). The superstructure matrix<\/p>\n $$M=\\left[\\begin{array}{cc}1.58&6.78\\\\ 4.39&1.32\\end{array}\\right]$$<\/p>\n was determined by LEED on a cleaved WSe2 bulk crystal (Extended Data Fig. 1d<\/a>) and is in good agreement with ref. 28<\/a> and the ARPES data (Extended Data Fig. 1e<\/a>).<\/p>\n Femtosecond momentum microscopy<\/p>\n All photoemission data were acquired with the G\u00f6ttingen in-house photoemission set-up17<\/a>,52<\/a> that combines a time-of-flight momentum microscope53<\/a> (Surface Concept GmBH, ToF-MM) with a 500-kHz high-harmonic generation beamline (26.5\u2009eV p-polarized, 20\u2009fs). For the time-resolved measurements, the photon energy of the s-polarized pump was tuned to h\u03bd\u2009=\u20091.7\u2009eV with 50\u2009\u00b1\u20095\u2009fs pulse duration using an optical parametric amplifier. The experimental set-up is described in detail in ref. 17<\/a>. The pump fluence was adjusted to 280\u2009\u00b1\u200920\u2009\u03bcJ\u2009cm\u22122, which results approximately in an initial K-exciton density of (5.4\u2009\u00b1\u20091.0)\u2009\u00d7\u20091012\u2009cm\u22122 (refs. 31<\/a>,48<\/a>). All experiments are performed with an energy, momentum and time resolution of 200\u2009\u00b1\u200930\u2009meV, 0.04\u2009\u00b1\u20090.01\u2009\u00c5\u22121 and 54\u2009\u00b1\u20097\u2009fs (refs. 17<\/a>,31<\/a>,54<\/a>). The static measurements (Fig. 1c<\/a>) were performed at T\u2009=\u200950\u2009K, while all pump\u2013probe delay-dependent measurements were performed at room temperature (300\u2009K).<\/p>\n Photoemission data processing<\/p>\n The time-of-flight momentum microscope enables the simultaneous measurement of the kinetic energy and both in-plane momenta of the photoemitted electrons53<\/a>. However, the acquired three-dimensional photoemission data are affected by various lens aberrations and other distortions, such as pump- and probe-induced space-charge effects and surface photovoltage55<\/a>,56<\/a>,57<\/a>. Therefore, the photoemission data needs to be preprocessed before further evaluation by (1) correcting a time-dependent rigid energy shift and (2) correcting for distortions that are induced by the projection and focal lens system.<\/p>\n First, the time-dependent energy shift was corrected by minimizing the variance between the momentum-integrated spectra for E\u2009\u2212\u2009EVBM\u2009<\u20091.8\u2009eV. Second, an additional measurement was performed with a grid inserted in the Fourier plane58<\/a>,59<\/a>. We then determined the parameters for an affine transformation that maps the measured data onto an undistorted and energy-independent grid. This transformation is applied to all datasets. Small remaining distortions induced by the first lens system were corrected by fitting the positions of the K-excitons and mapping them onto an equilateral hexagon. The same positions were used to perform the momentum calibration using the lattice constant of \\({{\\rm{WSe}}}_{2} \\, {a}_{{{\\rm{WSe}}}_{2}}=0.3297\\,\\rm{nm}\\). In addition, for each delay step, the data were momentum-wise normalized to the energy range between E\u2009\u2212\u2009EVBM\u2009=\u2009\u22121.8 and \u22123.8\u2009eV. This momentum-wise normalization accounts for potential changes in illumination due to possible instabilities during the long integration times of the time-resolved measurements.<\/p>\n Quantitative analysis of the exciton energies and dynamics<\/p>\n The EUV laser pulses fragment the Coulomb-bound electron\u2013hole pairs into their single-particle components. As this process conserves energy and momentum33<\/a>,60<\/a>,61<\/a>, the exciton energies Eexc can be extracted by fitting the delay-integrated (100\u2013500\u2009fs), background-substracted (see the non-negative matrix formalism (NMF) method below) and momentum-filtered EDCs shown in Fig. 2c<\/a> with either one (K- and \u03a3-excitons) or two (hX) Gaussian peaks Ip and an exponential background Ibg, that is,<\/p>\n $${I}_{{\\rm{p}}}(E\\;)=\\frac{A}{\\sigma \\sqrt{2\\uppi }}\\exp \\left\\{\\left(-\\frac{{(E-\\mu )}^{2}}{2{\\sigma }^{2}}\\right)\\right\\},$$<\/p>\n \n (1)\n <\/p>\n $${I}_{{\\rm{bg}}}(E\\;)={A}_{{\\rm{bg}}}\\exp \\left\\{\\left(-\\frac{E}{{\\tau }_{{\\rm{bg}}}}\\right)\\right\\}.$$<\/p>\n \n (2)\n <\/p>\n The extracted peak energies E\u2009\u2212\u2009EVBM of the K and \u03a3 excitons directly correspond to the exciton energies Eexc because the hole resides at the VBM of WSe2. For the hX, the peak energy of the higher lying peak at E\u2009\u2212\u2009EVBM\u2009=\u20091.57\u2009\u00b1\u20090.05\u2009eV directly corresponds to \\({{{E}}}_{{\\rm{exc}}}^{{\\rm{hX}}}\\), whereas the lower-energy peak at E\u2009\u2212\u2009EVBM\u2009=\u20090.38\u2009\u00b1\u20090.05\u2009eV has to be referenced to the HOMO at E\u2009\u2212\u2009EVBM\u2009=\u2009\u22121.2\u2009\u00b1\u20090.1\u2009eV, which results in the same exciton energy \\({{{E}}}_{{\\rm{exc}}}^{{\\rm{hX}}}\\)\u2009=\u20091.58\u2009\u00b1\u20090.1\u2009eV. In Extended Data Table 1<\/a>, the quantified exciton energies \\({{{E}}}_{{\\rm{exc}}}^{i}\\) are compared with the BSE@G0W0 calculations and ARPES and photoluminesence experiments on monolayer WSe2 (refs. 22<\/a>,23<\/a>,62<\/a>). The total error of the experimental values is estimated to be approximately 0.05\u2009eV, taking into account fitting errors and possible errors induced by the energy calibration and space charge effects.<\/p>\n To analyse the exciton dynamics of the K, \u03a3 and hX photoemission signatures (Fig. 5<\/a>), we filter the raw photoemission data, that is, without background subtraction, by their energy and momentum coordinate. The respective EDCs and the chosen region of interests are shown in Extended Data Fig. 4<\/a>. To quantify the rise time, we fitted the energy- and momentum-filtered time-resolved photoemission spectral weight traces with an error function<\/p>\n $$I(t)=\\frac{1}{2}\\left(1+\\,\\text{erf}\\,\\left(\\frac{t-{\\mu }_{{\\rm{onset}}}}{\\sqrt{2}{\\sigma }_{{\\rm{rise}}}}\\right)\\right)$$<\/p>\n \n (3)\n <\/p>\n in the delay regions \u2212200\u2009fs to 0\u2009fs, \u2212200\u2009fs to 20\u2009fs and \u2212200\u2009fs to 150\u2009fs, respectively (Extended Data Fig. 6<\/a>). Here, \u03bconset indicates the onset time, while \u03c3rise is directly related to the rise time.<\/p>\n Similarly, we fitted the decay of photoemission spectral weight with a bi-exponential decay (equation (4<\/a>)) between 0\u2009fs and 2,000\u2009fs and between 20\u2009fs and 2,000\u2009fs for the K- and \u03a3-excitons, respectively. The hX was fitted with a single exponential decay function (equation (5<\/a>)) between 200\u2009fs and 2,000\u2009fs:<\/p>\n $$I(t)=A\\left(\\frac{1}{1+f}\\exp \\left\\{\\left(-\\frac{t}{{\\tau }_{{\\rm{fast}}}}\\right)\\right\\}+\\frac{f}{1+f}\\left(-\\frac{x}{{\\tau }_{{\\rm{slow}}}}\\right)\\right),$$<\/p>\n \n (4)\n <\/p>\n $$I(t)=A\\exp \\left\\{\\left(-\\frac{t}{{\\tau }_{{\\rm{slow}}}}\\right)\\right\\}.$$<\/p>\n \n (5)\n <\/p>\n The relevant time constants are given in Extended Data Table 1<\/a>.<\/p>\n Non-negative matrix factorization<\/p>\n Due to the small size of the WSe2 monolayer and the small real-space selection aperture (10\u2009\u03bcm effective diameter), the measurement was susceptible to a time-independent background intensity. Therefore, the excited state momentum maps and momentum-filtered EDCs shown in Fig. 2<\/a> and in the insets of Extended Data Fig. 4<\/a> and 6<\/a> were background subtracted.<\/p>\n For the background determination, we used NMF, as implemented in the scikit-learn package for Python63<\/a>. NMF is a dimensionality-reduction method that so far remains unexplored in time-resolved ARPES, but it has recently found application in spatially resolved material science studies based on X-ray diffraction64<\/a> and also static ARPES experiments65<\/a>. Similar to principal component analysis, NMF is based on the numerical factorization of a given matrix X into two matrices W and H, with the additional condition that all matrices have only non-negative elements. In our case, X is given by the time-dependent raw dataset where we consider only excited state data above E\u2009\u2212\u2009EVBM\u2009=\u20090.15\u2009eV. In addition, we fix W by a static background and the four extracted time traces of the K, \u03a3, hX@1.6 eV and hX@0.4 eV photoemission signal plotted in Fig. 5<\/a>. The determined output then is the matrix H that consists of five components, each following one of the four given time dependencies and the time-independent background. Extended Data Fig. 3<\/a> shows extracted components integrated over the regions of interest in energy. Notably, the different components 1\u20134 can be assigned in reasonable agreement to the different excitonic photoemission signatures despite the strong overlap in time, energy and momentum (orange hexagons, grey circles and blue circle). Component 5 is time-independent and used for background substraction.<\/p>\n Photoemission orbital tomography<\/p>\n Using the plane-wave model of photoemission, the measured momentum-dependent photoemission intensity of electrons emitted from molecular orbital can be expressed as34<\/a><\/p>\n $$I({\\bf{k}})={\\left\\vert {\\bf{A}}\\cdot {\\bf{k}}\\right\\vert }^{2}{\\left\\vert {\\mathscr{F}}\\left(\\psi ({\\bf{r}})\\right)\\right\\vert }^{2}\\delta \\left({E}_{\\rm{b}}+{E}_{{\\rm{kin}}}+\\varPhi -h\\nu \\right),$$<\/p>\n \n (6)\n <\/p>\n where \u03c8(r) is the real-space electronic wavefunction, \\({\\mathscr{F}}\\) is the Fourier transform and \\({\\left\\vert {\\bf{A}}\\cdot {\\bf{k}}\\right\\vert }^{2}\\) is a polarization factor defined by the vector potential A of the incoming electromagnetic field. The Dirac \u03b4 function ensures energy conservation of the photoemission process, which includes the photon energy h\u03bd, the electron binding energy Eb, the work function \u03a6 and the kinetic energy Ekin of the emitted photoelectron. This model has been successfully applied to analyse orbital wavefunctions of transient excited states in PTCDA26<\/a> and extended to the description of the photoemission signature of excitons in C60 (ref. 36<\/a>). According to refs. 36<\/a>,37<\/a>, the photoemission signature of the hX with multiple hole contributions, but only a single electron contribution, must feature a two-peak structure, where the momentum distribution of both peaks resembles the Fourier transform of the LUMO of PTCDA as described by equation (6<\/a>). Based on this model, we calculate the expected momentum map of the HOMO and the hX considering all the different orientations of the PTCDA molecule28<\/a>. The real-space molecular orbitals calculated by DFT are extracted from ref. 66<\/a>. The results are plotted in Extended Data Fig. 2c,g<\/a>. We note that the theoretical momentum fingerprints were calculated for single-particle electrons (that is, using the Kohn\u2013Sham orbitals). In the near future, progress in the field of exciton photoemission orbital tomography36<\/a>,37<\/a> may enable the calculation of predicted momentum fingerprints also for excitonic states; however, such calculations are currently not possible for the present WSe2\/PTCDA structure.<\/p>\n Calculation of the electronic structure<\/p>\n A G0W0 treatment of the herringbone-type WSe2\/PTCDA heterostructure is beyond current computational possibilities. Instead, to meet the experimental conditions as closely as possible, we consider a configuration of a PTCDA molecule adsorbed on a 4\u2009\u00d7\u20094\u2009\u00d7\u20091 supercell of the pristine WSe2 structure with an in-plane lattice parameter of 3.317\u2009\u00c5 (Extended Data Fig. 5<\/a>; for a detailed discussion of the implemented supercell, refer to the Supplementary Information<\/a>). We optimize the atomic structure, consisting of 86 atoms, using the all-electron code \u2018FHI-aims\u201967<\/a> by minimizing the amplitude of the interatomic forces below a threshold value of 10\u22123\u2009eV\u2009\u00c5\u22121. For all species a tight basis is used. The resulting adsorption geometry is shown in Extended Data Fig. 5a<\/a>. The PTCDA molecule is slightly tilted, with the shortest and longest distance to the substrate being 2.87\u2009\u00c5 and 4.98\u2009\u00c5, respectively, measured from the top of the substrate.<\/p>\n The ground-state, G0W0 and BSE calculations are performed using the all electron full-potential code \u2018exciting\u201929<\/a>, which implements the family of linearized augmented plane wave plus local orbitals (LAPW+LO) methods. The muffin-tin spheres of the inorganic component are chosen to have equal radii of 2.2\u2009bohr. For PTCDA, the radii are 0.9\u2009bohr for hydrogen (H), 1.1 bohr for carbon (C) and 1.2 for oxygen (O). The electronic properties are calculated first using DFT with the generalized gradient approximation in the Perdew\u2013Burke\u2013Ernzerhof parametrization for the exchange\u2013correlation (xc) functional. The sampling of the BZ is carried out with a homogeneous 3\u2009\u00d7\u20093\u2009\u00d7\u20091 Monkhorst\u2013Pack k-point grid. To account for van der Waals forces and intermolecular interactions, we adopt the Tkatchenko\u2013Scheffler method68<\/a>. The quasi-particle (QP) energies are computed within the G0W0 approximation69<\/a>, where we include 200 empty states to compute the frequency-dependent dielectric screening within the random-phase approximation. A 2D truncation of the Coulomb potential in the out-of-plane direction z is used70<\/a>. The band structure is computed by using interpolation with maximally localized Wannier functions71<\/a> and Fourier interpolation (Extended Data Fig. 5b<\/a> and Extended Data Fig. 5c<\/a>, respectively). To keep the calculations feasible, SOC is not considered in this work. Although SOC leads to a splitting of the lowest-energy excitonic peak by approximately 450\u2009meV (ref. 72<\/a>), it would not alter the type-I level alignment of WSe2\/PTCDA. Importantly, the molecular states involved in the formation of the hX, that is, the HOMO and LUMO, would not be affected by the inclusion of SOC.<\/p>\n To allow a direct comparison with the experimental ARPES data (Fig. 1d<\/a>), we unfold the theoretical band structure by symmetry mapping of the Bloch-vector-dependent quantities defined in the supercell into the unit-cell calculations. Here, the wavefunctions are constructed in a uniform real-space grid of 120\u2009\u00d7\u2009120\u2009\u00d7\u2009120 and used to calculate the spectral function (Fig. 1f<\/a>).<\/p>\n The QP band gap of WSe2 in the heterostructure is in good agreement with that measured by STS28<\/a>; however, the PTCDA gap is underestimated, which is most evident in the level alignment of the HOMO. This discrepancy can be explained by the interplay of different effects such as SOC, the choice of the xc functional and its role as a starting point for the QP calculations, and the interlayer distance between PTCDA and WSe2. Uncertainties in the latter can be related to packing density, the xc functional and the treatment of van der Waals interactions, or temperature73<\/a>. In earlier studies on ZnO\/WSe2 it has been shown that increasing the interlayer distance leads to a noticeable increase in the QP gaps on both sides of the interface74<\/a>. Also in WSe2\/PTCDA, increasing (decreasing) the interlayer distance will decrease (increase) the mutual screening, leading to an increase (decrease) in the HOMO\u2013LUMO gap. This, in turn, would lead to an increase (decrease) in the VBM\u2013HOMO distance. Overall, there is an interplay of effects on the order of a few tenths of an electronvolt each, which can only be resolved through extensive future QP calculations.<\/p>\n To overcome this mismatch, we apply a scissors shift to the molecular levels. Shifting the LUMO by \u221250\u2009meV closer to the experimental value28<\/a> results in very good agreement of the excitonic spectrum with experiment (see below).<\/p>\n Calculation of the exciton spectrum<\/p>\n For the calculation of the exciton spectrum, we solve the BSE on top of the QP band structure, where the screened Coulomb potential is computed using 100 empty bands. In the construction and diagonalization of the BSE Hamiltonian, 16 occupied and 14 unoccupied bands are included, and a 12\u2009\u00d7\u200912\u2009\u00d7\u20091 shifted k-point mesh is adopted. Calculations are performed using the BSE module75<\/a> of the \u2018exciting\u2019 code.<\/p>\n Calculation of the correlation function<\/p>\n Following the definition in ref. 40<\/a>, we calculate the electron\u2013hole correlation function<\/p>\n $${F}^{i}(\\bf{r})={\\int}_{\\varOmega }{d}^{\\,3}{{\\bf{r}}}_{\\rm{e}}| {\\psi }_{i}({{\\bf{r}}}_{\\rm{h}}={{\\bf{r}}}_{\\rm{e}}+{\\bf{r}},{{\\bf{r}}}_{\\rm{e}}){| }^{2},$$<\/p>\n \n (7)\n <\/p>\n where Fi describes the probability of finding electron and hole separated by the vector r\u2009=\u2009rh\u2009\u2212\u2009re. We approximate this integral by a discrete sum over a finite number of fixed electron coordinates. For each electron position, the hole probability \u2223\u03c8i(rh\u2009=\u2009re\u2009+\u2009r, re)\u22232 is computed on an evenly spaced, dense grid of 100\u2009\u00d7\u2009100\u2009\u00d7\u2009100 sampling points, covering approximately 3\u2009\u00d7\u20093\u2009\u00d7\u20091 supercells. For the hX, we sampled 60 positions on the PTCDA molecule (0.5\u2009\u00c5\u22121 below and above the carbon and oxygen atoms) because its electronic contribution is almost entirely composed of the LUMO of PTCDA. Similarly, we calculated the electron\u2013hole correlation function of the K-exciton (Extended Data Fig. 7<\/a>), which is completely localized in the WSe2 layer. Here, we sampled 16 positions close to the W atoms where we expect a high probability of finding the electron.<\/p>\n For further analysis, the three-dimensional correlation function is split into its in-plane and out-of-plane components by integrating over the other direction (Fig. 4<\/a> and Extended Data Fig. 7<\/a>). Here, the intralayer (purple) and interlayer (green) in-plane distributions were extracted by integrating exclusively over the respective peak of F(r\u22a5), that is, r\u22a5\u2009=\u2009\u22123 to 4.5\u2009\u00c5 and r\u22a5\u2009=\u2009\u221210.5 to \u22123\u2009\u00c5 for the intra- and interlayer components, respectively. Notably, the in-plane component shows a distinct periodic pattern (see insets in Fig. 4b<\/a> and Extended Data Fig. 7b<\/a>). Thus, to extract the in-plane radial profile and the root mean square (RMS) radius, we first filter the data in Fourier space, thereby smoothing it in the real space.<\/p>\n Spatial analysis of the K-exciton and comparison with hX<\/p>\n In analogy to the analysis of the spatial structure of the hX in Fig. 4<\/a>, we analyse the K-exciton wavefunction (Extended Data Fig. 7<\/a>). From the BSE calculation, we find that FK(r\u22a5) is dominated by a single peaked feature centred around r\u22a5\u2009\u2248\u20090\u2009\u00c5, implying that the exciton is of pure intralayer character (Extended Data Fig. 7a<\/a>). Consequently, the probability density of finding K-excitons in the WSe2 layer is nearly 100%. Consistent with this, the two smaller side peaks located at a distance corresponding to the distance between the tungsten and selenium planes dW,Se, can be attributed to a residual probability of the electron and\/or hole being at the selenium atoms.<\/p>\n Due to its hydrogen-like structure, the in-plane electron\u2013hole probability distribution of the K-exciton can be directly reconstructed from the experimental photoemission momentum fingerprint via Fourier analysis31<\/a>,41<\/a>,42<\/a>,59<\/a>. This allows a direct comparison of FK(r\u2225) between theory and experiment (Extended Data Fig. 7b<\/a>). Both theory and experiment confirm the pure Wannier-like character of the K-exciton because the radial distribution is much larger than the WSe2 lattice constant. For a more quantitative analysis, we compare the extracted RMS radii to be (\\({{{r}}}_{{\\rm{K}}}^{{\\rm{BSE}}}\\)\u2009=\u200914\u2009\u00c5) (theory) and \\({{{r}}}_{{\\rm{K}}}^{\\exp }\\)\u2009=\u200910\u2009\u00b1\u20091\u2009\u00c5 (experiment), which are in excellent agreement. Note that, due to the finite momentum resolution of the photoemission signal, the derived RMS radius represents a lower limit of the true value. The distinct intralayer and Wannier-like character of the K-exciton can be further visualized by plotting the isosurface of a representative fixed electron and hole position (Extended Data Fig. 7c<\/a>), which stands in clear contrast to the isosurfaces of the hX (Fig. 4<\/a>).<\/p>\n After having identified the K-exciton as a Wannier exciton, we compare its in-plane correlation function directly to that of the inter- and intralayer contributions to the hX (Extended Data Fig. 8<\/a>). Here, we find that the interlayer contribution resembles the spatial distribution of the K-exciton, while the intralayer contribution is more localized and exhibits a stronger spatial modulation stemming from the molecular orbitals. This direct comparison confirms the previously assigned Frenkel- or Wannier-type character of the intra- and interlayer contribution of the hX.<\/p>\n