Quasiparticle electronic structure of phthalocyanine:TMD interfaces from first-principles GW

Interfaces formed between a pair of semiconductors feature intriguing electronic and optoelectronic properties due to the alignment of the energy levels of one component with respect to the other at the interface. The direction of the charge transfer in the exciton splitting process depends on type of heterojunction,¹ which, in turn, depends on the electronic structure of the interface. Mixed-dimensional heterostructures formed between 0D molecules and 2D substrates,^2,3 in particular, organic molecules deposited on monolayer transition metal dichalcogenides (TMDs),^4–6 have attracted much attention in recent years due to their promising applications in optoelectronic devices,⁷ photovoltaics,⁸ and photocatalysis.⁹ Among all molecular adsorbates, (metallo)phthalocyanines are notably interesting^10,11 by virtue of their structural planarity, large π-conjugation, photon-absorbing capability,¹² fruitful surface chemistry,¹³ as well as great tunability in electronic and optical properties via a change in the metal center.^14–16 Once the two photon absorbers—a (metallo)phthalocyanine molecule and a monolayer TMD—are placed together to form an interface, the underlying electronic structure, i.e., the relative alignment of energy levels of the two components, dictates the mechanism and direction of the charge transfer across the interface, giving rise to distinct optoelectronic properties. To quantitatively characterize the electronic structure at these interfaces, various experimental efforts have been put forward,^3,17–20 while benchmark results and a systematic account of the trends are still missing, which constitute the main goals of this paper.

Complementary to experimental techniques, first-principles calculations play a unique role in elucidating the electronic structure and structure–property relationship via modeling of atomistically well-defined systems. Notably, the energy levels at a heterogeneous interface pertinent to charge transfer are quasiparticle levels, whose accurate characterization, in principle, requires methods beyond the conventional density functional theory (DFT). Although many-body perturbation theory (MBPT), such as the GW formalism,^21–23 has been very successful^24–26 in resolving the gap problem of DFT,^27,28 its relatively high computational cost hinders its routine applications in large systems such as the interfaces formed between (metallo)phthalocyanines and monolayer TMDs. As a result, most prior computational studies of the phthalocyanine:TMD interfaces still employ DFT with different complexities such as semi-local functionals,^17,29 hybrid functionals,^3,30 or the DFT + U approach.³¹ 

The key ingredient in GW that is responsible for an accurate determination of interfacial energy level alignment is the so-called surface polarization^32,33 or, equivalently, dielectric screening due to the substrate. To effectively capture the dielectric screening while reducing the computational cost, the substrate screening GW was proposed in Ref. 34, which was shown to be accurate for weakly coupled interfaces with negligible orbital hybridization.^34–38 In this work, we apply this approach to a series of interfaces: metal-free phthalocyanine (H₂Pc) adsorbed on monolayer MX₂ (M = Mo, W; X = S, Se) and zinc phthalocyanine (ZnPc) adsorbed on MoX₂ (X = S, Se). We focus on the quasiparticle electronic structure, especially the frontier energy level alignment between the valence band maximum (VBM) and the conduction band minimum (CBM) of the monolayer TMD with the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of the H₂Pc or ZnPc. Our results reveal the structure–property relationship and provide GW-quality benchmark results for future studies.

Furthermore, in most experimental studies, a substrate is used to support the composite molecule:TMD interface system from the bottom, with a commonly used one being α-quartz (SiO₂).^7,20,39–41 Due to the dielectric screening effect, substrates could greatly affect the electronic properties of the adsorbate on top of it.^42–45 In this work, we use the dielectric embedding GW approach⁴⁶ to address the dielectric effect of the SiO₂ substrate on the H₂Pc:MoS₂ interface. We show that with additional screening from the SiO₂ underneath, the energy level alignment at the H₂Pc:TMD interface is modulated considerably. Our results on the embedded H₂Pc:MoS₂ system agree quantitatively with experimental measurements.⁴¹ 

This paper is organized as follows: In Sec. II, we detail the computational methodology and parameters. In Sec. III, we present our results in two aspects: the structure–property relationship across a series of systems and the dielectric effect of the SiO₂ substrate on the H₂Pc:MoS₂ interface. We then conclude in Sec. IV with brief remarks.  Appendixes A and  B are devoted to draw a quantitative connection between two GW-based methods—interface GW and projection GW— to supplement our discussion in Sec. III B.

As a first step, we relax the in-plane lattice parameter and atomic coordinates of each monolayer TMD unit cell using the vdw-DF-cx functional.⁴⁷ This is the functional that we will use to relax the structure of the interfaces, so we also employ it here to ensure consistency. The calculation uses a k-mesh of 18 × 18 × 1 and a kinetic energy cutoff of 100 Ry. All DFT relaxations employ the optimized norm-conserving Vanderbilt (ONCV) pseudopotentials^48,49 and the package.⁵⁰ The resulting lattice constants are 3.15, 3.29, 3.15, and 3.28 Å for monolayers MoS₂, MoSe₂, WS₂, and WSe₂, respectively. These results agree very well with experimental measurements, which yield 3.15,⁵¹ 3.30,⁵² 3.15,⁵³ and 3.28 Å⁵³ for the four systems, respectively.

After the monolayer unit cell relaxation, we build 6 × 6 supercells and place one H₂Pc molecule flat on each substrate and one ZnPc molecule flat on the MoS₂ and MoSe₂ substrates to form six interface systems. Each interface simulation cell is 30.0 Å along the c direction and includes about 23 Å of vacuum. During the relaxation of the interface, the atoms belonging to the substrate are kept fixed in their relaxed monolayer positions to ensure the exactness of the subsequent reciprocal-space folding of the noninteracting polarizability. The coordinates of the atoms belonging to the adsorbate molecule are fully relaxed until all residual forces are below 0.05 eV/Å. The relaxations are carried out using the vdw-DF-cx functional,⁴⁷ a k-mesh of 3 × 3 × 1, and a kinetic energy cutoff of 70 Ry. We found an adsorption height of about 3.0 Å for each system we study, similar to the result of a prior calculation³ (3.3 Å). In our relaxed structures, the center of the H₂Pc or ZnPc molecule is approximately at the bridge position of two S or Se atoms belonging to the top layer of TMD, which is a stable binding site for similar systems.³¹ We note that Ref. 54 reported that the binding energies and band alignments are roughly the same for different binding sites. Figure 1 shows the relaxed H₂Pc:MoS₂ structure in two different views.

Considering the large size of the interfaces and the associated computational cost of conventional GW, we apply the substrate screening GW approach³⁴ for all systems and have explicitly benchmarked this approach against direct GW calculations of the H₂Pc:MoS₂ and H₂Pc:MoSe₂ systems. All GW calculations are performed at the G₀W₀ level and employ a mean-field starting point using the Perdew–Burke–Ernzerhof (PBE) functional,⁵⁶ the Hybertsen–Louie generalized plasmon-pole model²³ for the frequency dependence of the dielectric function, the semiconductor screening for the treatment of the q → 0 limit, the slab Coulomb truncation⁵⁷ for the removal of spurious long-range interactions along the c direction, and the static remainder⁵⁸ in the self-energy calculation to improve convergence, as implemented in the package.⁵⁹ We note that our calculations do not include the spin–orbit coupling, which is known to cause a 0.4–0.5 eV splitting for WS₂ and WSe₂ and a 0.1–0.2 eV splitting for MoS₂ and MoSe₂ in the valence band.^54,60

Here, we list the computational parameters involved in the substrate screening GW calculations. For the calculation of the noninteracting polarizability of the substrate unit cell, $χsub0$⁠, we use a q-mesh of 18 × 18 × 1, a 5 Ry dielectric cutoff, and 200 bands in the summation. For the treatment of the q → 0 limit, 30 bands on a shifted q-grid are used. For the calculation of the noninteracting polarizability of the adsorbate molecule, $χmol0$⁠, we use a simulation cell that is of the same size along a and b as the interface but is much smaller in size (10 Å) along c. We use a q-mesh of 3 × 3 × 1, a 5 Ry dielectric cutoff, and 2400 bands in the summation. For the treatment of the q → 0 limit, 360 bands on a shifted q-grid are used. After that, the $χsub0$ is folded in the reciprocal space to a 6 × 6 supercell, and the $χmol0$ is mapped in the real space to the interface simulation cell, following Ref. 34. These quantities are then combined at each q-point to approximate $χtot0$⁠, the noninteracting polarizability of the interface. It is then inverted to generate the dielectric function in the interface simulation cell, after which the self-energies are computed (for the details of different approaches employed in this work, see below), where we use a k-mesh of 3 × 3 × 1, a 5 Ry dielectric cutoff, and 7200 bands of the interface in the summation for the Green’s function.

We consider two types of self-energy calculations for each interface, in line with our prior works. (i) “Interface GW” calculations, where we compute the expectation value of the self-energy operator using ϕ^tot, an orbital of the interface, i.e., ⟨ϕ^tot|Σ|ϕ^tot⟩. (ii) “Projection GW” calculations as proposed in Refs. 61 and 62, where we compute the expectation value of the self-energy operator using an orbital of the freestanding substrate (ϕ^sub) or that of the freestanding monolayer of adsorbate (ϕ^mol), i.e., ⟨ϕ^sub|Σ|ϕ^sub⟩ for the former and ⟨ϕ^mol|Σ|ϕ^mol⟩ for the latter. In both the approaches mentioned above, the self-energy operator Σ is calculated using iG^totW^tot, i.e., both G and W are from the interface system with W calculated using the substrate screening approximation. We compare the interface GW and projection GW results for every system and show that they agree very well except for the LUMO of H₂Pc when H₂Pc is adsorbed on a MoS₂ substrate. In  Appendixes A and  B, we show that this discrepancy is due to the strong orbital hybridization and explicitly establish a quantitative connection between the two.

At the end of this section, we clarify the terminologies we use in this paper. All results reported in Sec. III B are obtained using the “substrate screening GW”³⁴ approach, where we compute the self-energy operator using iG^totW^tot, together with the approximation $χtot0≈χmol0+χsub0$⁠. Within this framework, two types of self-energies are computed, “interface GW,” ⟨ϕ^tot|Σ|ϕ^tot⟩, and “projection GW,” ⟨ϕ^mol|Σ|ϕ^mol⟩ or ⟨ϕ^sub|Σ|ϕ^sub⟩, as discussed above. Note that the “substrate screening GW”³⁴ is physically identical to the “XAF-GW” approach,³⁶ while the difference between the two is only technical: the former first computes the $χmol0$ in a small cell and then performs a real-space mapping of this quantity to the interface simulation cell and hence is more computationally efficient than the latter.

All results reported in Sec. III C are obtained using the “dielectric embedding GW” approach⁴⁶ (but without the real-space truncation of $χsub0$ as done in Ref. 46). Here, we compute ⟨ϕ^ad|Σ|ϕ^ad⟩, where Σ = iG^adW^tot. One can see that the key difference here is the use of Green’s function of the adsorbate, G^ad, in the self-energy operator compared to the “projection GW” within the substrate screening framework. Because of this difference, the “dielectric embedding GW” further reduces the computational cost without sacrificing the accuracy for weakly coupled interfaces. Note that in Sec. III C, the “adsorbate” itself is already a H₂Pc:MoS₂ interface (i.e., both G^ad and ϕ^ad are for the H₂Pc:MoS₂), and only the SiO₂ substrate is included implicitly as a dielectric environment through W^tot such that the “weak coupling” condition is satisfied.

GW calculations are known to converge slowly. We describe our convergence study in this section to show that our calculations are reasonably well converged and the level alignment results are reliable.

For the monolayer MoS₂ unit cell, we have checked that compared to a 10 Ry dielectric cutoff and 1000 bands in the summation, our choice of parameters (5 Ry and 200 bands) leads to a convergence in the bandgap within 0.02 eV, although the individual quasiparticle energies are off by about 0.3 eV. Our prediction of the monolayer MoS₂ gap is 2.81 eV, in good agreement with prior GW calculations⁶³ (2.80 eV). We also note that for such low-dimensional materials, the nonuniform neck subsampling (NNS) method⁶⁴ leads to faster convergence. Although we do not use it here for the interface, our calculations of the monolayer MoS₂ unit cell achieve reasonably good agreement (within 0.1 eV in the gap) with NNS results with a 10 Ry dielectric cutoff.

For the interface systems, using H₂Pc:MoSe₂ as an example, we have checked that our choice of parameters (a 5 Ry dielectric cutoff and 7200 bands) leads to a convergence of the quasiparticle energies and energy level alignments at the interface within 0.05 eV compared to using a 7.5 Ry cutoff and 9000 bands in the summation. This choice of parameter is thus adopted for all other interface systems.

When a molecule is adsorbed on a semiconductor substrate, the interface could, in principle, exhibit the so-called type-I (straddling bandgap), type-II (staggered bandgap), or type-III (broken bandgap) energy level alignment.¹ Specific to the H₂Pc:TMD or ZnPc:TMD interfaces studied in this work, type-I and type-II heterostructures are possible based on our results below. Depending on the relative ordering between TMD and molecular levels, we further categorize the interfaces to type-Ia, type-Ib, type-IIa, and type-IIb, as schematically shown in Figs. 2(a)–2(d), respectively. All interfaces studied in this work have a direct bandgap at Γ, so all interface gaps and energy level alignments are reported for the Γ point.

In Fig. 2, we show frontier energy levels (bands) of the freestanding monolayer TMD and those of the freestanding molecular layer, together with their counterparts within the interface. Blue lines represent TMD levels, and red lines represent molecular levels. We discuss the following quantities that characterize the electronic structure of these interfaces: $ΔTMD0(Δmol0)$ is the bandgap of the TMD monolayer (H₂Pc or ZnPc molecular layer) in its freestanding form, while Δ_TMD (Δ_mol) is the gap of the TMD (H₂Pc or ZnPc) within the interface system. Δ_LL (Δ_HH) is the gap between the TMD CBM (VBM) and the LUMO (HOMO) of the molecular layer, which is of interest for both type-I and type-II heterostructures because the sign and magnitude of Δ_LL (Δ_HH) dictate the direction and barrier for electron (hole) transfer across the interface, respectively. For the type-II heterostructures in Figs. 2(c) and 2(d), we further consider Δ_HL, the fundamental (transport) gap of the entire interface, which is between the HOMO of the molecule (the VBM of the TMD) and the CBM of the TMD (the LUMO of the molecule) for type-IIa (type-IIb). $ΔTMD0$ is calculated at the K point of the Brillouin zone for the TMD unit cell, and all other quantities are calculated at the Γ point.

Table I shows the computed results for all quantities labeled in Fig. 2 from both DFT and substrate screening GW. To verify that the substrate screening approximation holds for the systems, we compare substrate screening GW results with direct GW calculations for two interfaces, H₂Pc:MoS₂ and H₂Pc:MoSe₂. This comparison shows that the substrate screening GW is very accurate: for all quantities reported in Table I, substrate screening GW leads to an agreement with direct GW results within 0.05 eV. With the same accuracy, the substrate screening GW approach costs only about 30% in computing time and 10% in memory (for the χ step) compared to the direct GW for the systems studied. We also note that our GW results on H₂Pc:MoS₂ and ZnPc:MoS₂ are largely in agreement with other calculations of the same systems (but with slightly different simulation cells) using range-separated hybrid functionals⁶⁵ and GW.⁵⁴ 

Figure 3 shows the GW interfacial energy level alignment for the six heterostructures, where we use different colors to represent different substrates or molecules, and all energy levels are measured with respect to a common vacuum. In Fig. 3, solid bars or lines are interface GW results (same as those reported in Table I), i.e., ⟨ϕ^tot|Σ[G^totW^tot]|ϕ^tot⟩, where the ϕ^tot is chosen as the interface orbital that mostly resembles the orbital of interest (HOMO, LUMO, VBM, or CBM) of the freestanding monolayer TMD or molecular layer, as quantitatively determined from orbital projections. Dashed lines are projection GW results, i.e., ⟨ϕ^sub|Σ[G^totW^tot]|ϕ^sub⟩ for the TMD, where ϕ^sub is the VBM or CBM of the freestanding TMD, and ⟨ϕ^mol|Σ[G^totW^tot]|ϕ^mol⟩ for the molecule, where ϕ^mol is the HOMO or LUMO of the freestanding molecular layer. For all cases, except for the H₂Pc LUMO on MoS₂ substrate, the interface GW results agree very well with projection GW results, which indicates negligible orbital hybridization upon formation of the interface such that ⟨ϕ^tot|ϕ^sub⟩ and ⟨ϕ^tot|ϕ^mol⟩ are close to unity. The special case of H₂Pc LUMO on the MoS₂ substrate is discussed in  Appendixes A and  B.

Below in this section, we discuss the results from three aspects: (i) the renormalization of gaps upon the formation of the interface, (ii) the qualitative difference between DFT and GW in the prediction of the type of some heterostructures, and (iii) the structure–property relationship across the six different systems.

For the freestanding monolayer TMD, our GW calculations yield bandgaps of 2.81, 2.40, 3.05, and 2.65 eV for MoS₂, MoSe₂, WS₂, and WSe₂, respectively. Our results are in good agreement with Ref. 63, where scGW₀ calculations of monolayer TMDs with similar lattice parameters (within 0.01 Å of our relaxed values) were performed, resulting in bandgaps of 2.80, 2.40, 3.11, and 2.68 eV, respectively, for the four materials. For the freestanding H₂Pc (ZnPc) molecular layer, our GW calculations yield a HOMO–LUMO gap of 3.86 (3.91 eV) and an ionization potential of 6.24 (6.16 eV), on par with Ref. 66, where GW calculation of an isolated H₂Pc molecule yielded a bandgap of 3.67 eV and an ionization potential of 6.08 eV. The difference between our results and Ref. 66 lies in the physical difference between a periodic molecular layer and an isolated molecule.

At the interface, the gap renormalization for the H₂Pc or ZnPc molecule is significant with about 1 eV decrease in the HOMO–LUMO gap when the molecule is brought in contact with the substrate. This is consistent with well-established understanding of the surface renormalization at molecule–substrate interfaces.^32,33 On the contrary, Table I shows that TMD bandgaps are still very similar to those of the freestanding phase, indicating negligible gap renormalization for the substrate, which we attribute to the relatively low coverage of the molecule on the substrate (see Fig. 1). This situation is different from a previous semiconductor–semiconductor interface that we studied before,³⁷ where we observed gap renormalization on both sides of an organic bulk heterojunction. For the H₂Pc:MoS₂:SiO₂ interface, the MoS₂ gap is significantly renormalized due to the dielectric screening effect of the extensive SiO₂ substrate underneath, as we will elaborate in Sec. III C. Needless to say, DFT clearly underestimates the relevant gaps and does not capture the gap renormalization at the interface.

A remarkable observation here is that for some systems, DFT and GW yield qualitatively different heterostructure types, as listed in Table I. To be specific, for H₂Pc:MoSe₂ and H₂Pc:WSe₂, DFT predicts a type-II, while GW predicts a type-I interface; for H₂Pc:WS₂, DFT predicts a type-I, while GW predicts a type-II interface; for ZnPc:MoSe₂, although both DFT and GW predict type-I, the relative ordering of TMD and molecular levels are opposite (type-Ia and type-Ib as shown in Fig. 2). In all cases, it is the self-energy correction to the LUMO of the molecular adsorbate that pushes this orbital upward in energy, causing a change in the heterostructure type (there is an additional effect in ZnPc:MoSe₂, where the self-energy correction to the HOMO of ZnPc shifts this orbital downward). We note that the experimental characterizations of these quasiparticle energy level alignments require techniques such as (inverse) photoemission spectroscopy, while most existing experiments^{3,7,18,20,67,68} focused on the description of excitons at such interfaces using, e.g., photoluminescence. Although we have not found direct experimental verification of most of the quasiparticle energy level alignments that we have computed (except for H₂Pc:MoS₂:SiO₂), we believe our work provides a reference point for future experiments.

Table I and Fig. 3 reveal the structure–property relationship that is helpful in materials design of similar systems. All sulfur-based interfaces form type-IIa heterostructures with the HOMO (LUMO) of the adsorbed molecule lying higher than VBM (CBM) of the TMD substrate. All selenium-based ones form type-Ia heterostructures with the HOMO (LUMO) of the adsorbed molecule lying lower (higher) than the VBM (CBM) of the TMD substrate. The gap renormalization of H₂Pc is larger on Mo-based substrates compared to the W-based counterparts. Moreover, comparing ZnPc with H₂Pc on the same substrates, we find that the gap renormalization of ZnPc is smaller, resulting in larger Δ_LL and Δ_HH values than H₂Pc-based interfaces. Since these values represent the charge transfer barrier across the interface in the exciton splitting process, this trend suggests that the electron or hole transfer rates across H₂Pc-based interfaces might be generally higher than those across ZnPc-based interfaces.

In typical experimental studies, the molecule:TMD interfaces are further supported by a substrate underneath the monolayer TMD.^41,42 It is well known that the bandgaps of TMDs are sensitive to the dielectric environment provided by other adjacent 2D materials or substrates.^43,44,69 Therefore, it is imperative to include the dielectric screening effects from any additional substrates to achieve quantitative agreement with experimental characterization of the electronic structure of the molecule:TMD interfaces. In this work, we employ the dielectric embedding GW approach as developed in Ref. 46 to include the effect of the substrate using the H₂Pc:MoS₂ system as an example. We focus on this system because an experimental measurement is available⁴¹ for a direct comparison.

A commonly used substrate in experimental studies of molecule:TMD interfaces is the α-quartz (SiO₂).^{3,7,20,39–41} To model the composite H₂Pc:MoS₂:SiO₂ system, we have applied a 5% compressive strain to the experimental lattice constant of SiO₂ to enforce a commensurate simulation cell with H₂Pc:MoS₂. Under this strain, the PBE bandgap of 5.93 eV for bulk SiO₂⁷⁰ increases by about 4% while other qualitative features of the band structure remain intact. Using the dielectric embedding GW approach,⁴⁶ we first compute the noninteracting polarizability of the SiO₂ substrate in its unit cell, then fold this quantity in reciprocal space to the supercell, and finally, combine the folded quantity with the noninteracting polarizability of the H₂Pc:MoS₂ (which, in turn, is calculated using the substrate screening GW³⁴ approach). As a result, the self-energy calculation is only explicitly performed for the H₂Pc:MoS₂ interface, where the W includes the dielectric effect of the SiO₂ substrate.

The dielectric embedding GW approach⁴⁶ applied here adds little extra computational cost compared to the GW calculation of H₂Pc:MoS₂ without the SiO₂ substrate. The inner products between the H₂Pc:MoS₂ orbitals of interest (the orbitals that represent the resonances of H₂Pc HOMO/LUMO and MoS₂ VBM/CBM) and the H₂Pc:MoS₂:SiO₂ orbitals are close to unity, suggesting negligible orbital hybridization between H₂Pc:MoS₂ and SiO₂ and the validity of the dielectric embedding GW approach. We note in passing that within this approach, the orbital hybridization within H₂Pc:MoS₂, i.e., between H₂Pc and MoS₂, is captured exactly.

Figure 4(a) shows the optimized structure of the H₂Pc:MoS₂:SiO₂ system, where we consider a two-layer SiO₂ substrate with silicon termination. Other terminations exist,^71,72 and we have checked explicitly that for a hydrogen-passivated SiO₂ surface, the band alignment results stay unchanged compared to what we will report below. The separation between the SiO₂ and MoS₂ is optimized using vdw-DF-cx to be 3.5 Å. Figure 4(b) shows the quasiparticle electronic structure of the H₂Pc:MoS₂ interface embedded in the dielectric environment of SiO₂, computed using the dielectric embedding GW approach,⁴⁶ with key energy gaps listed in Table I. Compared to the H₂Pc:MoS₂ system without the SiO₂ substrate, the MoS₂ bandgap is further renormalized to 2.07 eV due to the additional dielectric screening from SiO₂, which is in very good agreement with the experiment (2.10 eV as in Ref. 41). The H₂Pc gap is further renormalized only moderately, possibly due to its large distance to the SiO₂ substrate. Furthermore, the Δ_LL is considerably changed to 0.68 eV (compared to 0.14 eV without the SiO₂ substrate), and the Δ_HL is changed to 1.62 eV (compared to 2.33 eV without the SiO₂ substrate). Interestingly, the Δ_HH does not change compared to the case without the SiO₂ substrate.

As a direct comparison with the experiment, Ref. 41 characterized the energy level alignment of the H₂Pc:MoS₂ interface deposited on the SiO₂ substrate. Ultraviolet photoemission spectroscopy and inverse photoemission spectroscopy measurements were performed for pristine monolayer MoS₂ and the H₂Pc:MoS₂ interface, both on the SiO₂ substrate. By aligning the Fermi level of the two systems, Ref. 41 deduced that Δ_TMD = 2.1 eV, Δ_LL = 1.0 eV, Δ_HL = 1.2 eV, Δ_HH = 0.9 eV, and Δ_mol = 2.2 eV. We note that our results agree quantitatively with Ref. 41 in Δ_TMD and Δ_mol values, but our computed H₂Pc HOMO and LUMO levels are both lower by 0.3–0.4 eV than those reported in Ref. 41, resulting in lower Δ_LL, lower Δ_HH, and higher Δ_HL. We discuss two possible sources for the discrepancy. First, Ref. 41 aligned the Fermi level of the pristine MoS₂ (without the H₂Pc adsorbate) and that of the H₂Pc:MoS₂ interface. The precise position of the Fermi level within the bandgap might depend on the specific experimental condition, while it is not well-defined in first-principles calculations at zero temperature. This difference between an experiment and a computation affects how the MoS₂ and H₂Pc levels are relatively aligned without affecting the bandgap of each component. Second, the coverage of the H₂Pc adsorbate on the MoS₂ substrate might be different between Ref. 41 and our modeling, which might lead to different interface dipoles and consequently different Δ_LL, Δ_HL, and Δ_HH values.

The result is a clear indication of the dielectric screening effect of the substrate on molecule:TMD interfaces. Our dielectric embedding GW approach provides an accurate account of this effect without further increasing the computational cost compared to calculations that only involve molecule:TMD interfaces. Finally, we note in passing that our calculated GW bandgap is 8.69 eV for bulk SiO₂, in agreement with Ref. 73, which reports a GW gap of 8.77 eV. Our GW gap is 7.10 eV for the freestanding bilayer SiO₂ with Si termination and 6.64 eV when the bilayer SiO₂ is in contact with monolayer MoS₂, as shown in Fig. 4(b). These gaps are large enough such that the precise positions of SiO₂ band edges relative to TMD bands will unlikely have large effects on the energy level alignment at the H₂Pc:MoS₂ interface.

In this work, we have systematically and quantitatively characterized the quasiparticle electronic structure of a series of interfaces formed between (metallo)phthalocyanine molecules and monolayer TMDs using first-principles GW calculations. Besides the well-known gap renormalization of the adsorbate, we have found that in certain cases, GW and DFT yield qualitatively different heterostructure types, which is of interest for future experimental validation. Furthermore, we have elucidated the dielectric screening effect of the SiO₂ substrate on the electronic structure of H₂Pc:MoS₂, leading to quantitative agreement with existing experiments. Our findings provide a useful structure–property relationship for materials design, insight into charge transfer processes across such molecule:TMD interfaces, as well as benchmark data for future theoretical and experimental investigations.

Z.-F.L. acknowledges support from NSF CAREER Award No. DMR-2044552. O.A. acknowledges A. Paul and Carole C. Schaap Endowed Distinguished Graduate Award and Rumble Fellowship from Wayne State University. This research used computational resources at the Wayne State Grid and additionally via a user project at the Center for Nanoscale Materials at Argonne National Laboratory, an Office of Science User Facility, which was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. Furthermore, large-scale calculations used computational resources from the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We investigate the quantitative difference between projection GW calculations (dashed lines in Fig. 3) and interface GW calculations (solid lines in Fig. 3). They agree very well except for the H₂Pc LUMO on the MoS₂ substrate, which involves strong orbital hybridization. Here, we quantitatively connect the projection GW results with interface GW results.

To make the discussion self-contained in the appendix, we repeat here that in “projection GW,”^61,62 we compute ⟨ϕ^sub|Σ|ϕ^sub⟩ for the TMD substrate and ⟨ϕ^mol|Σ|ϕ^mol⟩ for the adsorbed molecule, where ϕ^sub (ϕ^mol) is an orbital of the freestanding substrate (molecule). In “interface GW,” we compute ⟨ϕ^tot|Σ|ϕ^tot⟩ with ϕ^tot, where ϕ^tot is an orbital of the interface that mostly resembles a substrate or molecular orbital (i.e., a resonance). In both cases, Σ is the same operator that involves the G and W of the entire interface.

We focus on the H₂Pc LUMO on the MoS₂ substrate and expand the H₂Pc LUMO (calculated for the freestanding H₂Pc molecular layer in the same simulation cell as the interface) in terms of the orbitals of the H₂Pc:MoS₂ interface,

Here, we have neglected expansion coefficients whose magnitude is below 0.05. The projection GW computes $⟨ϕLUMOmol|Σ|ϕLUMOmol⟩$⁠. Substituting Eq. (A1), this quantity is related to diagonal and off-diagonal matrix elements of Σ involving the three interface orbitals on the right-hand side of Eq. (A1). These matrix elements are from interface GW calculations, which are then connected to projection GW results that involve the left-hand side of Eq. (A1). For diagonal elements, we use the difference between the quasiparticle energies and the corresponding Kohn–Sham eigenvalues. For off-diagonal elements, we use the matrix elements of Σ − V_xc, with V_xc being the exchange-correlation potential operator.

The expansion coefficients and matrix elements (in eV) are

Using these values from an interface GW calculation, one can calculate from Eq. (A1) that $⟨ϕLUMOmol|Σ|ϕLUMOmol⟩$ = 1.011 eV. This is in very good agreement with the projection GW calculation of the H₂Pc LUMO on MoS₂, whose self-energy correction is 1.023 eV. The above analysis successfully explains the difference between the solid and dashed lines in Fig. 3: the solid line computes $⟨ϕCBM+2tot|Σ|ϕCBM+2tot⟩$⁠, which only involves the first term on the right-hand side of Eq. (A1), while the dashed line computes $⟨ϕLUMOmol|Σ|ϕLUMOmol⟩$⁠, which involves the left-hand side of Eq. (A1).

As a limiting case, if only one of the coefficients in Eq. (A1) is unity, then the isolated molecular orbital and its resonance at the interface are the same (the weak-coupling limit). The projection GW and the interface GW will then yield identical values in the quasiparticle energy. This is, in fact, a very good approximation for most interfaces. Specific to the systems studied in this work, the two approaches agree well for every resonance except for the H₂Pc LUMO on the MoS₂ substrate.

In a similar manner to the above analysis, we can derive interface GW results from projection GW calculations of the substrate and of the adsorbate, which can also explain the difference between the solid and dashed lines in Fig. 3. For the H₂Pc LUMO on the MoS₂ substrate, its resonance in the interface is the CBM+2 of the interface. We then consider the following expansion:

We again neglect expansion coefficients whose magnitude is below 0.05. The left-hand side of Eq. (B1) is the orbital used in an interface GW calculation, and the right-hand side consists of orbitals used in projection GW calculations of the substrate and of the adsorbate, respectively.

The expansion coefficients and the matrix elements (in eV) are

and the off-diagonal elements are all close to zero.

Using these values from projection GW calculations of the TMD and of the molecule, one can derive from Eq. (B1) that $⟨ϕCBM+2tot|Σ|ϕCBM+2tot⟩=0.654$ eV. This is in very good agreement with the interface GW calculation of the CBM + 2 of the H₂Pc:MoS₂ system (the resonance of H₂Pc LUMO), whose self-energy is 0.667 eV.