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

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 MoS2 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 MoS2 monolayer. Furthermore, we estimate the hole mobility to be relatively high (~106 cm2V−1s−1), which can be useful for ultra-fast devices at the nanoscale.


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 aand 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 pseudohexagonal (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 2 -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 2 -CdS heterostructure reported its enhanced photocatalytic activity [22]. Heterostructures of alternating SnS 2 -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 2+ , Sn 2+ or Sb 3+ and X = S. The H layer consists of octahedrons of disulfide compounds with the formula MX 2 , where M = Sn 4+ or Fe 2+ and X = S [16].
Following the success of fabricating 2D heterostructures from bulk franckeite, we here focus on 2D and bulk MoS 2 -PbS heterostructures that are inspired by, but simplified from another naturally occurring member of the cylindrite group, the relatively newly discovered mineral merelaniite (Mo 4 Pb 4 VSbS 15 ). Merelaniite is composed predominantly of a stacking of a MoS 2 -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 2 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 2 ) mineral was recently reported, thereby suggesting the existence of MoS 2 -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]. ≈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 2+ , Sn 2+ or Sb 3+ and X = S. The H layer consists of octahedrons of disulfide compounds with the formula MX2, where M = Sn 4+ or Fe 2+ and X = S [16].
Following the success of fabricating 2D heterostructures from bulk franckeite, we here focus on 2D and bulk MoS2-PbS heterostructures that are inspired by, but simplified from another naturally occurring member of the cylindrite group, the relatively newly discovered mineral merelaniite (Mo4Pb4VSbS15). Merelaniite is composed predominantly of a stacking of a MoS2-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 MoS2 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 (MoS2) mineral was recently reported, thereby suggesting the existence of MoS2-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].
The H layer is a rectangular cell of MoS 2 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 2 -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 −4 eV/Å, and the self-consistency tolerance was set to 10 −6 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].

Structural Properties
Considering the complexities of merelaniite's incommensurate structure and the chemical formula Mo 4 Pb 4 VSbS 15 [22], we focused on a highly simplified model of a 2D heterostructure that consists of one MoS 2 (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 2 -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 2 -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 2 -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 2 -PbS heterostructure are compared with those of the corresponding bulk MoS 2 -PbS material in Table 1.  The calculated results predict the stability of the 2D MoS 2 -PbS heterostructure with a binding energy of −0.75 eV compared to a value of −1.57 eV for the bulk MoS 2 -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 2 -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 2 -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]: For the MoS 2 -PbS heterostructure, the calculated elastic constants are C 11 = 134.3. C 22 = 94.5, C 33 = 36.1, and C 12 = 26.9 GPa, thus satisfying the above conditions for mechanical stability. Figure 2a shows the electronic band structures of the 2D MoS 2 -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 2 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).

Electronic Properties
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 2 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    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 MoS2 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 valency 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 MoS2-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 MoS2-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, 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 2 -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 2 -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 2 -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 2 -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 2 -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 2 monolayer (Table 2) is comparable to previously reported experimental values, which are within the 0.5-220 cm 2 V −1 s −1 range depending upon the different measuring techniques. Singlelayer MoS 2 fabricated using mechanical peeling or chemical exfoliation techniques and electrically contacted using electron-beam lithography on Si/SiO 2 substrate (which also acts as the gate dielectric) can show field-effect mobilities of up to tens of cm 2 V −1 s −1 . The mobility can be improved by depositing a high-κ dielectric material such as HfO 2 , which overpowers the Coulomb scattering due to the highκ environment and modifies the phonon dispersion giving high carrier mobility of~200 cm 2 V −1 s −1 [40][41][42][43].
The calculated results show the 2D MoS2-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 Ra = 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 MoS2-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 MoS2 monolayer (Table 2) is comparable to previously reported experimental values which are within the 0.5-220 cm 2 V −1 s −1 range depending upon the different measuring techniques. Singlelayer MoS2 fabricated using mechanical peeling or chemical exfoliation techniques and electrically contacted using electron-beam lithography on Si/SiO2 substrate (which also acts as the gate dielectric) can show field-effect mobilities of up to tens of cm 2 V −1 s −1 . The mobility can be improved by depositing a high-κ dielectric material such as HfO2, which overpowers the Coulomb scattering due to the high-κ environment and modifies the phonon dispersion giving high carrier mobility of ~200 cm 2 V −1 s −1 [40][41][42][43].   Table 2. Calculated in-plain stiffness (C 2D ), effective mass (m*), deformation potential (E 1 ), and mobility (µ) of 2D MoS 2 -PbS heterostructure and (free-standing) constituent layers;) and xh (yh), respectively, are electrons and holes along x-direction (y-direction).

Conclusions
In summary, the DFT calculations have shown that a merelaniite-inspired MoS 2 -PbS heterostructure can be formed from the corresponding bulk MoS 2 -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 6 cm 2 V −1 s −1 ), 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.