First-Principles Study of a MoS2-PbS van der Waals Heterostructure Inspired by Naturally Occurring Merelaniite

Abstract

Vertically stacked, layered van der Waals (vdW) heterostructures offer the possibility to design materials, within a range of chemistries and structures, to possess tailored properties. Inspired by the naturally occurring mineral merelaniite, this paper studies a vdW heterostructure composed of a MoS₂ monolayer and a PbS bilayer, using density functional theory. A commensurate 2D heterostructure film and the corresponding 3D periodic bulk structure are compared. The results find such a heterostructure to be stable and possess p-type semiconducting characteristics. Due to the heterostructure’s weak interlayer bonding, its carrier mobility is essentially governed by the constituent layers; the hole mobility is governed by the PbS bilayer, whereas the electron mobility is governed by the MoS₂ monolayer. Furthermore, we estimate the hole mobility to be relatively high (~10⁶ cm²V⁻¹s⁻¹), which can be useful for ultra-fast devices at the nanoscale.

1. Introduction

Two-dimensional (2D) materials are celebrated among the scientific community, due to the unparalleled exquisite properties compared to their bulk counterparts, arising from the quantum confinement effects. Some of the examples of 2D materials include graphene, h-BN, transition metal dichalcogenides (TMDs), phosphorene, MXenes, etc., which have been explored for novel device applications [1,2,3,4,5,6,7]. The possibility of synthesis of vertically stacked 2D structures, in which the consecutive layers interact through weak van der Waals forces called van der Waals (vdW) heterostructures, can give rise to nanomaterials with desired properties [8,9,10]. The vdW heterostructures present unique electronic properties due to the incommensurate and modulating nature of their structural configurations along the crystallographic a- and b-directions. Despite their fascinating and attractive properties, only a few groups have so far succeeded in the synthesis of such heterostructures [11,12,13,14,15,16,17].

Naturally occurring vdW structures such as cylindrite and franckeite belong to a family of complex, misfit, layered sulfide minerals. These materials are built from the vertically stacked alternating layers of pseudo-quadratic (referred to as Q) and pseudo-hexagonal (referred to as H) layers [18,19,20,21]. In general, the Q and H layers are a (100) slab of the rock-salt-type structure and a CdI₂-type layer, respectively, and the interlayer interaction is primarily dominated by the vdW interactions.

Theoretical and experimental investigations of the 2D heterostructures derived from these incommensurate layered sulfides are still at the stage of infancy. A recent theoretical study of MoS₂-CdS heterostructure reported its enhanced photocatalytic activity [22]. Heterostructures of alternating SnS₂-based (H) and PbS-based (Q) layers have been synthesized by exfoliation of bulk natural franckeite [15,16]. Spectroscopic measurements found the ultrathin franckeite nanosheet to be a p-type material with an estimated bandgap of ≈0.7 eV. Note that the PbS-based Q layers in the franckeite structure are with non-isometric unit-cells and thus structurally different from the ones directly derived from the PbS cubic bulk structure [23]. The Q layer is composed of four atomic layers of sulfide compounds with the formula MX, where M = Pb²⁺, Sn²⁺ or Sb³⁺ and X = S. The H layer consists of octahedrons of disulfide compounds with the formula MX₂, where M = Sn⁴⁺ or Fe²⁺ and X = S [16].

Following the success of fabricating 2D heterostructures from bulk franckeite, we here focus on 2D and bulk MoS₂-PbS heterostructures that are inspired by, but simplified from another naturally occurring member of the cylindrite group, the relatively newly discovered mineral merelaniite (Mo₄Pb₄VSbS₁₅). Merelaniite is composed predominantly of a stacking of a MoS₂-based monolayer (H) incommensurately aligned with a PbS-based bilayer (Q) [24]. Both the H and Q layers of merelaniite are triclinic with space group C1. The unit cell parameters for the H layer are: a = 5.547(9) Å; b = 3.156(4) Å; c = 11.91(1) Å; α = 89.52(9); β = 92.13(5); γ = 90.18(4). For the Q layer, the parameters are: a = 5.929(8) Å; b = 5.961(5) Å; c = 12.03(1) Å; α = 91.33(9); β = 90.88(5); γ = 91.79(4). Note that the a and b directions are parallel to the layers, and c is the stacking direction of the H and Q layers. In this study, the structure and chemistry are simplified to be coherently aligned H and Q layers, composed solely of MoS₂ and PbS, respectively, with a minimally sized commensurate periodic supercell in the a-b plane, as displayed in Figure 1. Specifically, we investigated the stability and electronic properties of this 2D vdW heterostructure employing density functional theory. The formation of such layered PbS precipitates along the (001) plane of molybdenite (MoS₂) mineral was recently reported, thereby suggesting the existence of MoS₂-PbS heterostructure [25]. We also calculate its carrier mobility to assess its applicability in next-generation devices based on 2D materials. We note that 2D vdW heterostructures based on the natural mineral franckeite have been proposed as candidate materials for optoelectronic devices [17,26].

2. Computational Model

Electronic structure calculations based on density functional theory (DFT) were performed employing the Perdew-Burke-Ernzerhof (PBE) functional for exchange and correlation [27]. To account for the vdW long-range interaction between the constituent layers, a posteriori corrected semi-empirical contribution 1/r6 Grimme dispersion correction (D2) was used [28]. The valence electrons were treated explicitly, whereas the core-electrons were treated using the Projector Augmented Wave (PAW) pseudopotentials [29], as implemented in the Quantum Espresso Plane Wave package version 6.4.1 [30] The complete electronic configuration of Mo is [Kr]4d⁵5s¹, of Pb is [Xe]6s²4f¹⁴5d¹⁰6p² and S is [Ne]3s²3p⁴. In the present work, the valence configurations of Mo, Pb, and S are taken to be 4s²5s¹4p⁶5p⁰4d⁵, 6s²6p²5d¹⁰, and of 3s²3p⁴, respectively.

The H layer is a rectangular cell of MoS₂ monolayer with 12 atoms per unit cell and the Q layer is the rectangular PbS bilayer with 8 atoms per unit cell. The 2D MoS₂-PbS heterostructure consists of 1 × 1 cells each of H and Q layers. For the constituent atoms of the heterostructure, the valence electrons were taken to be 14, 14, and 6 for Pb, Mo, and S atoms, respectively.

In our calculations, the wave function and density cut-offs were set to be 40 and 400 Ry, respectively. Atomic positions in the periodic unit cell were allowed to relax until the total forces on each atom were less than ~10⁻⁴ eV/Å, and the self-consistency tolerance was set to 10⁻⁶ Ry. A vacuum spacing of about ~15 Å was employed along the crystallographic c-axis. For the density of states (DOS) calculations, the integration of the Brillouin zone was made on a Monkhorst-Pack grid of (2 × 2 × 1) k-points in the reciprocal space [31]. The room-temperature electron and hole mobilities were calculated by applying a phonon-limited scattering model including the anisotropic characteristics of effective mass following the expression, ES1 (electronic supplementary information (ESI)) given by Bardeen and Shockley [32]. It is worth mentioning that the expression, ES1 (ESI) gives only an estimated value of the carrier mobility in a given material [33].

3. Results and Discussion

3.1. Structural Properties

Considering the complexities of merelaniite’s incommensurate structure and the chemical formula Mo₄Pb₄VSbS₁₅ [22], we focused on a highly simplified model of a 2D heterostructure that consists of one MoS₂ (H) and two PbS (Q) layers (Figure S1a). It is worth noting that a similar assumption was made in the computational investigation of the properties of the 2D franckeite-inspired heterostructures composed of SnS2-based (i.e., H) and PbS-based (i.e., Q) layers [15]. In our model, the 2D MoS₂-PbS heterostructure is simulated by a rectangular unit cell for which the lattice parameters, a and b, are calculated to be 5.64 and 6.40 Å, respectively at the PBE (DFT) + D2 level of theory. Note that the calculated lattice parameters, a, b, and c of the bulk MoS₂-PbS are 5.67, 6.22, and 12.73 Å, showing a good agreement with the corresponding experimental values of merelaniite [24]. In the 2D MoS₂-PbS heterostructure, the H-layer remains nearly strain-free while the Q layer becomes strained whereby its lattice constant contracts ≈8% along with a and expands ≈6% along b in the unit cell. Furthermore, the Q layer is not atomically flat and its non-planar configuration can be described by the buckling height (δ_z) and buckling angle (α) along the z-direction (i.e., perpendicular to the plane) as shown in Figure 1. The structural parameters of the 2D MoS₂-PbS heterostructure are compared with those of the corresponding bulk MoS₂-PbS material in Table 1.

The calculated results predict the stability of the 2D MoS₂-PbS heterostructure with a binding energy of −0.75 eV compared to a value of −1.57 eV for the bulk MoS₂-PbS, which is defined as (E_HQ−E_H–E_Q). Note that the heterostructure is not bound at the DFT level of the theory with binding energy of 0.03 eV. Thus, the inclusion of the D2 term appears to be essential to predict the stability of such 2D heterostructures.

The interplanar distance between the H and Q layer (R_S(H)–Pb(Q)) is 3.07 Å, which is about 7% smaller than that calculated in the bulk material (Table 1). For the Q layer, the buckling-height and buckling-angle values are slightly greater than the corresponding bulk values. This is expected due to the lack of periodicity along the z-axis for a 2D heterostructure. It is also worth mentioning that R_S(H)–Pb(Q) appears to be a representative interplanar distance of the constituent layers interacting via the vdW forces in the 2D MoS₂-PbS heterostructure. We note here that previous calculations on graphene/SnO heterostructures using either PBE-D2 or optB88-vdW functional forms find that the stability of such heterostructures does not sensitively depend on the way we describe the vdW interactions [34].

To evaluate the mechanical stability, we calculate the elastic constants for the 2D MoS₂-PbS heterostructure. The elastic constants for the rectangular symmetry must satisfy the following necessary and sufficient elastic stability conditions, also called Born stability conditions [35]:

C₁₁ > 0 & C₃₃ > 0 & C₁₁C₂₂ > C²₁₂(1)

or

(2)$$\begin{gathered} {\lambda }_I=\frac{1}{2}({\mathrm{C}}_{11}+{\mathrm{C}}_{22}+\sqrt{4{\mathrm{C}}_{12}^2-{({\mathrm{C}}_{11}-{\mathrm{C}}_{22})}^2} ) \\ {\lambda }_{II}=\frac{1}{2}({\mathrm{C}}_{11}+{\mathrm{C}}_{22}-\sqrt{4{\mathrm{C}}_{12}^2-{({\mathrm{C}}_{11}-{\mathrm{C}}_{22})}^2 }) \\ {\lambda }_{III}={\mathrm{C}}_{33}>0 \end{gathered}$$

For the MoS₂-PbS heterostructure, the calculated elastic constants are C₁₁ = 134.3. C₂₂ = 94.5, C₃₃ = 36.1, and C₁₂ = 26.9 GPa, thus satisfying the above conditions for mechanical stability.

3.2. Electronic Properties

Figure 2a shows the electronic band structures of the 2D MoS₂-PbS heterostructure, in which the top of the valence band is associated with the Q layer and the bottom of the conduction band is associated with the H layer. The calculated minimum energy gap is ≈0.2 eV and is predicted to be indirect, with the top of the valance band being at a point between Y and Γ and the bottom of the conduction band at a point between Γ and X. Interestingly, the bandgap of the free-standing H-layer (i.e., MoS₂ monolayer) is calculated to be 1.8 eV and that of the free-standing Q layer (i.e., PbS bilayer) is calculated to be 1.4 eV (Figure S2a,b, (ESI)). This is what we expected because when two layers with different work functions come into contact with each other, the charge redistribution and band bending occur to equalize the work functions between the layers in a heterostructure. Moreover, the weak interlayer coupling also leads to changes in band positions of the constituent monolayers of the heterostructures [36] (Figure S1b).

The interaction of constituent layers via weak vdW forces is affirmed by the partial density of states, which shows that the PbS valence band, constituting the Pb-d and S-p orbitals, move upward towards the Fermi level while the MoS₂ conduction band, comprised of Mo-d and S-p orbitals, moves downward towards the Fermi level for the 2D heterostructure (Figure S2c, (ESI) and Figure 2b). On the other hand, S-p and Mo-d orbitals appear to form a band at about 2 eV below the Fermi energy (Figure 2b). Figure S2c shows the upshifting of S-p and Pb-d orbitals (dominant valence band orbitals) and the downshifting of the Mo-d orbitals (dominant conduction band orbitals) in the heterostructure relative to the individual monolayers.

We should stress here that our calculations did not take account of the spin-orbit coupling (SOC) effects for Pb atoms. However, it has previously been demonstrated that the inclusion of SOC effects does not significantly change the band structure of the 2D franckeite-inspired heterostructure [15].

Since the interlayer interaction between the H and Q layers is dominated by the vdW interactions, we expect a negligible charge transfer between the two constituents. This is affirmed by the calculated Lowdin’s charges, which suggest that ≈0.05 electron transfers from the H layer to the Q layer in the 2D heterostructure. Interestingly, such a small charge transfer still induces an electrostatic potential at the interface (Figure S3 (ESI)) that acts as a barrier between the H and Q layers, thereby promoting separation of electron-hole pairs in the heterostructure. The 2D MoS₂-PbS heterostructure can therefore be considered for optoelectronic applications owing to its narrow bandgap of about 0.2 eV calculated at the PBE (DFT)-D2 level of theory. It is worth mentioning that calculations using a relatively accurate exchange and correlational functional form, the hybrid nonlocal exchange-correlation functional yield the bandgap to be about 0.45 eV (Figure S4, ESI) [28,37]. Such narrow bandgap materials are usually good for thermoelectric applications [38,39].

Carrier transport properties were investigated using the expression ES1 (ESI), which depends upon the in-plane stiffness, carrier effective mass, and deformation potential values listed in Table 2. These values were calculated by applying a strain ranging from −2% to +2% in steps of 0.5% to the 2D MoS₂-PbS heterostructure. Note that the deformation potential of the electrons (holes) can be obtained by linear fitting of the curve of the conduction band edge (valance band edge) versus applied strain for electrons (holes). On the other hand, the in-plane stiffness can be obtained by the quadratic fitting of the expression, ES2 (ESI). We find that the in-plane stiffness for the 2D MoS₂-PbS heterostructure exhibits anisotropic characteristics that are likely due to the Q. layer (Table 1). We estimate the effective mass of electrons or holes by the quadratic fitting of bands near the conduction band minimum (CBM) or valence band maximum (VBM), respectively. A smaller value of a carrier effective mass leads to higher carrier mobility. We find that the carrier effective masses to be anisotropic with holes being heavier along the x-direction and electrons being heavier along the y-direction.

The calculated results show the 2D MoS₂-PbS heterostructure to be p-type with the hole mobility being significantly higher than the electron mobility. Moreover, the mobility anisotropy, which is calculated as R_a = max(μ_x, μ_y)/min(μ_x, μ_y), is predicted to be 3.5 and 634.4 for electrons and holes, respectively in the heterostructure. It is interesting to note that the carrier mobility of the heterostructure is essentially governed by the constituent layers; the hole mobility is governed by the Q layer whereas the electron mobility is governed by the H layer, as shown in Figure 3. This is what we expected, since the constituent layers in the 2D MoS₂-PbS heterostructure are bonded by the weak vdW forces, which do not significantly alter the dispersion of the valence or conduction bands of the constituents (Figure S5, ESI). It is to be noted that the calculated electron mobility of the MoS₂ monolayer (Table 2) is comparable to previously reported experimental values, which are within the 0.5–220 cm²V⁻¹s⁻¹ range depending upon the different measuring techniques. Single-layer MoS₂ fabricated using mechanical peeling or chemical exfoliation techniques and electrically contacted using electron-beam lithography on Si/SiO₂ substrate (which also acts as the gate dielectric) can show field-effect mobilities of up to tens of cm²V⁻¹s⁻¹. The mobility can be improved by depositing a high-κ dielectric material such as HfO₂, which overpowers the Coulomb scattering due to the high- κ environment and modifies the phonon dispersion giving high carrier mobility of ~200 cm²V⁻¹s⁻¹ [40,41,42,43].

4. Conclusions

In summary, the DFT calculations have shown that a merelaniite-inspired MoS₂-PbS heterostructure can be formed from the corresponding bulk MoS₂-PbS material. The heterostructure is predicted to be p-type with a bandgap of about 0.45 eV obtained at the HSE06-DFT level of theory. We found that the carrier mobility in the heterostructure is governed by the constituent layers as the heterostructure is weakly bonded by vdW forces. Moreover, the hole mobility is predicted to be relatively high (~10⁶ cm²V⁻¹s⁻¹), which can be useful for nanoscale devices. Considering the recent fabrication and characterization of 2D franckeite-inspired heterostructures, we expect our theoretical work to motivate experimentalists in fabricating the 2D merelaniite heterostructures for energy-related applications.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/ma14071649/s1, Figure S1. (a) The unit cell of a 2D heterostructure of MoS₂-PbS, and (b) the charge density difference plot for the MoS₂-PbS heterostructure. (blue for depletion and red for accumulation). Isovalue 0.01 e/ų., Figure S2. Calculated band structure of the (free standing) (a) H (i.e. MoS₂ monolayer) and (b) Q (i.e. PbS bilayer) layers. (c) The partial density of states for the H layer (Red), Q layer (Green), and the 2D MoS₂/PbS heterostructure (Blue); (c)-(i) p-orbitals of the S atoms, (c)-(ii) d-orbitals of the Mo atoms and (c)-(iii) d-orbitals of the Pb atoms., Figure S3. The calculated profile of the planar averaged self-consistent electrostatic potential of the 2D merelaniite heterostructure as a function of position., Figure S4. (a) Calculated PBE (DFT)+D2 band structure and (b) HSE06 band structure obtained at the PBE (DFT)+D2 structural configuration. Figure S5. Band energy (E_edge) of CBM and VBM as a function of lattice dilation e along x and y directions for 2D MoS₂-PbS heterostructure.

Author Contributions

Conceptualization, J.A.J. and R.P.; methodology, S.K. and G.D.D.; software, G.D.D. and S.K.; validation, G.D.D., S.K., R.P. and J.A.J.; formal analysis, G.D.D., R.P. and S.K.; investigation, G.D.D., S.K., J.A.J. and R.P.; resources, J.A.J. and R.P.; data curation, G.D.D. and S.K.; writing—original draft preparation, S.K. and G.D.D.; writing—review and editing, G.D.D., S.K., R.P. and J.A.J.; visualization, S.K. and G.D.D.; supervision, R.P. and J.A.J.; project administration, R.P. and J.A.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article or Supplementary Material. The data presented in this study are available in the Supplementary section and can also be obtained from the corresponding author.

Acknowledgments

Helpful discussions with Luca Bindi are acknowledged. The RAMA and Superior High-Performance Computing Clusters at Michigan Technological University, Houghton, MI, USA, were used to obtain the results presented in this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Side view (c-axis vertical), (b) top view of the MoS2-PbS heterostructure, composed of H (i.e., MoS2 monolayer), and Q (i.e., PbS bilayer) layers. (c) The buckling parameters of the Q (PbS-like bilayer) layer. Color code: Mo (grey), Pb (orange), and S (green).

Figure 2. Calculated (a) band structure and (b) partial density of states of the 2D QH MoS2-PbS heterostructure. The conduction band minimum (CBM) and valence band maximum (VBM) are indicated in (a).

Figure 3. Calculated electron (e) and hole (h) mobilities of the 2D MoS2-PbS heterostructure and the free-standing PbS and MoS2 monolayers.

Table 1

Calculated bond distances (in Å) and buckling parameters of the 2D MoS₂-PbS heterostructure consisting of a MoS₂ monolayer (H) and a PbS bilayer (Q). The corresponding values of the bulk MoS₂-PbS material are also listed.

Computed Parameter	Bond Distances (Å)
	2D MoS2-PbS Heterostructure	Bulk MoS2-PbS
RS(H)–Pb(Q)-interplanar	3.07	3.30
RS(H)–Mo(H)-intraplanar	2.41	2.41
RS(Q)–Pb(Q)-intraplanar	2.89	2.88
RS(Q)–Pb(Q)-interplanar	2.73	3.16
Q-layer buckling height (δz)	0.71	0.51
Q-layer buckling angle (α)	13.5◦	9.8◦

Table 2

Calculated in-plain stiffness (C_2D), effective mass (m*), deformation potential (E₁), and mobility (µ) of 2D MoS₂-PbS heterostructure and (free-standing) constituent layers;) and xh (yh), respectively, are electrons and holes along x-direction (y-direction).

		C2D(N/m)	m*(me)	E1(eV)	µ(103 cm2V−1s−1)
H-Qheterostructure	xe	276.3	0.36	5.83	0.90
	xh	276.3	0.91	3.48	1.86
	ye	190.1	0.83	5.50	0.30
	yh	190.1	0.09	0.36	1180.0
(free standing) H layer(MoS2 monolayer)	xe	167.5	0.42	7.8	0.28
	xh	167.5	0.42	2.5	1.10
	ye	164.1	0.83	6.9	0.15
	yh	164.1	3.9	2.0	0.18
(free standing) Q layer(PbS bilayer)	xe	112.0	0.16	3.00	5.6
	xh	112.0	2.4	4.00	0.14
	ye	26.6	0.56	1.45	1.6
	yh	26.6	0.09	0.66	31.0