Electronic structures and enhanced optical properties of blue phosphorene/transition metal dichalcogenides van der Waals heterostructures

As a fast emerging topic, van der Waals (vdW) heterostructures have been proposed to modify two-dimensional layered materials with desired properties, thus greatly extending the applications of these materials. In this work, the stacking characteristics, electronic structures, band edge alignments, charge density distributions and optical properties of blue phosphorene/transition metal dichalcogenides (BlueP/TMDs) vdW heterostructures were systematically studied based on vdW corrected density functional theory. Interestingly, the valence band maximum and conduction band minimum are located in different parts of BlueP/MoSe2, BlueP/WS2 and BlueP/WSe2 heterostructures. The MoSe2, WS2 or WSe2 layer can be used as the electron donor and the BlueP layer can be used as the electron acceptor. We further found that the optical properties under visible-light irradiation of BlueP/TMDs vdW heterostructures are significantly improved. In particular, the predicted upper limit energy conversion efficiencies of BlueP/MoS2 and BlueP/MoSe2 heterostructures reach as large as 1.16% and 0.98%, respectively, suggesting their potential applications in efficient thin-film solar cells and optoelectronic devices.


Results
The lattice constants of monolayer BlueP, MoS 2 , MoSe 2 , WS 2 and WSe 2 were fully optimized and the values are 3.268, 3.164, 3.295, 3.165, 3.295 Å, respectively, which are in consistent with the reported data listed in Supplementary Table S1 24,26,27,[30][31][32] . The interlayer lattice mismatches between BlueP and MoS 2 , MoSe 2 , WS 2 and WSe 2 were evaluated to be + 3.18%, − 0.82%, + 3.15% and − 0.82%, respectively, which are all in an acceptable range and accessible in experimental synthesis. For each BlueP/TMDs heterostructure, six most possible stacking configurations with six special rotation angles between the adjacent layers were explored, as shown in Supplementary Fig. S1. The rotation angles of BlueP monolayer with respect to TMDs are 0°, 60°, 120°, 180°, 240° and 300°, respectively. All systems are geometrically optimized for getting stable atomic configuration. The difference of total energy between the six different stacking configurations compared with the most stable one, the interlayer distance between monolayer BlueP and TMDs, as well as the M − X (M = Mo, W and X = S, Se) and P-P bond lengths for the BlueP/TMDs heterostructures are listed in Supplementary Table S2. According to the total energy data for each configuration and Eq. (1), configuration (a) is the most stable structure for all the heterostructures. Interestingly, with the rotation of the BlueP monolayer with respect to TMDs, we can find in Supplementary Table S2 that the total energy and interlayer distance of the heterostructures increase and maximize when the rotation angle reaches 180° (configuration (d)) for all the four heterostructures.
In order to investigate the thermodynamic stability of the BlueP/TMDs heterostructures, the formation energies according to Eq. (2) were calculated to be − 165.6, − 210.5, − 161.2 and − 209.7 meV/unit-cell for BlueP/ MoS 2 , BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures, respectively. The negative formation energies indicate thermodynamic stability and the possibility to obtain the BlueP/TMDs heterostructures experimentally. On the other hand, the binding energies 33 between BlueP and TMDs monolayers in BlueP/MoS 2 , BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures are 16.11, 19.46, 15.69 and 19.38 meV/Å 2 , respectively, which are in consistent with the typical vdW binding energy of around 20 meV/Å 2 obtained by the advanced DFT calculations 34 . Therefore, the above BlueP/TMDs heterostructures belong to vdW heterostructures. Considering the fact that each P atoms occupies ~5 Å 2 space of the vdW interface, the vdW binding between the BlueP and TMDs layers are slightly stronger than that of bilayer black phosphorus from the Quantum Monte Carlo study 4 .
The calculated electronic band structures of BlueP and TMDs monolayers are presented in Supplementary  Fig. S2. It is seen that single-layer BlueP is an indirect gap semiconductor with the conduction band minimum (CBM) located between the Γ (0, 0, 0) and Μ (0, 1 2 , 0) points while valence band maximum (VBM) located between the Κ ( 1 3 , 2 3 , 0) and Γ points. The magnitude of the band gap of BlueP is 2.02 eV using optB86b-vdW, which agrees well with the literature values using GGA and DFT-D2 methods (see Supplementary Table S1) 26,27 . The band structures for TMDs monolayers are shown in Supplementary Fig. S2(c-f). It is seen that all TMDs monolayers exhibit direct gap characteristics because all of the CBM and VBM are located at the K point of the hexagonal first Brillouin zone. Furthermore, both VBM and CBM of TMDs mainly originate from the d orbital electrons of the transition metals. The calculated band gap values of monolayer MoS 2 , MoSe 2 , WS 2 and WSe 2 using optB86b-vdW are 1.77, 1.50, 1.89 and 1.63 eV, respectively, which agree well with previous work (see Supplementary Table S1) 32,[35][36][37][38] .
To understand the band offset nature of BlueP/TMDs vdW heterostructures, the band edge alignments, the total density of states (DOS) and orbital-resolved partial DOS were studied systematically, as illustrated in Fig. 1 and Supplementary Fig. S3. For getting self-consistent results, we evaluated the accurate band gap values by the hybrid-DFT method with 8% Hartree-Fock exchange energy based on the optB86b-vdW approach for all the monolayers and heterostructures, which exactly reproduces the experimental gap of monolayer MoS 2 . The calculated band gap using hybrid functional calculations for BlueP, MoS 2 , MoSe 2 , WS 2 and WSe 2 monolayers are 2.39, 1.90, 1.65, 2.08 and 1.75 eV, respectively. These results are in good agreement with previous experimental data 11,30,38 and theoretically reported values, as summarized in Supplementary Table S1 24,32,35,37 . As seen in Fig. 1 heterostructures. In contrast, the CBM primarily originates from the P-3p electrons for BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures; while for BlueP/MoS 2 heterostructure, the CBM is mainly contributed by 4d electrons of Mo atoms. It is consistent to the band edge alignments analysis. As a result, the Fermi level (E F ) of BlueP/TMDs vdW heterostructures shifts and appears between the CBM of BlueP and VBM of TMDs. The work function of a material is a critical parameter commonly used as an intrinsic reference for band alignment 37 . It is seen from Fig. 1 that the predicted work functions of BlueP/MoS 2 heterostructures are smaller than pristine BlueP and MoS 2 monolayers; whereas, the work functions of BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures are smaller significantly than pristine BlueP but larger slightly than TMDs. Owing to the energy of the Se-4p orbital higher than S-3p orbital, the work functions of BlueP/MSe 2 (M = Mo, W) heterostructures are smaller than BlueP/MS 2 heterostructures.
From the calculated band structures in Fig. 2, it is seen that all the BlueP/TMDs vdW heterostructures are identified as indirect band gap semiconductors owing to the indirect band gap nature of BlueP and the interlayer coupling effect in the heterostructures. The BlueP/MS 2 heterostructures can be devided into two groups. Group I includes the BlueP/MoS 2 and BlueP/WS 2 heterostructures (showing in Fig. 2(a,c)), where VBM is located at the Γ point and CBM is located at the K point. This is due to the fact that the BlueP monolayer shares similar valence band level at the Γ point and conduction band level at the K point with single layer MoS 2 or WS 2 . The band interactions coming from the formations of the heterostructures pull up the energy level of the valence band at the Γ point and push down the energy level of the conduction band at the K point, which can also explain the variations of the band alignment as shown in Fig. 1(a,c). Group II involves the BlueP/MoSe 2 and BlueP/WSe 2 heterostructures (showing in Fig. 2(b,d)), where CBM is located between the Γ and Μ points. Although VBM is located at the Γ point for the BlueP/MoSe 2 heterostructure, and which is located at the K point for the BlueP/WSe 2 heterostructure, it is worth noting that the energy level of VBM at the Γ and K points are very close in these cases. For these heterostructures, the BlueP monolayer present deeper conduction and energy band levels than the MoSe 2 or WSe 2 layers. As a result, the interactions of the valence and conduction bands are weak between the BlueP monolayer and MoSe 2 or WSe 2 single layer. Hence the VBM keeps the main features of the single layer MoSe 2 or WSe 2 , and CBM represents the features of the BlueP monolayer between the Γ and Μ points, which agrees well with the band alignment results in Fig. 1(b,d).
Considering the fact that there is a significant spin-orbit coupling (SOC) effect in TMDs 39 , we have calculated the band structures for all the five 2D monolayer materials and four heterostructures including the effect of SOC (illustrated in Supplementary Figs S4 and S5). It is seen that there is no significant SOC effect on band structure of BlueP. Otherwise, SOC splits off the double-degenerated transition metal (Mo or W) d electron occupied valence band state at the K point for the TMDs monolayers and BlueP/TMDs vdW heterostructures. As a results of the band splitting, the VBM shift to the K point for the BlueP/MoSe 2 heterostructure. Apart from this, SOC slightly influences the predicted band structures without changing the characteristics and shapes of the band eigenvalues, which does not change the main results and conclusions of this work.
The calculated band-decomposed charge density is shown in Fig. 3, which provides us a vividly picture to understand the role of the constituent layers in BlueP/TMDs vdW heterostructures. The decomposed charge density was calculated for the lowest unoccupied molecular orbital (LUMO) and the highest occupied molecular orbital (HOMO). The HOMOs are mainly consists of the Mo-4d or W-5d electrons, while the LUMOs are mainly contributed by the P-3p electrons in BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures, except for BlueP/MoS 2 heterostructure (its LUMO is mainly contributed by the Mo-4d electrons), which are in consistent  with the analysis of the electronic structures. The results reveal that CBM and VBM are localized in different monolayers of BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures, resulting in different spatial distribution of the lowest energy electron-hole pairs. Therefore, the MoSe 2 , WS 2 or WSe 2 layer can be potentially used as the electron donor and the BlueP layer as the electron acceptor in the corresponding heterostructures. The free electrons and holes can be effectively separated in the BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures, indicating the potential applications of these vdW heterostructures for optoelectronics and solar energy conversion. Based on the distribution of charge densities, we found that the electrons are confined in the BlueP layer, while the holes in the TMDs side. Therefore, it is expected that spontaneous charge separation occurs when excitons diffuse to the BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures, which is beneficial for the optoelectronic and photovoltaic applications 32 . We have calculated the frontier states to further study the interlayer interactions and hybridizations of the BlueP/MS 2 heterostructures (see Fig. S6), which are beneficial for the microscopic understanding of the conduction channels 33 . We found that the HOMO/LOMO orbitals and frontier states show very similar characters. Comparing to the strong interlayer hybridizations in phosphorene/ SiS heterostructures 40 , the interactions between the BlueP and MS 2 layers are much weaker, which is expected in the vdW heterostructures. Very recently, Zhang et al. found that the indirect band gap nature in BlueP/MoS 2 heterostructure can be tuned to direct by applying external electric fields 35 . We have confirmed the appearance of the indirect to direct band gap transition in the BlueP/MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures under − 0.47, + 0.53 and − 0.47 V/Å external electric fields, respectively (see Fig. 4). Hence, the potential applications of these BlueP/TMDs vdW heterostructures are not restricted by their indirect band gap features.
The calculated plane-averaged electron density differences of BlueP/TMDs vdW heterostructures are shown in Supplementary Fig. S7 to better understand the nature of the bonding mechanism and the charge transfer between the BlueP and TMDs layers, where the magenta regions represent the charge accumulation, and the cyan regions represent depletion in the combined system relative to the two isolated components. As a result of the interlayer coupling effect, there is an obvious space charge accumulation region at the BlueP/TMDs heterointerfaces. In addition, the charge redistribution at the TMDs layer indicates that the formation of vdW heterostructures may affect the M − X bonds (M = Mo, W and X = S, Se).
To understand the effect of quantum confinement on the optical properties of BlueP/TMDs vdW heterostructures, the real (ε 1 ) and imaginary (ε 2 ) parts of the dielectric function were examined and compared to the dielectric function of corresponding single-layer BlueP and TMDs. Figure 5 presents the real parts ε 1 and imaginary parts ε 2 of dielectric function of the hybrid and single systems. For all the BlueP/TMDs vdW heterostructures, both the real and imaginary parts of the dielectric functions show larger values than the individual single systems (see Fig. 5). There are two distinguished peaks in UV light regions (200∼390 nm) and visible light regions (390∼ 770 nm) from the imaginary parts of dielectric functions, indicating enhanced optical properties of all the BlueP/TMDs vdW heterostructures. Figure 6 shows the optical refractive index n(λ), extinction coefficient κ(λ), reflectivity coefficient R(λ) and absorption coefficient α(λ) of BlueP/TMDs vdW heterostructures from the dielectric functions according to Eqs (6)(7)(8)(9). Since all the vdW heterostructures show very similar tendency of the optical properties, we present the BlueP/MoSe 2 heterostructure as an example in detail (see Fig. 6(b)). The refractive index n(λ) of BlueP/MoSe 2 below 300 nm is very close to that of BlueP, but is much higher than that of MoSe 2 above 500 nm. It is seen that extinction coefficient κ(λ) of BlueP/MoSe 2 heterostructure is generally enhanced, especially in UV light ranges.  Moreover, the reflectivity coefficient R(λ) of BlueP/MoSe 2 heterostructure is increased twice compared with the individual monolayers. More importantly, BlueP/MoSe 2 heterostructure shows substantial absorption in both visible light and UV light areas, while BlueP monolayer shows only UV light absorption and MoSe 2 monolayer shows absorption mainly in the visible light region. In addition, its absorption region around 300 and 500 nm is clearly enhanced. Generally speaking, it clearly shows that BlueP/TMDs vdW heterostructures are able to harvest the visible light and UV light more efficiently than the individual single systems. The similar phenomenon can be found in the hybrid GeH/graphene nanocomposite 41 . The heterostructure formation can not only improve the charge separation of photo-induced electron-hole pairs but also broaden the light absorption range 42 . Therefore, the BlueP/TMDs vdW heterostructures possess potential applications in efficient solar cells and ultrathin optoelectronic devices.

Discussion
To go a step further, we introduced the estimation method by Shi et al. 43,44 according to Eqs (10)(11)(12)(13) to calculate the theoretical upper limit conversion efficiency of sunlight for BlueP, TMDs monolayers and BlueP/TMDs vdW heterostructures. The corresponding results are listed in Table 1. It is reported that WSe 2 monolayer has a strong optical absorbance and can be used to fabricate solar cells with light-to-electricity conversion efficiencies of ~0.5% 12 . As seen in Table 1 that our estimated efficiency for WSe 2 with 0.32% agrees well with the experimental result, which verifies the reliability and authenticity of this work. Theoretical calculations also predict that the BlueP/TMDs vdW heterostructures have consistently larger conversion efficiencies than the corresponding individual single systems due to the larger overlap between their absorbance and the solar spectrum. In particular, the predicted upper limit energy conversion efficiencies of BlueP/MoS 2 and BlueP/MoSe 2 heterostructures reach as large as 1.16% and 0.98%, respectively, which are comparable to the graphene/MoS 2 ultrathin heterostructure solar cell with ~1% energy conversion efficiency 45 . The semiconductor solar cell, as a green and efficient technology, has gained wide attention because of its promising applications in solar energy utilization and pollutant elimination. The proposed BlueP/TMDs vdW heterostructures not only exhibit higher energy conversion performance than the individual 2D TMDs monolayers, but also present many distinguished properties, such as good stability, moderate band gap, indirect-direct band gap transition, electronic-hole separation as well as fascinating light adsorption. Therefore, these vdW heterostructures break the limitation of the individual 2D materials, which provide more choices for the scientist and engineer to design the next generation solar cell devices.

Conclusion
In summary, we have systematically studied the stacking configurations induced electronic structures, band edge alignments, charge density distributions and optical properties of BlueP/TMDs vdW heterostructures based on vdW corrected density functional theory. The negative formation energies as well as the good lattice match protect the thermodynamic stability of the BlueP/TMDs vdW heterostructures. All the BlueP/TMDs vdW heterostructures exhibit indirect gap characteristics. The VBM and CBM are localized in different parts of the BlueP/ MoSe 2 , BlueP/WS 2 and BlueP/WSe 2 heterostructures. The MoSe 2 , WS 2 or WSe 2 layer can be used as the electron donor and the BlueP layer can be used as the electron acceptor in the corresponding heterostructures. The BlueP/ TMDs vdW heterostructures exhibit unusually stronger optical absorbance and energy conversion efficiency in the visible range. The results in this work provide a fundamental understanding and guideline for designing high performance optoelectronic and photovoltaic devices.

Methods
Density functional theory calculations. Our calculations were performed based on the density functional theory (DFT) in conjunction with the projector-augmented-wave (PAW) potential as implemented in the Vienna ab initio Simulation Package (VASP) 46,47 . For the exchange-correlation energy, the Perdue-Burke-Ernzerhof (PBE) version of the generalized gradient approximation (GGA) was used 48 . The van der Waals density functional (vdW-DF) of optB86b were considered for all the calculations 49,50 . The valence electron configurations for P, S, Se, Mo and W were 3s 2 3p 3 , 3s 2 3p 4 , 4s 2 4p 4 , 4d 5 5s 1 and 5d 4 6s 2 , respectively. A plane wave cutoff energy of 600 eV was used for the plane-wave expansion of the wave function. The first Brillouin zone was sampled with a fine grid of 7 × 7 × 1 for structure optimization and 15 × 15 × 1 for static calculation 51 . A vacuum of 30 Å along the z direction (perpendicular to the BlueP/TMDs vdW heterostructures layers) was constructed to eliminate the interaction with spurious replica images. The atomic positions were fully relaxed until satisfying an energy convergence of 10 −5 eV and force convergence of 0.01 eV/Å. Due to the fact that the PBE functional underestimates the band gaps of semiconductors, the accurate band gap values and optical properties were further calculated by the hybrid functional approach. As seen in Table S1 and Fig. S8, the standard optB86b-vdW functional underestimates the band gap of well-known monolayer MoS 2 by ~0.13 eV. Comparing to the standard HSE06 hybrid functional with 25% Hartree-Fock exchange energy 52 , we found that the 8% short range Hartree-Fock exchange energy based on optB86b-vdW can reproduce the experimental band gap of monolayer MoS 2 and other TMDs without much deviation of the lattice parameter. Considering the fact that there is no experimental results available for the band gap of monolayer BlueP, we choose monolayer MoS 2 as the reference to predict the band gap of the BlueP/TMDs vdW heterostructures. It is worth noting that the proportion of the short range Hartree-Fock exchange energy influences the absolute values of the predicted band gap only, without changing the variation tendency of the formation of the vdW heterostructures. To further check our results, we have calculated the quasiparticle band gaps of the monolayers and heterostructures using the G 0 W 0 approximation 53 , which are listed in Tables S1 and S3 as well. As seen, our results agrees well with previously scGW 0 results 54 . It is interesting to find in Fig. S8, the G 0 W 0 approximation shows very similar band gap variation pattern to the hybrid functional calculations.
To determine the configuration of the heterostructure stacking, the energy difference Δ E was written as the following equation: where E 0 is the total energy of the most stable configuration, and E i is the total energy of each configuration.
To exam the stability of the heterostructures, the formation energy was obtained according to the following equation: 24 is the sum of the total energy of mutually independent single-layered BlueP and TMDs fixed in the corresponding heterostructure lattice. The work function W was obtained by the following expression: where E vac is the energy of a stationary electron in the vacuum nearby the surface, E F is the Fermi energy level. In addition, the plane-averaged electron density difference along the perpendicular direction to the interface was calculated according to the following equation: Optical properties. Optical properties are described by the photon wavelength dependent dielectric function ε(λ) = ε 1 (λ) + iε 2 (λ), which is mainly contributed from the electronic structures 55 . The refractive index n(λ ) and the extinction coefficient κ(λ ) are determined from ε 1 and ε 2 using 43 λ ε λ ε λ ε λ = + + n( ) 1 2 ( ) ( ) ( ) where λ is the photon wavelength. The reflectivity and absorption coefficient were calculated by 43 and α λ πε λ = ( ) 2 (9) 2 The estimated upper limit to the converted power P is based on the overlap between the solar spectrum and the absorbance: 44 0 0 max where λ max is the longest wavelength that can be absorbed by ultrathin materials and is determined by the lowest-exciton energy E g (i.e., the minimum band gap of the material), where α is the absorption coefficient, and d is the thickness of the simulation cell perpendicular to the layers. The term C(λ) is the conversion factor to account for the fraction of the photon energy converted to lowest-exciton energy (i.e., the thermalization loss). g