The cell structures of monolayer P, GaN and GeC are presented in Fig. 1(a-c). The optimized in-plane lattice constants of these monolayers are a = b = 3.27 Å, 3.24 Å and 3.25 Å respectively, which are closely aligned with earlier reports26,35,41,60. Obviously, the lattice misfit between monolayer P and GaN or GeC is less than 1%, which is an ideal condition for the epitaxial formation of outstanding vdWHs. Besides the lattice mismatch, the local atomic stacking arrangement will also affect the interfacial characteristics. Here, we configured P/GaN and P/GeC vdWHs with six different stacking arrangements, as presented in the supplementary information (SI) Fig. S1. To find the most stable stacking arrangements, we evaluated the binding energy (*E**b*) which is defined as *E**b* = *E**vdWH* - *E**P* - *E**GaN/GeC*. Where *E**vdWH*, *E**P*, and *E**GaN/GeC* represent the total energies of the vdWHs, P, and GaN/GeC monolayers, respectively. As listed in table S1, the stacking-A of P/GaN and stacking-D of P/GeC vdWHs having the binding energy of -0.67 eV and − 0.64 eV are the most stable configurations, respectively. Therefore, in the subsequent study, we will focus on the most stable vdWHs, as depicted in Fig. 1(d, e) for P/GaN and P/GeC.

To verify the dynamic stabilities of the P/GaN and P/GeC vdWHs, we first calculated their phonon band spectra. As shown in Fig. 2(a, b), the vdWHs hold a good dynamic stability because they do not have any negative frequencies in their phonon band spectra. Both vdWHs unit cells are composed of four atoms, resulting in 12 distinct vibrational modes, including three acoustic modes and nine optical modes. In Fig. 2(a, b), we highlight the three acoustic modes, i.e. ZA (red), TA (blue), and LA (purple), which have a notable impact on thermal conductivity. The lowest-frequency acoustic mode (ZA mode) displays a parabolic shape in the phonon spectrum. Because the frequency *ω* of ZA mode is a quadratic function of the wave vector *k*, i.e. \(\omega ={\left( {\frac{\kappa }{\rho }} \right)^{1/2}}{k^2}\). This type of acoustic branch is a prevalent characteristic of 2D systems61. It has been reported phonon frequencies are correlated with their atomic masses62. In the P/GaN and P/GeC vdWHs, the lightest N/C and heaviest Ga/Ge atoms primarily act at high and low frequencies. Whereas the intermediate atomic mass P atoms are involved throughout the phonon frequency range. In addition, to evaluate the thermal stabilities, we also performed AIMD simulations of P/GaN and P/GeC vdWHs at 700 K for 18 ps. As displayed in Fig. 2(c, d), both vdWHs are thermally stable at 700K. Because the energy fluctuations are minimal during the whole 18 ps-AIMD, which implies there are no notable structural deformation. Therefore, exploring the thermoelectric characteristics at such a high temperature is also justifiable.

The electronic structure can reflect the fundamental features of electrical transport. The band structure and density of states (DOS) of parent monolayers are shown in Fig. S2. The P and GaN monolayers possess indirect semiconductive band gaps of 1.94/2.75 eV and 1.99/3.26 eV using PBE/HSE06 functional, respectively. Whereas the GeC monolayer has a 2.06/2.81 eV direct band gap using PBE/HSE06 functional. These electronic properties are well consistent with earlier literatures26,35,41,60. Figure 3(a, b) displays the band structures of P/GaN and P/GeC vdWHs attained by the PBE (red) and HSE (blue) functionals, respectively. The estimated energy gaps are presented in Table S2. Both vdWHs exhibit an indirect band gap, as evidenced by the presence of conduction band minima (CBM) located at the M-point and valence band minima situated at the K-point of the irreducible Brillouin zone. Further, we also investigated the DOS of these vdWHs, as illustrated in Fig. 3(c, d). The DOS shows that both vdWHs retain the Type-II band alignment. The VBM of vdWHs arises due to the main involvement of P atom orbitals, whereas the CBM appears due to the primary contribution of C/N atom orbitals. In both vdWHs, band degeneracy is observable in the lower conduction band (CB) region, suggesting that n-type doping exhibits enhanced transport capabilities compared to p-type doping15. Additionally, the observed differences in band gap between heterostructure and their respective monolayers indicate that vdWHs have the potential to modulate the transport capabilities.

The assessment of TE capability primarily relies on two variables: carrier mobility (*µ*) and carrier relaxation time (*τ*). Based on the deformation potential (DP)63 theory, we calculated the relaxation time and carrier mobility for vdWHs10,64. The formulation for the carrier mobility and relaxation time can be written as follows:

$${\mu _{2D}}=\frac{{2e{\hbar ^3}{C_{2D}}}}{{3{K_B}T|{m^*}{|^2}E_{{DP}}^{2}}}$$

4

$$\tau =\frac{{{m^*}\mu }}{e}$$

5

In the preceding expression, *e* stands for the elementary charge, *ħ* is the reduced Planck constant, *m** is the effective mass, *K**B* is the Boltzmann constant and *T* designates the temperature. \({C_{2D}}={1 \mathord{\left/ {\vphantom {1 {{S_0}}}} \right. \kern-0pt} {{S_0}}}({\partial ^2}{E_{total}}/\partial {\varepsilon ^2})\) is the in-plane stiffness of the 2D system. Where *S**0* is the unstrained system surface area. \(E=\partial {E_{total}}/\partial \varepsilon\) is the DP constant, representing a shift of the band edge site upon uniaxial strain. This constant describes the influence of lattice potential field modulation on carriers around the Fermi level, and quantifies the intensity of electro-acoustic coupling. Table S3 illustrates the coefficients linked with the DP theory and the carrier mobility at 300 K for the P/GaN and P/GeC vdWHs. For both vdWHs, the electron mobility is marginally higher than that of hole. The carrier relaxation time at different temperatures is presented in Fig. 4(a). According to the Eqs. (4) and (5), the relaxation time monotonically decreases as temperature increases.

To judge the TE performance of these vdWHs, we will focus on their electrical conductivity (*σ*), Seebeck coefficient (*S*), and electronic thermal conductivity (*κ**e*). The rigid band approximation is employed to examine the impact of doping on electrical transport parameters. Negative and positive values of chemical potential (*µ*) are used to describe the p- and n-types doping, respectively. The *S* of P/GaN and P/GeC vdWHs as a function of *µ* at numerous temperatures are presented in Fig. 5(a, d). The p-type and n-type doping at both sides near the fermi level result in more pronounced and powerful peaks of the *S* at 300 K. The best Seebeck coefficients for P/GaN and P/GeC vdWHs at 300 K are 1793 µV/K and 1210 µV/K with n-type doping, respectively. As reported previously, the Seebeck coefficient for a talented thermoelectric device often exceeds 200 µV/K65,66. We can observe a distinct response to temperature wherein the parameter *S* decreases when subjected to heat. This is a common characteristic of many good TE materials due to their energy-sensitive DOS, which exhibits significant changes (*∂n(ε)/∂ε*) near the Fermi level. As shown in Fig. 5(a, d), the *S* of both vdWHs display peaks in lower *µ* region. The Mott equation can be better to define *S* of thermoelectric materials67:

\(S=\frac{{{\pi ^2}k_{B}^{2}T}}{{3q}}{\left\{ {\frac{1}{n}\frac{{\partial n(E)}}{{\partial E}}+\frac{1}{\mu }\frac{{\partial \mu (E)}}{{\partial E}}} \right\}_{E={E_f}}}\) (6) In the given equation, *q* represents the charge of the carrier, *E**f* denotes the Fermi energy, *µ(E)* stands for mobility, and *n(E)* is associated with the DOS. *S* can be enhanced through two approaches: optimizing localized DOS and adjusting carrier concentration. A more energy-sensitive DOS distribution can lead to an increase in *∂n(ε)/∂ε*. Figure 3 (c, d) illustrates comparable high rates of DOS for P/GaN and P/GeC vdWHs near the Fermi level, which aligns with the order of their *S* coefficients. The monolayer P and GaN/GeC interface may create an energy barrier that filters or scatters low-energy carriers, while high-energy electrons are allowed to pass24. The energy barrier and its carriers scattering effect play a crucial role in the enhancement of *S*. Thus, vdWHs could increase *S* and provide a large *S**2**σ*, which further increases *ZT* as a result10. Figure 5(b, e) presents the *σ/τ* as a function of *µ* for P/GaN and P/GeC vdWHs. In contrast to *S*, the *σ/τ* exhibits a slight variation for both vdWHs and is less affected by the change in temperature. It is evident that the n-type doping of both vdWHs exhibits a higher slope of (*σ/τ*) vs. *µ*, especially for *µ* in the range of -2.0 to 2.0 eV. This observation aligns with the outcomes mentioned above. The conduction electrons adhere to Fermi-Dirac statistics. So, the slopes of *σ/τ* gradually become flattened in the low *µ* region, even when the temperature rises. In order to accurately predict the TE performance of the materials, it is essential to consider electron thermal conductivity (*κ**e**/𝜏*). The *κ**e**/𝜏* of both vdWHs shown in Fig. 5(c, f) are acquired using the Wiedemann-Franz law57. The *κ**e**/𝜏* curves are similar to those of *σ/𝜏*, but are highly responsive to changes in temperature. The *κ**e* gradually increases as the *µ* moves away from the Fermi level (*µ* = 0). The transport coefficients presented above demonstrate the good electronic transport characteristics of the P/GaN and P/GeC vdWHs with n-type doping.

Next, we investigated the lattice thermal conductivity (*κ**l*) of these vdWHs. The lattice thermal conductance is mainly achieved by phonons. A good thermoelectric material should have an ineffective phonon thermal transport in order to preserve a uniform temperature gradient. The group velocity (*ν*), Grüneisen parameter (*γ*), and relaxation time of phonons are essential variables for calculating *κ**l*. The phonon velocity is a key physical parameter for assessing thermal transport. We can determine the phonon group velocities by:

$${\nu _i}{\text{(q)}}=\frac{{\partial {\omega _i}{\text{(q)}}}}{{\partial {\text{q}}}}$$

7

The *∂ω**i*(q) represents phonon frequency. Figure 6(a, d) shows the group velocities of P/GaN and P/GeC vdWHs. Regarding both vdWHs, distinct colors are used to represent the acoustic phonon group velocities for the ZA, TA, and LA modes. Due to the relatively lower acoustic group velocities, especially for the ZA and TA modes, *κ**l* will exhibit a declining tendency. According to Fig. 6(a, d), the velocities in the low-frequency region (about 2 ~ 7 THz) are significantly greater, while they are quite minimal in the high-frequency zone. Meanwhile, there is a slight difference in the group velocity between both vdWHs. Figure 6(b, e) displays the Grüneisen parameters (*γ*) and phonon relaxation time. The Grüneisen parameter's describes the anharmonic interactions among lattice atoms. The Grüneisen parameter can be written as follows:

$${\gamma _i}{\text{(q)}}= - \frac{V}{{{\omega _i}{\text{(q)}}}}\frac{{\partial {\omega _i}{\text{(q)}}}}{{\partial V}}$$

8

The *V* factor in the above equation represents the absolute volume of the crystal. In general, a large |*γ*| implies a significant phonon-phonon anharmonic scattering68, which reduces the *κ**l*, especially at high temperatures. Here, both vdWHs hold a moderate value of *γ*, suggesting a notable anharmonicity. When two phonons collide, anharmonicity (a marker of an imbalance in phonon transport) occurs69. In the low-frequency region, the large |*γ*| strongly suppresses heat transmission and results in low *κ**l*. At low temperatures, anharmonicity is more significant because the nonbonding electrons in the form of lone-pairs usually couple with the valence electrons of the adjacent atoms70. Then the anharmonic third-order IFCs are used to compute the phonon relaxation time, which is a key factor in *κ**l*. Figure 6(c, e) shows that both vdWHs have significantly shorter phonon relaxation times, indicating sufficient phonon scattering, which also infers a low *κ**l*.

Based on the harmonic and anharmonic force constants, we can obtain the temperature dependent lattice thermal conductivity (*κ**l*) of P/GaN and P/GeC vdWHs by

$${\kappa _{l,i}}=\sum\limits_{v} {\sum\limits_{{\text{q}}} {{C_{ph}}v_{{g,i}}^{2}({\text{q,}}v)\tau ({\text{q}},v)} }$$

9

The variable *C**ph* represents the specific heat per mode, \(\nu _{{g,i}}^{2}\) denotes the phonon group velocity and 𝜏(q,*ν*) represents the relaxation time of the phonon mode with wave vector q and dispersion branch *ν*. As illustrated in Fig. 4(b), at 300 K, the *κ**l* values for the P/GaN and P/GeC vdWHs are 6.18 and 8.47 WK−1m−1, which are significantly lower than those of GaSe/SnS2 (14.61 WK−1m−1)62, WSe2/SnS2 (22.69 WK−1m−1)71 and MoS2/WS2 (68.44 WK−1m−1)72 vdWHs. The lattice thermal conductivity declines as the temperature increases, following the expression \({\kappa _l} \propto {1 \mathord{\left/ {\vphantom {1 T}} \right. \kern-0pt} T}\), due to the more intense phonon thermal movement that reduces phonon mean free path and remarkable phonon-phonon scattering upon high temperature73. Acoustical modes play a significant role in controlling the behavior of monolayer and heterostructure systems74. These low-frequency acoustic modes may contribute low phonon group velocities, resulting in low *κ**l*. Ultimately, creating a heterostructure has a notable impact on the transmission of phonon heat.

Finally, based on the electronic and lattice transport parameters above, we calculated the the *ZT* values of P/GaN and P/GeC vdWHs at different temperatures (300 K, 500 K, and 700 K). Figure 7(a, b) demonstrates a better TE performance of n-type doping compared with that of p-type doping. This anisotropic feature can be attributed to *S* and *σ* as shown in Fig. 5. For both vdWHs, they exhibit a larger *ZT* upon higher temperature. Meanwhile, the peak positions of *ZT* also shift towards µ = 0 as the temperature rises. More importantly, on the effective temperature scale, the *ZT* values for P/GaN and P/GeC vdWHs reached remarkable values of 5.07 and 4.67, respectively. The *ZT* of our designed vdWHs is superior to those of previously proposed vdWHs and single layers, as displayed in Fig. 8(a)10,26–28,45,62,75–77. The thermoelectric conversion efficiency (*η**TEG*) of a TEG can be estimated using the following formula78:

$${\eta _{TEG}}=\frac{{{T_h} - {T_c}}}{{{T_h}}}\frac{{\sqrt {1+Z{T_{avg}}} - 1}}{{\sqrt {1+Z{T_{avg}}} +{{{T_c}} \mathord{\left/ {\vphantom {{{T_c}} {{T_h}}}} \right. \kern-0pt} {{T_h}}}}}$$

10

Where *ZT**avg* represents the average *ZT* over the temperature range. *T**h* and *Tc* are the temperatures of the material on hot and cold surfaces, respectively. TE materials with *ZT**avg* > 1 have the potential to be practically used in a wide temperature range 79. Remarkably, the *ZT**avg* for P/GaN and P/GeC vdWHs are 3.91 and 3.44, respectively, which are much larger compared with the well-known SiPGaS/As, SnSe/SnS, SnSe and SnS TE systems10,28. In addition, based on the *ZT**avg*, we calculated the thermoelectric conversion efficiency (*η*TEG) by Eq. (10). The P/GaN and P/GeC vdWHs have much better conversion efficiency, up to 26% and 25% at 700 K, respectively. Compared with P/GeC, the P/GaN vdWH has a larger band gap, resulting in a larger *S* and a lower *κ**l*. This makes P/GaN a potential superior thermoelectric material. Our analysis offers valuable insights into thermoelectric transport, which can be utilized to design devices using vdWHs, thus driving advancements in thermoelectric technology.