Effective control of the interfacial Schottky barrier height (SBH) in nanoelectronic devices is crucial for enhancing device performance^1–6. However, conventional experimental measurements and theoretical calculations often require substantial resource investment and time, which significantly limits the rapid screening and optimization of emerging two-dimensional material systems. Therefore, developing efficient theoretical and computational methods that can rapidly and accurately predict and modulate SBH has become one of the key challenges in current nanoelectronics research^7–14.

Among the wide array of 2D semiconductors, transition metal dichalcogenides (MX₂; M = transition metal; X = S, Se, Te; e.g., MoS₂, MoSe₂, WTe₂) have attracted considerable attention due to their layer-dependent tunable bandgaps and outstanding optoelectronic properties. These features make MX₂ promising candidates for high-performance van der Waals heterostructures. In particular, heterostructures composed of graphene (C) and MX₂ exhibit favorable lattice compatibility and reasonable interfacial binding energies, resulting in structurally stable systems with desirable electronic characteristics^11,15–17. Nonetheless, the existence of Fermi level pinning effects and the C-MX₂’s non-negligible SBH, which hinder efficient charge injection and limit the potential of these heterostructures in device applications^18,19. Precise tuning of the interfacial SBH is therefore essential for realizing optimal device performance. Various modulation strategies have been explored, including interlayer spacing engineering, strain modulation, chemical doping, and the application of external electric fields^6,8–11. Within these, electric field modulation stands out for its non-destructive nature, fast switching capability, and ease of integration into practical devices. Notably, recent studies have demonstrated that vertical electric fields can significantly reduce the SBH in C-MoSe₂ and C-WSe₂ heterostructures, even enabling a transition from Schottky to Ohmic contact^6,11. However, a comprehensive understanding of the intrinsic mechanisms by which electric fields modulate SBH in C-MX₂ heterostructures remains elusive. Furthermore, there is still a lack of robust, quantitative models capable of accurately predicting SBH variations under different field conditions. Addressing these challenges is essential for the rational design of next-generation 2D electronic and optoelectronic devices.

The emergence of machine learning (ML) techniques, particularly symbolic regression approaches such as the SISSO, has provided a novel theoretical framework and methodological foundation for addressing the aforementioned challenges²⁰. Compared with conventional ML methods, SISSO not only enables efficient identification of key descriptors from high-dimensional feature spaces, but also yields physically interpretable predictive models. This has led to remarkable progress in various domains, including the prediction of catalytic activity, electrocatalytic performance, and Schottky barrier heights in two-dimensional materials^7,8,21,22. Therefore, integrating density functional theory (DFT) with physically interpretable ML approaches like SISSO for the prediction and modulation of SBH in 2D heterostructures represents a highly innovative and promising research direction.

Motivated by the above background, this study aims to systematically investigate the precise modulation mechanisms of both n-type and p-type Schottky barrier heights (SBH) in five C-MX₂ heterostructures (CrSe₂, MoS₂, MoSe₂, WS₂, and WSe₂) under external electric fields. In addition, we integrate the SISSO machine learning method to construct an efficient and interpretable predictive model for SBH. The overall research workflow is illustrated in Fig. 1, and includes the following steps: (1) Screening suitable candidate materials from the Computational 2D Materials Database (C2DB) and constructing corresponding supercells; (2) Building vertical van der Waals heterostructure models of graphene and C-MX₂; (3) Performing first-principles (DFT) calculations to obtain key parameters such as SBH, work function, and interfacial charge transfer under varying external electric fields; (4) Constructing a dataset based on the calculated results and simple physical quantities (SPQ) of the each C-MX₂’s constituent elements and elemental forms; (5) Applying SISSO and other machine learning techniques to identify key descriptors and develop physically interpretable, efficient, and reliable SBH prediction models. This study not only provides theoretical guidance for interfacial engineering of 2D material heterostructures, but also lays a solid theoretical and methodological foundation for the design and development of high-performance nanoelectronic devices in future integrated circuits.

To further analyze the effect of the external electric field on the SBHs of C-MX₂ heterostructures, an out-of-plane electric field ranging from − 0.5 to 0.5 V/Å (with a step size of 0.1 V/Å) was applied along the Z-axis, and the band structures under different field strengths were calculated. The corresponding results are shown in Fig. S2. The calculations reveal that the C-MX₂ heterostructures exhibit significant responses to the applied electric field. For instance, C-CrSe₂, C-MoS₂, and C-WS₂ transition from Schottky contact to quasi-Ohmic contact within the field ranges of − 0.3 to 0.0 V/Å, − 0.4 to 0.0 V/Å, and − 0.3 to 0.1 V/Å, respectively, while beyond these ranges, they exhibit Ohmic contact behavior. Similarly, C-MoSe₂ and C-WSe₂ undergo a transition to quasi-Ohmic contact within the ranges of − 0.4 to 0.4 V/Å and − 0.3 to 0.4 V/Å, respectively. The specific threshold values for contact type transitions are summarized in Table S2. These results indicate that the SBHs of C–MX₂ heterostructures are highly sensitive to external electric fields, and contact type transitions can be achieved within a relatively narrow field range. However, the number of data points currently available is insufficient for robust machine learning–based SBH prediction. Therefore, to enlarge the dataset and enable more detailed evaluation of field-induced modulation, the electric field intervals were further refined using the quasi-Ohmic transition field as the interval boundary. As a result, a total of 77 samples were generated. The specific electric field values corresponding to each sample in the subdivided intervals are listed in Table S3. The results show that as the electric field increases from negative to positive, the Фₙ decreases while the Ф_p increases. This trend is approximately linear overall, reflecting the high sensitivity of SBH to external electric fields. Such behavior is highly desirable for realizing efficient electrical control and multi-state switching functionalities in nanoelectronic devices. What’s more, the SBH values Ф_n and Ф_p for each C-MX₂ heterostructure under different electric field strengths were extracted and plotted in Figs. 4a-e, showing the variation trends with respect to the applied electric field. Figure 4f illustrates the field thresholds at which contact type transitions occur for each system. In all heterostructures, Ф_n decreases monotonically while Ф_p increases monotonically as the electric field strength increases. Near the quasi-Ohmic region, both trends exhibit some nonlinearity, whereas in the Schottky regime, the relationships are nearly linear. Additionally, E_g of each C-MX₂ remains almost constant under different electric field conditions, which allows Ф_n to be directly inferred from E_g and Ф_p. Given the current challenges in tuning p-type field-effect transistors (FETs), the subsequent machine learning model development will primarily focus on the prediction and modulation of Ф_p.

A deeper understanding of the mechanism by which external electric fields modulate the SBH in C-MX₂ heterostructures can be achieved by introducing the pinning factor S, which quantifies the Fermi-level pinning effect and is defined as follows³²:

where W_C denotes the work function of the graphene layer in the C-MX₂ heterostructure under an applied electric field. The pinning factor S reflects the extent to which the Fermi-level responds to external tuning. When S approaches 1, the Fermi level is highly tunable and there is minimal pinning; when S approaches 0, the Fermi level is strongly pinned and SBH modulation becomes ineffective. Figure 5a presents the relationship between Фₙ and the graphene work function W_C for five C-MX₂ heterostructures. The results exhibit strong linear correlations between Фₙ and W_C, with the slopes of the fitted lines representing the pinning factors S. The specific values of S are labeled in the Fig. Experimentally, Sotthewes et al.³³ used conductive atomic force microscopy (C-AFM) and scanning tunneling microscopy (STM) to analyze Fermi-level pinning behavior at both pristine and defect sites by forming Schottky contacts between various probes (e.g., Pt, PtSi, and n-doped Si) and MX₂ materials (M = Mo, W; X = S, Se, Te). The reported pinning factors for pristine regions were: MoS₂ (0.30 ± 0.03), MoSe₂ (0.19 ± 0.03), WS₂ (0.21 ± 0.03), and WSe₂ (0.28 ± 0.03). In the present work, the constructed C-MX₂ heterostructures do not contain obvious defects such as point or line defects, apart from minor lattice distortions. The computed S values range from 0.17 to 0.30, which are consistent with the experimentally reported values for pristine MX₂ regions. This indicates that C-MX₂ heterostructures do exhibit a certain degree of Fermi-level pinning, while still maintaining good tunability under external electric fields. Moreover, a comparison between Fig. S4 and Fig. 4 shows that the work function W_C of monolayer graphene responds more significantly to electric field modulation than the Φ_n, further confirming that external fields primarily modulate the SBH by tuning the W_C.

Traditional experimental and theoretical approaches for determining the SBH often require significant resource investment and time, which severely limits their applicability for large-scale materials screening. Against this backdrop, machine learning (ML) has emerged as an efficient and promising alternative. Compared with conventional ML models, the SISSO method demonstrates distinct advantages. On the one hand, it efficiently selects key features from a vast pool of candidate descriptors in a high-dimensional feature space; on the other hand, SISSO can derive analytical expressions with clear physical meaning, thus overcoming the interpretability limitations typically associated with traditional "black-box" models. As a result, SISSO-based descriptors not only enable rapid prediction of SBH in heterojunctions but also provide valuable theoretical insights for the rational design and optimization of materials. In building the model, particular attention was paid to ensuring both physical interpretability and computational efficiency of the descriptors. An ideal SISSO expression should satisfy two key criteria: (1) it should rely as much as possible on the intrinsic physical properties of constituent elements; (2) it should minimize dependence on intermediate variables derived from DFT calculations. Following these principles, the constructed feature space includes the following key parameters: applied electric field strength (E), lattice constant of the heterojunction (a, b, c), work function of graphene (W_C), work function of the heterojunction (W_h), interlayer charge transfer (ΔQ_C), carrier tunneling probability (P_TB), and basic physical properties of the constituent elements and their elemental forms (see Table S7 for details). The initial dataset consists of 77 samples with 47 feature dimensions, 83% of which (39 features) can be directly obtained from databases.

To ensure both the accuracy and generalization ability of the model, all samples were first randomly shuffled. The top 70% (i.e., 54 samples) were then selected as the training set for model development. As a result, a concise descriptor with clear physical significance was obtained:

This descriptor profoundly elucidates the physical and chemical mechanisms governing the modulation of the Schottky barrier. Firstly, the electric field strength E, as the dominant factor, directly regulates interfacial charge distribution and band alignment. Its exponential form, e^E, reflects the non-linear enhancement effect of the electric field on the SBH, which is consistent with the monotonic increase of Φ_p with the electric field observed in Figs. 4, S2, and S3. It is noteworthy that the SISSO algorithm automatically discarded redundant electronic structure features (such as P_TB and ΔQ_C) during the optimization process, indicating that the electric field E alone is sufficient to capture the core physics underlying SBH variations, demonstrating the model's exceptional feature compression capability.

From the perspective of interface chemistry, the lattice constant a serves as a structural descriptor that effectively distinguishes the lattice backgrounds of different C-MX₂ systems, reflecting the constraint of heterojunction interface matching degree on electronic coupling. The incorporation of the chalcogen atomic radius (R_m2) from the transition metal dichalcogenide highlights the critical influence of the interfacial chemical environment on Schottky behavior. In the C-MX₂ heterostructure, chalcogen atoms (X) are in direct contact with the graphene layer, and their atomic size modulates the van der Waals gap and orbital overlap between layers, thereby influencing charge transfer and the formation of interface dipole moments. This finding aligns with the reported view that chalcogen elements play a dominant role in interfacial electronic structure³⁴, revealing the intrinsic relationship between chemical composition and interfacial physical properties at the atomic scale. As shown in Fig. 6a, a highly linear relationship exists between the descriptor δ and Φ_p (R² = 0.944), indicating the excellent predictive capability of this descriptor for the p-type Schottky barrier height.

For further evaluation of the generalization and physical transferability of the SISSO model, telluride-based systems (C-MoTe₂ and C-WTe₂) were selected for validation. To ensure the lattice mismatch remains within 5% (4.77% for C-MoTe₂ and 4.80% for C-WTe₂), the corresponding lattice constants differ significantly from those of the previously studied sulfide and selenide systems (see Fig. S7). In terms of electronic structure, telluride systems typically exhibit stronger interlayer charge transfer and narrower band gaps (0.855 eV for MoTe₂ and 0.906 eV for WTe₂)^35,36, posing a greater challenge for descriptor-based prediction. Similar to Fig. 6a, Fig. 6b shows the relationship between descriptor δ and Φ_p for C-MoTe₂ and C-WTe₂, respectively. Both exhibit strong linear trends, with R² values of 0.997 and 0.995. Despite notable differences between Te-based and S/Se-based materials—such as Dirac point alignment, carrier effective mass, and structural geometry—the SISSO model maintains high predictive accuracy. This further confirms the generalizable physical mechanism extracted by the SISSO method, namely the synergistic effects of chalcogen atomic radius, heterostructure lattice constant, and applied electric field in modulating the Schottky barrier height. These findings not only enhance the credibility of the model but also reveal common principles in interfacial engineering of 2D heterostructures, providing theoretical support for developing predictive frameworks for SBH. Looking forward, the inclusion of more diverse 2D material systems (e.g., phosphorene, hexagonal boron nitride) is expected to further extend the applicability of the model and accelerate the rational design of novel nanoelectronic devices.

In this study, SISSO-based machine learning has been carried out to systematically investigate the electric-field modulation of Schottky barrier heights in C-MX₂ van der Waals heterostructures. Our structural analyses confirm that all constructed heterostructures exhibit low lattice mismatch, appropriate interlayer distances, and binding energies characteristic of stable van der Waals systems, providing a reliable foundation for subsequent electronic property exploration. We further demonstrate that external vertical electric fields effectively modulate the SBH, enabling transitions between Schottky and Ohmic contacts. This tunability arises primarily from the field-induced changes in the work function of graphene and the resultant charge redistribution at the interface, as quantified by the Fermi-level pinning factor S and planar-averaged charge density difference analysis. Notably, the carrier tunneling probability, though modest, responds monotonically to the applied field, underscoring the potential for electrical control of carrier injection in nanoelectronic devices. A key achievement of this work is the development of a physically interpretable descriptor δ for predicting p-type SBH, derived via SISSO from a minimal set of features: the applied electric field strength, the heterostructure lattice constant, and the chalcogen atomic radius. This descriptor not only achieves high predictive accuracy (R² = 0.944) on the training set but also generalizes effectively to telluride-based heterostructures (C-MoTe₂ and C-WTe₂), despite their distinct electronic properties. The success of δ highlights the critical roles of interfacial chemical environment and structural matching in governing Schottky behavior. Our findings provide a robust theoretical framework and a efficient predictive tool for designing and optimizing 2D material–based heterojunctions. The methodology presented here—combining high-throughput DFT simulations with interpretable machine learning—can be extended to other 2D systems and interface properties, accelerating the discovery and development of next-generation electronic and optoelectronic devices.