Temperature-dependent electronic ground state charge transfer in van der Waals heterostructures

Electronic charge rearrangement between components of a heterostructure is the fundamental principle to reach the electronic ground state. It is acknowledged that the density of states distribution of the components governs the amount of charge transfer, but a notable dependence on temperature has not yet been considered, particularly for weakly interacting systems. Here, we experimentally observe that the amount of ground state charge transfer in a van der Waals heterostructure formed by monolayer MoS2 sandwiched between graphite and a molecular electron acceptor layer increases by a factor of three when going from 7 K to room temperature. State-of-the-art electronic structure calculations of the full heterostructure that account for nuclear thermal fluctuations reveal intra-component electron-phonon coupling and inter-component electronic coupling as the key factors determining the amount of charge transfer. This conclusion is rationalized by a model applicable to multi-component van der Waals heterostructures.


Introduction
The discovery that atomically thin layers can exist at room temperature in air marked the launch of research on two-dimensional (2D) materials. [1,2] A single layer of a 2D material may be insulating (e.g., hexagonal boron nitride), semiconducting (e.g., MoS2), or conducting (e.g., graphene), and thus all electrical material properties required for the construction of an electronic or optoelectronic device are available in the monolayer limit. [3][4][5][6] Stacks of such monolayers are typically bound by weak van der Waals (vdW) interlayer interactions, and vdW heterostructures have emerged as prime candidates for realizing electronic and optoelectronic functions in the smallest possible volume. The type and efficiency of the functionality that can be achieved depends critically on the electronic energy level alignment across the heterostructure, [7] and substantial charge density re-arrangement upon contact can occur, depending on the electronic structure of each component. Consequently, charge rearrangement and also charge transfer (CT) phenomena that define the electronic ground state of vdW heterostructures must be thoroughly understood for knowledge-guided device design. Beyond layered inorganic 2D materials, organic molecular semiconductors are also attractive as a component in vdW heterostructures, because they extend the range of available energy gap values and feature strong light-matter interaction. [8,9] Furthermore, it has been recognized that strong molecular electron acceptors and donors can be employed as dopants for numerous 2D semiconductors, particularly for transition metal dichalcogenides (TMDCs). [10][11][12][13][14] This doping, on the one hand, allows controlling the Fermi level (EF) position and mobile carrier density in the semiconductor, and, on the other hand, enables manipulation of optical quasi-particles, such as charged excitons (positive or negative trions) in TMDC monolayers. [15][16][17][18][19] Such vdW heterostructures can be ideal model systems for exploring intriguing physical phenomena, [20,21] and they can pave the way for a broader class of device applications. Despite many opportunities offered by vdW heterostructures, the current understanding of contactinduced charge density re-arrangement phenomena and inter-layer ground state CT is not well established for such weakly bound stacks, mostly due to their challenging complexity. For instance, organic semiconductors feature vast spatial degrees of freedom, and 2D semiconductors experience substantial band structure renormalization as function of the thermodynamic conditions [22,23] and of the surrounding environment, including a supporting substrate. [11] A recent study on vdW heterostructures comprising the molecular electron acceptor 1,3,4,5,7,8-hexafluoro-tetra-cyano-naphthoquinodimethane (F6TCNNQ) deposited on monolayer (ML) MoS2 revealed the key role of the substrate. [11] When insulating sapphire was used as substrate for ML-MoS2, electrons were transferred from n-type gap states (induced by native sulfur-vacancies) to molecular acceptors. In contrast, when sapphire was replaced by highly oriented pyrolytic graphite (HOPG), the n-type gap states of MoS2 were emptied into the charge reservoir of the conductive substrate and the electron transfer took place directly from HOPG to F6TCNNQ, as schematically shown in Figure 1. This scenario was discussed in analogy to a parallel 4 plate capacitor, where positive and negative charges reside on the two "plates", HOPG and F6TCNNQ, respectively. The ML-MoS2 then represents a dielectric layer experiencing the electric field between the plates, and in this model the frontier energy levels of the semiconductor are not involved in CT.
Considering only this picture, the amount of charge transferred across the vdW heterostructure to reach electronic equilibrium is determined by the density of states (DOS) distribution near EF in HOPG and that of the lowest unoccupied molecular orbital (LUMO) level manifold in F6TCNNQ, as well as their energy offset. However, conjugated molecules exhibit strong electron-vibrational coupling, i.e., the electronic energy levels change upon intramolecular bond length alterations. This translates into a variation of the ground state electronic levels' energy when vibrational modes are populated. [24][25][26] Therefore, the temperature-dependent population of vibrational modes is expected to alter significantly the electron distribution near EF and hence the magnitude of CT in vdW heterostructures. This dependence of CT on temperature has not been investigated in vdW heterostructures to date. An indication of the relevance of temperature on CT is the observation that the electrical characteristics of ML-MoS2-based field-effect transistors depend notably on temperature when a molecular acceptor is deposited onto the semiconductor. [14] For molecular acceptors from the tetracarboxylic-dianhydride family deposited on silver single crystals it was found that decreasing the temperature significantly increased the amount of CT, due to the shortening of the molecule-metal distance at lower temperature. [27,28] However, these systems are characterized by covalent interactions and pronounced orbital hybridization, which makes a direct relation to vdW heterostructures uncertain. In fact, for the vdW heterostructure F6TCNNQ/HOPG we find that the amount of CT is virtually independent of temperature (see Supplementary Figure S3).
Here, we demonstrate with angle-resolved photoelectron spectroscopy (ARPES) that the amount of ground state CT in the prototypical vdW heterostructure F6TCNNQ/ML-MoS2/HOPG increases with increasing temperature, gradually and reversibly. This unexpected phenomenon is rationalized via high level electronic structure calculations of the full heterostructure, which account for electron-phonon coupling in the adiabatic limit through sampling of stochastic phonon displacements at different temperatures. The calculations show a shift to lower energies and broadening of the acceptor's LUMO manifold with increasing temperature that enhances CT. Within a simple model, it is shown that even in the absence of large-scale structural changes, the energy-gap renormalization of MoS2 and F6TCNNQ with temperature can contribute to the CT enhancement. Therefore, the earlier suggested capacitor model falls short of correctly describing such systems, because it only considers electrostatic interaction.

Results
Temperature-dependent charge transfer in the F6TCNNQ/ML-MoS2/HOPG vdW heterostructure ARPES allows tracking the valence band structure of ML-MoS2 and the characteristic valence fingerprints of negatively charged F6TCNNQ molecules (anions). [29] For a vdW heterostructure consisting of ca. 0.5 ML F6TCNNQ on ML-MoS2/HOPG the corresponding energy distribution curves (EDCs) in the near-EF region from ARPES measurements are shown in Figure 2a, for different temperatures. At room temperature (300 K), this energy region is dominated by two features centered at 0.2 eV and 0.8 eV binding energy (see Supplementary Figure S5), as the valence band onset of MoS2 is at a binding energy higher than 1.5 eV (see Supplementary Figure S6). These can readily be assigned to emission from the singly occupied former LUMO level of the F6TCNNQ anion (denoted as L*) and its relaxed highest occupied molecular orbital (HOMO) level (denoted as H*), due to electron transfer from HOPG. [11] Spectra recorded at lower temperature exhibit the same features, but with incrementally reduced intensity. We note that the same molecular anion emission features can be identified for F6TCNNQ deposited directly onto HOPG, as shown in Supplementary Figure S3, but no intensity variation as function of temperature was observed. The acceptor molecules adopt a face-on lying orientation on MoS2 supported by HOPG from room temperature to 77 K, as determined from X-ray absorption measurements (see Supplementary Figure S7), and it is reasonable to assume that this configuration persists even at lower temperatures. Consequently, the area of the emission features L* and H* is representative of the fraction of molecular anions on the surface. Plotting this area x, normalized to the signal at 300 K, as function of temperature (Figure 2b), reveals a gradual decrease in the range between room temperature and 7 K. This is a direct indication that the fraction of charged F6TCNNQ molecules at 7 K is only about one third of that found at 300 K. The fraction of charged 6 molecules at 300K was calculated from the transferred electron density divided by the surface molecule density of F6TCNNQ and found to be 21.7 %, while that at 7 K was 7.2 % (for details see Supplementary Section 3). Note that the temperature-dependent change of x is fully reversible, as the temperature sequence of this experiment started at 300 K, followed by cooling to 7 K, and subsequent warming-up steps until reaching 300 K again. To provide further support for this finding of smaller CT amount at lower temperature, we analyze the corresponding MoS2 valence band maximum (VBM) positions at the K and Γ points of the ML's Brillouin zone in Figure 2c   Energy distribution curves (EDCs) of MoS2/HOPG near the Fermi-level as a function of temperature, before (black) and after (colored) 0.5 ML F6TCNNQ deposition. The additional density of states after F6TCNNQ deposition is due to negatively charged acceptors (F6TCNNQ -) only. (b) Relative change in charged F6TCNNQ (ξ) as a function of temperature, where the temperature-axis also reflects the experimental sequence, i.e., starting at 300 K, cooling to 7 K, followed by incremental return to 300 K, to highlight the reversibility. (c) Valence band maximum (VBM) of MoS2 at the Γ and K point plotted as a function of temperature with/without absorption of F6TCNNQ. L * and H * denote the partially filled LUMO and relaxed HOMO of charged F6TCNNQ. [11] These observations are consistent with the overall CT mechanism shown in Figure 1, i.e., electron transfer from HOPG to F6TCNNQ, whereas ML-MoS2 experiences the electric field created by the separated charges and its energy level shift accordingly. The electric field is directly proportional to the amount of CT in the capacitor model, and the same temperature-dependent trend for x and the VBM shifts is expected, and verified here by comparison of Figure 2b and 2c. For the vdW heterostructure investigated here, the position of EF in the semiconductor gap can be tuned by ca. 260 meV simply by varying the temperature, as the different CT amount changes the electric field from approximately 0.26 V/nm at 7 K to 0.57 V/nm at 300 K (estimated with an effective distance between HOPG and the molecular layer of 9 Å). The origin of the temperature-dependent CT amount, however, remains to be elucidated. If lower temperature was to simply decrease nuclear motion and thus reduce the average vertical inter-atomic distances in the heterostructure, the CT amount would be expected to increase with decreasing temperature in the capacitor model; the opposite is observed in experiment. Therefore, to shed light on this observation, elaborate first-principles modeling is employed to obtain more insight into how temperature influences the density of states (DOS) of F6TCNNQ/ML-MoS2/HOPG.

Charge transfer in vdW heterostructures
We searched for stable structures of the vdW heterostructure composed of either one (dilute case) or two F6TCNNQ molecules on a 4×8 ML-MoS2 unit cell, which corresponds to almost the same molecular density as in experiments (see SI Section 3). We employed density-functional theory with a range-separated hybrid functional (see computational details in Methods). Globally, the acceptor molecules adopt a flat-lying orientation on ML-MoS2 with a small corrugation of adsorption energy (maximum 150 meV) with respect to lateral displacements (see Supplementary Figure S8 and Table   S2). The most stable structures are displayed in Figure 3a Figure S10). We note that the orbital shows non-vanishing hybridization with ML-MoS2. This is a direct manifestation of the finite electronic coupling between the two components, even if their interaction can be classified as van der Waals type.
In the dilute case, the LUMO is partially occupied by 0.3 e, whereas the partially occupied LUMO orbitals of the two non-equivalent molecules in the densely-packed case are not equally charged, i.e., we find that one molecule presents a partial occupation of 0.25 e on its partially filled LUMO and the other 0.1 e. The fractional charge transfer observed here (including the disproportional CT for two molecules in the unit cell) obtained with a range-separated hybrid functional, in contrast to the experimental evidence of integer charge transfer, is likely a consequence of the remaining electron delocalization error present in this functional, allied to the need of much larger supercells that are not computationally feasible to simulate integer CT. [30,31] With this static description (i.e., without temperature-induced effects from nuclear motion) of the vdW heterostructure we obtain adequate qualitative agreement between experiment and theory. In the following, we continue by including quantum zero-point motion and temperature-dependent effects in our calculations.

Effect of temperature on charge transfer in the vdW heterostructure
We use the total electronic DOS (TDOS) and the PDOS to analyze the possible impact of finitetemperature nuclear motion on the CT of our vdW heterostructure, for the dilute case. We pick this structure because it is sufficient to grasp the CT, as shown in the previous section. The two aspects that can change the position of electronic states, namely thermal population of vibrational modes and thermal lattice expansion, are investigated separately here. In particular, thermal fluctuations are captured through thermodynamic averages on stochastically sampled nuclear configurations consistent with the distribution of quantum harmonic phonons at different temperatures [32,33] (details given in 10 The effect of thermal lattice expansion, which is comparably pronounced for MoS2, [34,35] was investigated with static calculations. The results shown in Supplementary Figure S11 show that the effect of the expansion on the electronic orbitals is very small, but nevertheless shifts the partially occupied LUMO to lower energies with increasing temperature. Consequently, nuclear fluctuations and lattice expansion work hand in hand to promote CT in the vdW heterostructure with increasing temperature, as found in our experiment.

Discussion
Considering the above results for our specific vdW heterostructure consisting of the molecular electron acceptor F6TCNNQ, the 2D semiconductor MoS2, and the conductor HOPG, one reason for the temperature-dependent amount of charge transfer can be rationalized as consequence of the electronphonon coupling that is particularly pronounced for ML-MoS2 and the molecules. The molecules' LUMO level, which receives electrons from HOPG and slightly hybridizes with MoS2, is shifted to lower energies and broadened with increasing temperature, resulting in a larger fraction of molecules available for CT. In addition, it is well known that increasing temperature also induces a band-gap renormalization for MoS2 [36,37] (see Supplementary Figure S12 and S13). These phenomena, however, cannot be captured by the capacitor model discussed in the introduction, because the amount of CT is

Conclusions
The 14 Methods Sample preparation. ML-MoS2 were grown on sapphire via chemical vapor deposition (CVD) and transferred onto an HOPG substrate by conventional poly(methyl methacrylate) (PMMA) transfer method. [38][39][40] Sample cleaning step of 300-350 °C in situ annealing in preparation chamber (base pressure of 10 -10 mbar) was done to get rid of unwanted carbon based contamination including residual PMMA.
F6TCNNQ (Novaled), as molecular acceptor, was deposited onto the clean MoS2/HOPG in the same preparation chamber while monitoring the nominal mass-thickness by using a quartz crystal microbalance.
Photoemission measurements. ARPES spectra were measured at the beamline BL7U (UVSOR, Japan), the IMS (MBS-A1 lab system, Japan) and the Humboldt University (lab system, Germany) using a hemispherical electron analyzer with monochromatic light source of 21 eV and 9 eV, respectively. The determination of total resolution of 100 meV at 300 K were performed using clean Au (111) single crystal. The energy calibration with respect to Fermi level were carried out by measuring the electrically grounded clean Au (111) single crystal.
Theoretical modeling. The calculations presented here were based on density functional theory (DFT), [41,42] as implemented in the all-electron, numeric atom-centered electronic-structure package FHI-aims. [43,44] The Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation [45] functional was employed for the atomic geometry relaxation, while the hybrid Heyd-Scuseria-Ernzerhof (HSE06) was used for the electronic properties to grasp the correct charge transfer and level alignment. [46] The "tight" basis sets were used for PBE calculations, while the "intermediate" settings were used for HSE06. All the calculations include the Tkatchenko-Scheffer van der Waals (TS-vdW) interactions [47] and non-self-consistent spin-orbital coupling (SOC). [48] Comparison of TS-vdW with the many body dispersion of Ref. [49] shows that the corrugation of the potential energy surface is similarly described by these methods (details in Supplementary Table S2). A vacuum region of 100 Å and a dipole correction were employed to avoid the spurious electrostatic interactions between periodic supercells. 47 A 20×20×1 k-point sampling was used for the (1×1) ML-MoS2 (0001) and graphite (001) surface to attain the correct electronic properties. Four layers of graphite (001) surface were considered to guarantee the convergence of work function, which is also reported in a previous experiment. [50] The atoms in the two bottom layers were fixed while other atoms are allowed to fully relax. The minimum force on atoms is guaranteed to be below 0.001 eV/Å.
We initially performed a grid search to find the most stable geometry of one molecule adsorbed on freestanding MoS2 with lying-down and standing-up geometries at different packing densities. At the low packing density under the supercell of (6×6), 41 models for lying-down were simulated by PBE functionals, and 10 models were selected to perform the electronic property calculations. For the higher packing densities, 22 models were considered. This method is described in detail in the previous work. [25] To simulate MoS2 (0001)/graphite (001), a typical commensurate superstructure with (4×4) MoS2 and (5×5) graphite was considered, to reduce the lattice mismatch. Based on this moiré superstructure, 6 models with one molecule lying-down on the supercell graphite (5×10)/MoS2 (4×8) and 3 models with molecule lying-down, short-tilted, and long-tilted on the supercell graphite (5×5)/MoS2 (4×4) were considered to perform the geometry optimization with the PBE functional. We then used a random structure search strategy [51] to find geometries of two molecules lying flat at the (4×8) MoS2 supercells. The most stable structure was selected for consideration with the graphite substrate. In the theory part of this paper, HOPG refers to a graphite (001) (5×10) supercell. All calculations have been performed with a Gaussian occupation smearing with a width of 150 meV to ensure convergence of the self-consistent cycles. The (P)DOS shown in the manuscript have been shifted such that this occupation function assumes values below 0.09 above the zero of energy.

Section 5. Determination of orientation of F6TCNNQ on top of ML-MoS2 with temperature
Supplementary Figure S7 | Temperature dependent x-ray absorption spectroscopy (XAS).

Section 6. Calculated properties of F6TCNNQ on the free-standing ML-MoS2
Supplementary Figure    To understand the impact of temperature on the charge transfer between F6TCNNQ and HOPG alone, ARPES measurement were carried out as shown in Figure S3 and it shows only a sharpening of the F6TCNNQ features without a change of spectral weight (area). This shows that the ratio of neutral and charged F6TCNNQ molecules is not influenced by temperature.

Section 3. Estimation of the fraction of charged F6TCNNQ at 300 K.
To quantitatively estimate the amount of charge transfer, the fraction of charged F6TCNNQ molecules is calculated from the transferred electron density divided by the density of F6TCNNQ on top of the monolayer (ML) MoS2 surface (ρe/ρF6TCNNQ), based on the assumption that one electron can be transferred into only one molecule.
(1) Density of F6TCNNQ on top of ML-MoS2 surface In the F6TCNNQ deposition step, the nominal thickness was monitored by a quartz crystal microbalance (QCM), which enables us to calculate the density of F6TCNNQ on top of ML-MoS2 using the following equation: From that, we now can estimate the fraction of charged F6TCNNQ at 300 K to be 21.71 % (ρe/ρF6TCNNQ), and that at 7 K to be 7.24 % in the proportional way using ξ at 7 K as defined in the main manuscript text (ξ@7K ·ρe/ρF6TCNNQ).
Supplementary Table S1 | Impact of temperature on the amount of charge transfer and the fraction of charged molecules, obtained according to the procedure described above.  In order to address the effect of temperature on the band structure of ML-MoS2 without F6TCNNQ, ARPES spectra of clean ML-MoS2/HOPG were measured with increasing temperature from 7 K to 300 K in Figure S6 left dashed line frame. Figure S6a and b show that ARPES spectra of ML-MoS2 around the Γ point of BZ for the selected temperatures. The local valence band maximum (VBM) at Γ point is estimated by 1.94 eV, which is in good agreement with previous reported values. [1,2] In Figure S6e and f, we found two split bands with local valley (K-K') around K point due to the lack of inversion symmetry and spin-orbit coupling. [3,4] The overlap of the bands along the two high symmetry directions Γ-M and Γ-K is an intrinsic feature of 2H phase ML-MoS2 in azimuthally disordered sample. [5] It is noteworthy that the spectra of the two extreme temperatures (7 K, 300 K) are nearly identical except for the spectral broadening originating from thermal fluctuation of the atoms. From that, it can be 9 concluded that the band structure of ML-MoS2 itself does not change by thermal fluctuations within the experimental resolution.
Upon deposition of a nominal 0.5 ML of F6TCNNQ, the band structure of ML-MoS2 moves closer to the Fermi level by ca. 0.11 eV at both Γ and K points of BZ at 7 K as shown in Figure S6c and g. With increasing temperature up to 300 K, the band structure of ML-MoS2 is gradually shifted towards the Fermi level resulting in a lower binding energy of VBM of 1.72 eV at Γ and 1.54 eV at K point, respectively. Therefore, all measured spectra show a gradual and rigid shift toward the Fermi level in the whole BZ upon stepwise increase of the temperature as summarized in Figure 2c in main manuscript.
In addition, we also found a significant broadening of spectra originating from inelastic scattering of emitted photoelectrons from ML-MoS2 by the adsorbed molecules.
In short, we did not find any evidence for any change in the electronic structure of ML-MoS2 itself as a function of temperature. On the other hand, it is observed that the amount of ML-MoS2 energy shift when molecules are adsorbed depends notably on temperature.

Section 5. Determination of orientation of F6TCNNQ on top of ML-MoS2 with temperature
Supplementary Figure  In order to determine the orientation of F6TCNNQ on top of HOPG at different temperatures, x-ray absorption spectroscopy was performed at the beamline PM4 (Bessy II, Germany). First, ML-MoS2 /HOPG was annealed at 300°C in the preparation-chamber to obtain a clean surface. The XAS spectra of clean ML-MoS2 /HOPG in N-K edge energy range were measured by varying the temperature from 300 K to 71 K and the angle between incident beam and surface normal from 15° to 90°, respectively.
The spectra of ML-MoS2/HOPG measured at different angles and temperatures are identical, therefore, only two selected extreme cases are shown in Figure S7 bottom two spectra. Since there were no nitrogen atoms present, only the Mo-M2 and M3 edge features originating from ML-MoS2 are observed without temperature and angle dependency. This enables to subtract the overlapping of ML-MoS2/HOPG features easily from the spectra of F6TCNNQ/ML-MoS2/HOPG as shown in Figure S7b and explained below.
As shown in the middle two spectra in Figure S7, varying the angle between the sample and incident light significantly changes the N-K edge feature coming from the F6TCNNQ. On the other hand, the change of temperature does not impact the N-K edge feature except for its sharpness. This is a clear evidence that F6TCNNQ on ML-MoS2 adopts the same orientation with respect to the surface normal regardless of temperature.
As plotted in Figure S7b and c, we then subtract the contribution from ML-MoS2 from the XAS spectra of F6TCNNQ/ML-MoS2/HOPG in order to further analyze the N-K edge of F6TCNNQ. In Figure S7b and c, XAS spectra are further deconvoluted by one red and three violet Gaussians, according to previous reports as shown in Figure S7d (xy-plane for red peak and z-axis consisting of three blue peaks), respectively. [6] From this assignment, the orientation of F6TCNNQ is found to be lying on top of ML-MoS2 surface regardless of the temperature.

Section 6. Calculated properties of F6TCNNQ on the free-standing ML-MoS2
Supplementary Figure

Section 7. Temperature dependence of electronic levels
We here address temperature-dependence, caused by electron-phonon coupling, of the electronic density of states obtained from the single-particle Kohn-Sham levels 3 >,N calculated with the HSE06 functional. Following the discussion in Ref. 7, we separate this temperature dependence as follows.
Assume that a variation in the electronic density of states O(9) as a function of volume V and temperature T is given by Then we can write that its variation with temperature at a constant pressure is In Eq. 2, we identify that the full temperature dependence is composed by the variation of D with temperature at a fixed volume plus its variation with volume at a fixed temperature, multiplied by the We address the last term of Eq. 2 by taking the experimental thermal expansion coefficients for graphite and ML-MoS2 from references [8,9] and scaling the ground-state lattice constants by this factor at the temperatures of 50, 150 and 300 K. Since there is a mismatch in the in-plane expansion coefficient of graphite and ML-MoS2, we approximated these lattice constants by their predicted average at each 16 temperature, in order to minimize strain. We then calculated the ground-state O(9) at each volume, consistent with T = 50, 150, 300 K, which is shown in Figure S11. The first term of Eq. 2 describing the electron-phonon effect on D, is more challenging to address, and we detail the procedure in the following section.
Supplementary Figure S11 | The PDOS of F6TCNNQ in the system considering only lattice expansion at different temperatures. Projected density of states has been consistently shifted by 5.4 eV, see main text.

Section 8. Stochastic sampling of the vibrational space
In the framework of the Williams-Lax theory [10,11] and harmonic approximation, the temperature dependence of the electronic density of states D is obtained as a multidimensional Gaussian integral [12,13] :  [14] : v a,W = w a x2y a,W + 1, (5) where w a = x ℏ/2| } n a is the zero-point vibrational amplitude, and y a,W = [7 ℏ d Ä Å f − 1] r= is the Bose-Einstein distribution. Additionally, | É and n a denote the mass of proton and the frequency of the νth normal mode. Taking Eq. 5 together with the set of frequencies and eigenmodes obtained with phonopy, [15] we can create a list of atomic displacements at a given temperature. In the following, the method of Sobol low-discrepancy sequences [16] was used to sample efficiently the normal coordinates. [12] Finally, 35, 50 and 35 configurations were used to calculate the electronic density of states (EDOS) of the F6TCNNQ/ML-MoS2(4×8)/graphite (5×10) [522 atoms] for temperatures 50 K, 150 K, and 300 K, respectively. This amount of sampling ensured sufficient statistical convergence for the effects discussed in this manuscript. In the same way with the density of sates in the Eq. 3, we also calculated the temperature-dependent band gap for the free-standing ML-MoS2 presented in Figure S12 and Table   S3. The band gap becomes smaller as the temperature increases, with VBM and CBM being shifted up and down in energy, respectively, and the shift we observe agrees with the literature. [13,17] We have also calculated the renormalization of the levels of the F6TCNNQ molecule as shown in Table S4. Position of VBM and CBM of ML-MoS2 as a function of temperature corresponding to Figure S12 (HSE06 functional, ZORA, SOC). In parenthesis, relative difference to static, in meV.

Section 9. Minimal model to understand T dependence of level occupations
We write the following single-particle donor-bridge-acceptor Hamiltonian in the basis of the isolated systems (commonly called the diabatic basis) Ñ = Ö V |O⟩⟨O| + ∑ Ö â ä |ã å ⟩⟨ã å | h å + Ö ç |*⟩⟨*| + ∑ p Vâ ä |O⟩⟨ã å | h å + ∑ p çâ é |*⟩èã ê | h ê + p Vç |O⟩⟨*|, (6) where Ö refer to the energy levels of the isolated systems and δ to the diabatic couplings between different levels. Taking the energy Ö V to set the chemical potential in this system, we are interested in knowing, upon interaction, what will be the final population of the A state, i.e., what is the CT from D to A via B. Upon solution of this model, the eigenvalues Ö ê and eigenvectors |ë⟩ of the interacting system are obtained. We are then interested in calculating the population of the original states after interaction. In a single-particle picture, this can be obtained by bó é Ä ï f ⁄ é = ∑ èòôëöèëôòö é ∑ .
bó é Ä ï f ⁄ é , (7) where X is the state of interest. In this expression, the chemical potential is implicitly assumed to be at zero of energy. One can thus, e.g., solve the model by placing the Fermi level of graphite (D) at 0 eV and expressing the relative position of the other three levels with respect to it, taking their values at the potential energy surface, as shown in Table S5. We assume that p Vâ ä = 0.5, p çâ ä = 0.2 and p Vç ä = 0.01 7R. These are empirical values but reasonable for the vdW bonded systems regarded here and the distances involved between the components. [18] A first-principles calculation of these coupling parameters is desirable but could not be achieved yet for these systems with the electronic-structure codes we employ. We vary T for the along the temperatures considered in this study.
We then let the position of the VBM of ML-MoS2 and the acceptor LUMO vary in the range we have calculated in Tables S3 and S4 (we found that varying the CBM made no difference to the results). The variation of the population of the molecular level is shown in Figure S14. This allows us to exemplarily visualize how the population of acceptor state A (LUMO of F6TCNNQ) depends on the energy Ö ç and the energy of the VBM of ML-MoS2 (Ö â= ), which both vary as function of temperature as shown in Tables S3 and S4. We can conclude the following: (i) The direction of the renormalization of the state energies can lead to an appreciable increase of the molecular state population. This is due to the coupling between the different parts of the system and the electronic level renormalization. (ii) The quantitative amount of variation of the population of the molecular levels will depend on the actual values of the coupling.

21
However, varying the coupling constants within a sensible range (0.1-0.8 eV) does not change the qualitative picture shown in this paper.
Supplementary Figure S14 | Population of the molecular level upon a linear variation of the relative VBM and LUMO energies. The blue dot corresponds to the values at the potential energy surface and the red dot corresponds to population with the renormalized levels at 300 K.