Electronic and vibrational properties of van der Waals heterostructures of vertically stacked few-layer atomically thin MoS2 and BP

Materials Today Communications 19 (2019) 383–392

Contents lists available at ScienceDirect

Materials Today Communications journal homepage: www.elsevier.com/locate/mtcomm

Electronic and vibrational properties of van der Waals heterostructures of vertically stacked few-layer atomically thin MoS2 and BP C.E. Ekuma

⁎,1

T

, S. Najmaei, M. Dubey

Sensors and Electron Devices Directorate, United States Army Research Laboratory, Adelphi, MD 20783, United States

A R T I C LE I N FO

A B S T R A C T

Keywords: 2D Materials Van der Waals Heterostructure First-principles Calculations Band- and Defect-engineering Nanoelectromechanical applications

The state-of-the-art heterostructure-based devices often involve stacks of epilayers of few nanometer thick crystals. However, the ultimate limit would be a hitherto single-atomic-layer structure. Using material-by-design approach, the flexibility of two-dimensional (2D) materials could be explored to develop a multilayered artificial van der Waals (vdW) heterostructure where the synergistic coupling between the individual materials and the underlying interface creates new and novel functionalities. Herein, based on first-principles calculations, we report that vertically stacked few-layer vdW heterostructure of monolayer molybdenum disulfide (MoS2) and hexagonal boron phosphide (BP) is a strong electrically narrow direct bandgap material. More intriguingly, fewlayer vdW heterostructure of monolayer MoS2 and BP exhibits superior mechanical properties. In striking contrast to heterostructures of MoS2/graphene, MoS2/WS2, and few-layer graphene, where the 2D moduli are lower than the sum of the 2D modulus of each nanosheet, the mechanical strength of few-layer hybrid vdW heterostructure of monolayer MoS2 and BP rather increased with increasing thickness. We attribute this difference to the unique interlayer and interface coupling, which affected the vibrational properties of the heterostructures including the emergence of shear and breathing phonon modes as well as the transformation of flexural phonon modes. Such superior mechanical properties could be explored for nanoelectromechanical device applications, e.g., as a nanoresonator. We demonstrate that the electronic structure could as well be tuned with both increasing numbers of monolayer stacks and defect-engineering promising for low-power applications.

1. Introduction The stacking of thin films of at least two different materials to form a heterostructure is at the center of modern semiconductor devices. Devices built with heterostructures have been shown to generally led to improved properties [1,2]. Unlike conventional materials, the properties of two-dimensional materials (2DMs) are easily tunable through rational design of the material structure [3,4]. Hence, exploring the flexibility in 2DMs could lead to better device performance and aid the further development of heterostructure-based devices at the ultimate limit of a discrete atomic layer. Unique to heterostructures of 2DMs is the rather weak van der Waal (vdW) interactions between the adjacent layers in contrast to the strong covalent or ionic bonds that mediate the interactions between adjacent heterostructures constructed with conventional materials. At the heart of this technological revolution is the monolayers of graphite (graphene) and transition metal dichalcogenides. In recent years, the heterostructures of the above-mentioned ones and many other 2DMs have been fabricated and extensively studied (see e.g., Ref. [2]). For example, the Schottky barrier in MoS2-metal

⁎ 1

interface was tuned by adopting BN (or graphene)-molybdenum disulfide (MoS2) heterostructure [5]. The work of Tsai et al. [6] reported a power conversion efficiency of ∼5.23% in the heterostructure of monolayer MoS2/Si-based solar cells and a photogain > 108 was demonstrated in MoS2-graphene heterostructure [7]. Falin et al. [8] reported superior mechanical properties of few-layer BN nanosheets with an almost thickness independent 2D mechanical strength. Many other improved device performances have been demonstrated in heterostructures of 2DMs [1,9,10]. We have selected the combination of monolayer MoS2 and hexagonal boron phosphide (BP) due to their interesting nanoelectromechanical properties, specifically, MoS2, which have been studied for diverse applications and BP a family of the boronbased compounds that is largely unexplored. Motivated in part by these emerging functionalities in 2DMs hybrid structures, we explore the stacking effects of nanosheets of MoS2 and BP on the mechanical, vibrational, electronic, and the impact of interfacial sulfur vacancy. In particular, we aim to explore the possibility of designing a hybrid vdW heterostructure formed from monolayer stacks of MoS2 and BP to create a material with superior properties that will

Corresponding author. Current address: Department of Physics Lehigh University Bethlehem PA 18015 United States.

https://doi.org/10.1016/j.mtcomm.2019.03.005 Received 4 December 2018; Received in revised form 7 March 2019; Accepted 15 March 2019 Available online 16 March 2019 2352-4928/ © 2019 Elsevier Ltd. All rights reserved.

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

Table 1 The calculated crystal lattice parameters, electronic and mechanical properties, and elastic constants of pristine MoS2, BP, and the various vdW hybrid structuresa. Pristine materials C11 C12

E2D

Name

a

SG

Eg

me(kx/ky)

mh(kx/ky)

MoS2 BP BM B2M B2M2 B5M5

3.183 3.213 6.376 6.376 6.376 6.376

D3h[187] D3h[187] C3v [156] D3h[187] C3v [156] C3v [156]

1.70 0.90 0.89 0.85 0.79 0.76

0.56/0.55 0.21/0.21 0.20/0.18 0.21/0.18 0.20/0.25 0.21/0.18

BM B2M

12.732 12.732

C3v [156] C3v [156]

0.44 0.79

3.39/3.32 0.20/0.25

0.60/0.62 132.22 33.88 123.54 0.19/0.20 147.16 39.90 136.34 0.19/0.17 287.15 78.21 265.85 0.18/0.17 441.96 122.37 408.08 0.18/0.23 573.35 155.46 531.20 0.18/0.16 1492.28 393.31 1388.61 Disordered vdw heterostructuresb 0.16/0.17 278.83 72.64 259.91 0.18/0.23 431.22 115.77 400.14

G = C66

γ2D

D

ϵ′c

ν

49.11 53.47 104.47 159.78 209.17 520.02

83.05 93.53 182.68 282.17 364.42 942.80

6.75 0.54 85.24 331.57 1041.58 22725.18

34.57 2.50 202.81 513.92 1240.23 9632.24

0.26 0.27 0.27 0.28 0.27 0.27

100.91 156.18

175.73 273.50

83.12 321.78

202.29 508.65

0.26 0.27

a The lattice constant a (Å), crystal symmetry SG (space group in brackets), energy bandgap Eg (eV), electron me and hole Mh effective masses in units of free electron mass m0, components (C11, C12, and C66) of the elastic tensor in units of N/m, Young's modulus E2D (N/m), shear modulus G (N/m), 2D layer modulus γ2D (N/m), bending modulus D (eV), the critical buckling strain ϵ c = −ϵ′c / L2 , and the Poisson ration ν, where L is the length of the nanosheet in Å. Note, C11 = (C11 + C22)/2 and technically, for hexagonal two-dimensional materials, C66 ≈ (C11 − C12)/2. b Defect is induced by removing one interfacial sulfur atom from the vdW heterostructure.

proceeding to the discussion of our main results, we present in Fig. 1(f)(k) the calculated electronic properties of the fully relaxed pristine monolayer MoS2 and BP. The obtained band gap of 1.70 and 0.90 eV for MoS2 and BP are in good agreement with previous results [16,17]. We also show in Table 1 their predicted effective masses along the kx and ky directions, lattice and elastic constants, which are also in good agreement with previous results [18]. Fig. 1(h & k) are only intend as a reference to demonstrate the topological variance of the surface states. Elastic properties. The elastic properties not only provide the information on the mechanical stability of a material but also serve as a guide on the potential device applications. In this regards, we have systematically determined the elastic properties of the vdW heterostructures and their parent materials. For a 2D material with a hexagonal crystal lattice, only the C11 ≈ C22, C12, and C66 of the elastic tensor Cij are the most relevant due to the lattice symmetry. As such, we can express the 2D in-plane elastic matrix as [19]

mitigate any negative properties of the individual materials and as such, enhance their potential applications. Herein, we studied vertically integrated hybrid structures formed by interchanging monolayers of MoS2 and BP up to ten layers. For brevity, the various stacks studied are between the parenthesis BP/MoS2 (BM), BP/MoS2/BP (B2M), BP/ MoS2/BP/MoS2 (B2M2), and (BP/MoS2)|1,..,(BP/MoS2)|5 (B5M5). Even number of stacks exhibit C3v (space group No. 156) symmetry while odd number of stacks exhibit D3h (space group No. 187) symmetry (Table 1). Our first-principles calculations reveal that the bandgap could be tuned from 0.89 to as low as 0.76 eV by increasing the number of interchanging BP and MoS2 monolayer stacks in the hybrid structure. In all cases, the bandgap remained direct at the K-point of the high symmetry zone. The predicted properties could be explored for applications in low-power devices. Our calculations further reveal that sulfur vacancy leads to bandgap narrowing but with in-gap impurity levels induced below the conduction band minimum. We show that the impurity states are predominantly uncompensated Mo-d states adjacent to the vacant site and the nearest-neighbor S site to the vacancy in the MoS2 layer with some density of pz states due to charge transfer from B and P, respectively. Such defect-induced levels could act as local scattering sites that will diminish electron-hole lifetimes leading to nonradiative recombination of photoexcited states, which could adversely affect device performance [11]. The defect-induced levels could as well be beneficial for intermediate band photovoltaics applications [12,13] whereby they could be delocalized via insulator-metal-transitions [13]. More intriguingly, the vdW hybrid heterostructures show superior mechanical properties. The 2D moduli of the vdW heterostructures are higher than the sum of the 2D modulus of the individual layers, which rather increased with increasing thickness. This behavior is in contrast to heterostructures of MoS2/graphene, MoS2/WS2 and few-layer graphene [14,15,8] that suffer from interlayer slippage that diminishes their mechanical strength as thickness is increased. Our data show that the superior properties of the vdW heterostructures could be due to the unique interlayer and interface coupling that abound between the stacks. This is further supported by the rather large frequency gap between the optical and acoustic phonon branches including the transformation of the flexural phonon mode and the emergence of shear and breathing phonon modes in the vibrational and phonon spectra of the vdW heterostructures. Our results not only provide a deep insight into understanding diverse properties of hybrid vdW heterostructure of monolayer MoS2 and BP but also the unique interlayer interactions in 2D van der Waals structures and potentially guiding the engineering of their nanoelectromechanical properties as desired.

C11 C12 0 ⎤ ⎡ ξ1 ⎤ σ ⎡ σ1 ⎤ ⎡ 0 ⎥ ⎢ ξ2 ⎥ C = ⎢ 2 12 C11 ⎢ ⎥ ⎥⎢ ⎥ 0 0 C ⎣ σ6 ⎦ ⎢ 66 ⎦ ⎢ ξ6 ⎥ ⎣ ⎣ ⎦

(1)

where σi and ξi is the stress and strain matrix, respectively. A polynomial fit of Eq. (1) led to the C11, C12, and C66 values listed in Table 1. Using the predicted independent elastic constants, we obtain the 2 2 − C12 )/ C11 and Poisson ratio ν = C12/C11, in-plane Young's E = (C11 shear modulus G = C66 (Table 1). The predicted elastic properties of monolayer MoS2 and BP are in good agreement with previous results [20–22,18,15]. For the hybrid structures, while the Poisson ratio is generally constant (changes between ± 0.10), our calculations reveal that the elastic constants increased rather linearly as a function of increasing number of layers reaching an in-plane Young's (Shear) modulus of 1388.61 (520.02 N/m) for B5M5. Our results show that the strength of the vdW heterostructures largely increased as the thickness is increased. We attribute the high elasticity to the very strong (covalent) bonding within the basal plane (mainly from BP) and the rather weak (van der Waals) bonding between the interlayer planes. The predicted high linear elasticity could be further enhanced by increasing the number of stacks to engineer a hybrid material with a desired mechanical property for applications in, e.g., structural composites, fibers, and protective coatings [23–25]. Such a superior mechanical property could also be explored for nanoelectromechanical device applications, such as in pressure sensing and fabrication of nanoresonators [26]. Our elastic constants data further show that the heterostructures do not suffer from interlayer slippage, which is an undesirable nanomechanical failure. Interlayer slippage occurs due to an abrupt loss of bonding in materials under strain and it has been reported in few-layer

2. Results and discussion Electronic properties of monolayer MoS2 and BP. Before 384

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

Fig. 1. Morphology of the van der Waals heterostructures. (a) The optimized lattice of monolayer MoS2 (above) and BP (below) showing the predicted lattice constants and the Mo-S and B-P bond lengths. hb denotes thickness. (b)-(e) The optimized crystal structure of the various stacking configurations of monolayers of MoS2 and BP plotted with the difference of the electron localization function. (f)–(h) The plot of the calculated electronic properties of monolayer MoS2 for the energy-lattice constant profile, band structure, and the surface states in that order. (i)–(k) The same plots as in (f)-(h) for BP. Observe the non-trivial topological nature of the surface states as shown in (h) and (k) for MoS2 and BP, respectively. 4π 2 D

strain ϵc = − 2D 2 , where hb and L (both in units of Å) are the thickness E L and length of the nanosheet, respectively. The layer modulus measures the resistance of a nanosheet to stretching [27]. Our data show that γ2D increases with increasing number of stacks in the vdW heterostructures (Table 1). On the other hand, the bending modulus defines the out-ofplane bending movement of atomically thin material [28,29]. The accurate determination of D depends sensitively on thickness hb, which in turn depends on the bending rigidity under deformation of the material. Due to this unavoidable uncertainty, D could vary over a wide range of values for a given material [28,30]. For our current study, we adopt the absolute thickness of the nanosheet, which provides a qualitative comparison of the lower limit of the bending moduli. For the monolayers of MoS2 and BP, we obtain hb of 3.13 and 0.84 Å, respectively. We note that hbMoS 2 is close to the value of ∼3.116 Å reported in Ref. [28,31] while hbBP is closer to the 0.80 Å reported for graphene, which shares the same planar structure [32]. For the vdW heterostructures, hb is estimated as the sum of hbMoS 2 and hbBP with

graphene, MoS2/graphene and MoS2/WS2 heterostructures [14,15,8]. We verify this for the vdW heterostructures by noting that the Young's modulus of the vdW heterostructures is always larger than the sum of those from individual components, i.e., EB2Di Mj > ∑i, j EB2Di Mj , where i and j are the number of sheets of BP and MoS2, respectively. For example, the E2D of B5M5 is greater than the sum of five stacks each of monolayer MoS2 and BP, respectively by more than 89.21 N/m. The absence of interlayer slippage in the vdW heterostructures could be due to the nature of the interlayer interactions despite their analogous structure to those of few-layer graphene, MoS2/graphene and MoS2/WS2 heterostructures [15]. Even in the presence of disorder due to interfacial sulfur vacancy, the vdW heterostructures still show superior mechanical properties. For example, the 2D Young's modulus and the associated elastic parameters only decreased by about 2% (Table 1). We further obtained the 2D layer modulus γ2D = (C11 + C12)/2, the bending or flexural modulus D =

2D 2 1 E hb , 12 1 − ν 2

and the critical buckling

385

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

the Brillouin zone center (Γ) could be expressed as [35–37] Γ = A1′ (R) + A2′ ′ (IR) + E′ (R + IR) + E′ ′ (R) , where R and IR denote Raman- and infrared-active mode, respectively. Since the primitive cell of monolayer MoS2 (BP) contains three (two) atoms, there will be a total of nine (six) normal vibrational modes at the Γ point as shown in Fig. 2(a & b). For monolayer MoS2, the three acoustic modes are characterized by transverse acoustic (TA), longitudinal acoustic (LA), and the out-of-plane transverse acoustic or flexural mode (ZA) while the optical branch could be grouped into two out-of-plane optical (ZO1/ ZO2), two longitudinal optical (LO1/LO2) and two transverse optical (TO1/TO2) branches. From the irreducible point group symmetry, each of A2′ ′ and E′ are acoustic modes, another A2′ ′ (ZO1) is IR-active, A1′ (ZO2) and E′ ′ (LO1/TO1) are R-active, and another E′ (LO2/TO2) are both R- and IR-active as denoted in Fig. 2(a). In the long wavelength regime, due to the slightly polar nature of MoS2, the relative displacement of Mo and S atoms induces intrinsic macroscopic electric field, which couples with the lattice leading to the LO-TO splitting of the IR-active phonon modes (E′) [38,39]. The predicted characteristic frequencies including the frequency gap between the optical and acoustic phonon branches ΔΩTO−LA ∼49.50 cm−1 are in good agreement with previous results [38,40,35,39,41]. The phonon spectra and the corresponding total and projected phonon density of states for BP is presented in Fig. 2(b) in good agreement with few reported phonon spectra of BP [20,18]. Below we note specific observations unique to our calculations. The ZA branch of the acoustic mode has a zero frequency while both the LA and TA modes are slightly finite with degenerate frequency at the zone center. In the optical branch, we have one ZO at ∼ 307 cm−1, one each of LO and TO at 965 and 983 cm−1. Applying the irreducible point group symmetry at the Γ, the A1′ (ZO) is R-active while E′ (TO/LO) is both Rand IR-active. The latter two optical modes E′ correspond to the inplane counter-phase motions of nearest-neighbor atoms. We expect characteristic features in the Raman spectra of BP within the above phonon frequencies. We also predicted a large ΔΩTO−LA ∼288.81 cm−1, which we attribute to the significant difference between the atomic mass of B and P (B:P∼10.8:30.97) and an interatomic bonding mediate almost by pure covalent bonding, which prohibits scattering across the LO-TA channel thus resulting in very long phonon relaxation time [42]. We present in Fig. 2(c & d) the phonon spectra and vibrational properties of two vdW hybrid structures: BM Fig. 2(c) and B2M Fig. 2(d). The phonon spectra and vibrational properties are important on many levels. Aside from providing insights for fundamental research, it is essential in the understanding of the observed mechanical, electronic, and thermal properties in the vdW heterostructures. Observe that the phonon and vibrational properties of the hybrid structures are distinctly different from either MoS2 or BP. This is due to the lattice dynamics induced by nearby layers especially in the low-frequency outof-plane modes, which we describe in details below. The vdW hybrid structures show relatively the same phonon spectra and vibrational density of states. As could be seen, the optical phonon modes are dominated by BP while the acoustic modes are formed by a strong hybridization between monolayer MoS2 and BP leading to a characteristic frequency gap between the optical and acoustic phonon branches of 288.81 and 291.30 cm−1 for BM and B2M, respectively in the pristine structures. The disordered spectra remains practically unchanged [Fig. 3(a & b)] aside from the renormalization of the characteristic frequency gap between the optical and acoustic phonon branches, which increased slightly to 299.20 and 298.15 cm−1 for BM and B2M, respectively. We attribute the rather large phonon bandgap to be partly due to the diversity in the masses of the basis atoms (Mo:S:B:P∼95.94:32.07:10.8:30.97) coupled with a strong interatomic interaction in the basal plane. Such a wide phonon bandgap means a higher quality factor of the resonator. Because of the large phonon frequency gap in the hybrid structures, the phase space for scattering

respect to the number of nanosheets and the number of vdW (interlayer) distances hbvdW [see Fig. 1(a)]. We obtain hbBM ∼ 7.56, hbB2M ∼ 12.0 , hbB2M2 ∼ 18.70 , and hbB5M5 ∼ 52.10 Å. The predicted bending modulus are listed in Table 1 . Observe that monolayer BP has a significantly smaller D ∼ 0.54 eV meaning that it is extremely soft in the out-of-plane direction probably due to the planar nature of the crystal structure. Our simulation reveal that the bending modulus significantly increased in monolayer MoS2 that could be rationalized from the trilayer atomic structure, which inhibits bending motion. Our predicted D ∼ 6.75 eV for MoS2 is in good agreement with previous results, which are in the range of 6.62-9.6 eV [33,28,30,31]. In the case of the pristine vdW heterostructures, we observe a systematic increase in the bending modulus reaching ∼ 22.73 keV in the B5M5 vdW heterostructure. This value is orders of magnitude higher than that of the parent materials. The larger D in the vdW heterostructures implies strong flexural rigidity. We ascribe the observed systematic increased in D in the vdW heterostructures to the increase in the number of layers as the stacking size is increased, which in turn, increases the strength of the intralayer covalent bonding inhibiting any bending motion. Using our predicted mechanical quantities, we could estimate the critical buckling strain for both the parent materials and the various vdW heterostructures. The critical buckling strain is an important parameter in structural design as it measures the maximum load above which a material suffers from a lateral deflection [34]. The calculated ϵc are presented in Table 1. Observe that for a sample of the same length L, our predicted ϵc for monolayer MoS2 is ∼14 times larger than that of monolayer BP while the largest vdW heterostructure we studied is 3853 (279) times larger than that of monolayer MoS2 (BP). We note that while monolayer BP has a larger Young's modulus than monolayer MoS2, the latter has a higher bending modulus and does not buckle easily under external compression. However, all the vdW heterostructures show higher Young's modulus with a corresponding bending modulus and critical buckling strain that is orders of magnitude larger than that of either monolayer MoS2 or BP, which highlights their superior mechanical stability that is essential in nanoscale devices due to their high surface-to-volume ratio. Most importantly, both the bending modulus and the critical buckling strain are practically unaffected by disorder. As could be seen from Table 1, D and ϵc decreased by less than 1%. This is significant as it demonstrates that the ubiquitous intrinsic native defects, e.g., sulfur vacancies prevalent in MoS2-based materials have insignificant impact on the structural and mechanical properties of MoS2-based vdW heterostructures. Phonon spectra and vibrational properties. As presented above, all the hybrid structures show a negative formation energy, which suggests that they are stable or metastable. The rather high elastic constants confirm their mechanical stability. In order to determine their dynamical and vibrational stability, we calculated their phonon spectra. An instability of the phonon mode in any K in the first Brillouin zone (BZ) will manifest as an imaginary frequency. Such phonon modes characterized by negative phonon modes cannot generate the needed restoring force in the event of lattice vibrations leading to the deviation of the material from its original configuration. In Fig. 2, we present the phonon spectra for the pristine structures along several high symmetry directions for the parent monolayer materials (MoS2 and BP) and two representative hybrid structures (BM and B2M) along with their total and projected phonon density of states. The corresponding phonon spectra and vibrational properties for the disordered vdW heterostructures are shown in Fig. 3. The phonon spectra show the absence of imaginary frequency in the whole BZ indicating that the hybrid structures are dynamically stable, and as a result could be experimentally synthesized if the necessary precursors are used with appropriate growth conditions. Both monolayers of MoS2 and BP exhibit D3h point group symmetry due to the presence of a horizontal reflection plane that passes through the anion. So, the irreducible representation of the vibrational modes at 386

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

Fig. 2. The phonon, vibrational, and thermal properties of pristine van der Waals heterostructures. Phonon spectra and corresponding density of states of (a) MoS2, (b) BP, (c) BM, (d) B2M. We highlight in the phonon spectra of the vdW heterostructures [Fig. 2(c & d)], the three acoustic modes (ZA, TA, and LA) shown with blue solid-line and the interlayer phonon modes, i.e., two shear modes (red solid-line) and N − 1 breathing mode (green solid-line). (e) Schematics of the breathing and shear interlayer phonon modes. (f) Heat capacity at constant volume C v = |∂E / ∂T|v . The horizontal dashed-line is the Dulong-Petit limit ∼3nNAκB.

Fig. 3. The phonon, vibrational, and thermal properties of disordered van der Waals heterostructures. Phonon spectra and corresponding density of states of (a) BM and (b) B2M. We highlight in the phonon spectra of the vdW heterostructures [Fig. 2(c & d)], the three acoustic modes (ZA, TA, and LA) shown with blue solid-line and the interlayer phonon modes, i.e., two shear modes (red solid-line) and N − 1 breathing mode (green solid-line). (c) Heat capacity at constant volume C v = |∂E / ∂T|v . The horizontal dashed-line is the Dulong-Petit limit ∼3nNAκB while the dashed-black line corresponds to the pristine C v of the same size as the disordered lattice. 387

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

and B2M, respectively. To ascertain the origin of the rather anomalous heat capacity, we did a calculation for the pristine sample of the same density (size) and compared both the density of vibrational states (not shown) and the specific heat [dashed-black line in Fig. 3(c)]. Aside from the broadening of some states, the heat capacity from the pristine and disordered vdW heterostructures are practically the same. Hence, the observed anomalous C v in the disordered vdW heterostructures is be due to the difference in the phonon density of states similar to the observations made in Ref. [51]. At low temperature, the specific heat data generally follow the Debye T3 law while at high temperature, it approaches the classical Dulong-Petit limit (3nNAκB), which is denoted with a horizontal dashedline in Fig. 2(f). The observed trend of C v in the vdW heterostructures is consistent with the superior thermal and mechanical properties discussed above. The large heat capacity coupled with the outstanding cycling stability in the MoS2-based heterostructure is promising for exploring potential applications as anode materials in batteries [52,53]. We could gain further insight into the thermal vibrations by computing the average acoustic Debye temperature ΘD, which is a measure of the temperature above which all modes begin to be excited and below which modes begin to be frozen out [54]. We calculated the 1 1 average acoustic Debye temperature as Θ3D = 3 ∑i Θ3 , where

across the optical-acoustic channel will be decreased due to loss of conservation of momentum [43,42]. As a result, one would expect weak phonon-phonon scattering that could lead to high thermal conductivity. Such observation has been reported in diamond and related materials [43,42]. The calculated phonon spectra of the vdW heterostructures show superior mechanical properties because of the transformation of the flexural modes and the appearance of shear and breathing modes [44,45]. These vibrational modes are due to the relative displacement of MoS2 sheet with respect to BP sheet and vice-versa in the perpendicular (out-of-plane) and tangential (in-plane) directions [Fig. 2(e)]. Materials that exhibit flexural modes are attractive because of their distinctive properties and their potential for novel applications, which seem to be a characteristic of layered 2DMs [46–48,29,45,49,50] that have been predicted in some multilayered 2DMs [50,46–48,29], where a set of different low-energy interlayer vibrational phonon modes associated with each layer thickness yields a rich spectrum of electronic properties distinct from the parent materials. Focusing on this observation, we could assign the three lowest frequency phonon branches to the ZA, TA, and LA acoustic modes, which are shown with blue solidline in Fig. 2(c & d) for the pristine and Fig. 3(a & b) for the disordered vdW heterostructures. Around the zone center, we also observe the interlayer phonon modes, i.e., two shear modes (red solid-line) that are blue-shifted with increasing N and N − 1 breathing mode (green solidline) that are red-shifted with increasing N, where N is the number of stacks. The shear modes are almost degenerate at the zone center in both the pristine and the disordered vdW heterostructures. We note that in the disordered BM vdW heterostructure, we only observed one each of the shear and breathing modes, respectively. The breathing phonon modes have been extensively studied because of their scientific and practical significance. They tend to be highly sensitive to interlayer vdW interactions and depending on their count N − 1, which generally increases with increasing number of stacks N, their zone center characteristic frequency could be IR-active or R-active [36,47]. We note that the presence of disorder did not induce any significant change to the phonon spectra and vibrational properties of the vdW heterostructures. From the projected phonon density of states [Fig. 2(c & d) for the pristine and Fig. 3(a & b) for the disordered vdW heterostructures], observe that the frequency scale of the optical (acoustic) phonon mode is dominated by atoms with smaller (large) masses similar to a diatomic linear chain model. Specifically, all the optical modes are derived from a strong hybridization of B and P atoms while the uppermost frequency regime of the acoustic mode around 575 cm−1 at Γ, which is composed mainly of the LA branch is exclusively from B atom. We also observed large phonon density of states of B atom in the ZA branch of the acoustic mode with a complex hybridization between Mo, S, B, and P for the TA and at the intermediate frequency region of the LA branch. Such large contribution from the flexural acoustic mode to the phonon density of states has been extensively studied and are reported to play a significant role in the observed high thermal conductance (see, e.g., Ref. [29,49]) and an almost ballistic transport [50] in some layered materials. Thermal properties. The isomeric heat capacity C v = |∂E / ∂T|v obtained from our phonon data are shown in Fig. 2(f) for the parent materials and the pristine structure of two of the vdW heterostructures. Our data show that the heat capacity of monolayer MoS2 is almost twice as large as that of monolayer BP. The predicted C v for monolayer MoS2 is in good agreement with reported values [40]. We are not aware of any reported C v for monolayer BP. The heterostructures exhibit very large heat capacity orders of magnitude greater than that of the constituent monolayers. At room temperature (300 K) we could extra a heat capacity of 63.10, 31.75, 379.56, and 506.84 J/K/mol for MoS2, BP, BM, and B2M, respectively. In the presence of disorder [Fig. 3(c)], the heat capacity of the vdW heterostructures increased dramatically, more than thrice in magnitude of the pristine counterpart. For instance, at room temperature, we obtain C v∼1.52 and 2.02 kJ/K/mol for BM,

Γ, i

ΘΓ,i = ℏ ΩΓ,i/κB and ΘΓ,i corresponds to the frequency at the zone boundary of the ith acoustic mode [55]. The calculated ΘD for MoS2, BP, BM, and B2M are 400.59, 586.18, 498.03, and 515.27 K, respectively. For the disordered vdW heterostructures, we obtained a slightly higher ΘD of 509.72 (529.49 K) for BM (B2M). The average acoustic Debye temperature seems to decrease with increasing average atomic mass. We, however, note that ΘD could as well be sensitive to the large frequency gap observed between the acoustic and optic phonons. The unique characteristics coupled with the multiple crossing phonon bands in the vdW heterostructures calls for a further investigation, e.g., calculation of the lattice thermal conductivity. Electronic properties and band-engineering in vdW heterostructures. In order to study the electronic structures, we present in Fig. 4(a)–(d) the calculated electronic band structure unfolded unto the first Brillouin zone of the hexagonal lattice and the corresponding density of states for the various pristine hybrid structures of BP and MoS2 heterostructures. Regarding the magnetic states of the pristine hybrid structures, our calculations yield spin unpolarized ground state with negligible spin-orbit interaction effects. All the various stacks of BP and MoS2 exhibit a direct band gap Eg < 1.0 eV. Specifically, we predict a band gap of 0.89, 0.85, 0.79, and 0.76 eV for BM, B2M, B2M2, and B5M5 hybrid structures, respectively. Comparing Fig. 4 with the band structure of pristine MoS2 [Fig. 1(g)] and BP [Fig. 1(k)], it is apparent that the spectra of the hybrid structures are a combination of that of monolayers of MoS2 and BP. In order to explore the origin of the observed large band gap renormalization, we checked the effects of lattice mismatch in the hybrid structures. At equilibrium geometry, the lattice constant of the 2 × 2 supercell is calculated to be 6.38 Å. This indicates that the lattice constant of monolayer MoS2 increased by 0.16% while that of monolayer BP decreased by 0.78%. This rather small strain in the lattice constant of MoS2 and BP monolayers only led to a ± 0.03 eV change in their band gaps. Such infinitesimal change is too small to explain such large renormalization in the band gap of the studied hybrid structures. We further investigated the origin of the significant reduction in the band gap of the hybrid structures by carrying out charge transfer calculations. In building the various pristine heterostructures, we observed an average charge transfer of 0.084|e|, 0.081|e|, 0.089|e|, and 0.091|e| per interlayer configuration for BM, B2M, B2M2, and B5M5 hybrid structures from BP to MoS2 layer. Such a charge transfer corresponds to a charge density in the neighbor of ∼1013e/cm2 leading to electrondoping in the hybrid structures. The BP monolayer bonds to the MoS2 monolayer through van der Waals interactions. At the equilibrium 388

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

Fig. 4. The electronic band structure of the pristine van der Waals heterostructures unfolded unto the primitive-cell Brillouin zone and the corresponding density of states for (a) BM, (b) B2M, (c) B2M2, and (d) B5M5. A direct bandgap of ∼0.89, 0.85, 0.79, and 0.76 eV in that order is predicted at the K-point of the reciprocal space. The horizontal dashed lines is the Fermi level denoted as EF is set at zero of the energy scale.

geometry, the MoS2 structure lies above the BP for the BM stacking and in-between the BP sheet for other hybrid structures at ≈ 3.59Å with no chemical bonding at the BP-MoS2 interface. We attribute the reduction in the predicted band gap and hence, changes in the spectra of MoS2 and BP hybrid structures to the observed charge transfer mediated by vdW interactions. The BP sheet acts as a passivation layer, which modified the local dielectric screening leading to a markedly lower energy bandgap. In order to obtain a more detailed information about the character of the states in the proximity of the bandgap, we explore the projected density of states (not shown). The valence band is formed by a strong hybridization of MoS2 and BP states. The density of states of the valence band maximum (VBM) is predominantly formed by hybridized Mo-(d z 2 & dxy) and pz-states of S, B, and P, respectively. The conduction band is dominated by Mo-d states. There is a significant contribution from the p-states of S, B, and P. Our data show that the conduction band minimum (CBM) exhibits a strong antibonding character of Mo-d z 2 and S-py states with some density of pz states from B and P, respectively. The character of the decomposed density of states is similar in the pristine hybrid structures studied herein. To gain insight into the transport properties of the hybrid structures, we calculated their carrier effective mass. The band effective mass mb is ℏ2 →T → obtained from the band structure by fitting a parabola Ek = 2m k A k 0 to the states around the band extremum (CBM and VBM), where k = (kx, ky) is the in-plane k-point measured from the band extremum, and m0 is the free electron mass. In Table 1, we show the calculated electron and hole effective masses along the kx (ky) direction. While mb are slightly different in the kx and ky directions, the transport could be said to be largely isotropic. Defect-engineering in vdW heterostructures. Next, we explore the effects of atomic defects in some of the hybrid structures. In our MoS2-BP hybrid example, we focus on sulfur vacancies, as the Mo vacancies could only be stable in specialized experimental conditions due to the formation of MoS3 [56]. The presence of sulfur vacancies has been extensively studied in monolayer MoS2 and characterized using various experimental, e.g., high-resolution transmission electron microscopy and computational approaches [56–60]. To model the sulfur

vacancies, we removed one sulfur atom from the interface of the MoS2 and BP. After optimizing the structure, the corresponding S on the other side of the S-Mo-S structure moved slightly towards the vacant site by ∼0.05Å with an average bond length of Mo-S in the proximity of the vacant site changing by about 0.85% due to ∼0.36% change in the lattice constant of MoS2. We also observed a slight shift from the next nearest neighbor sites in the direction of the vacant site. The average interlayer distance between the BP and MoS2 also increased by about 0.03 Å. With such small changes in the lattice parameters, we do not anticipate any localized electrons due to symmetry breaking. In the BP layer, the changes are relatively small with an average B-P bond length difference of less than 0.04% and the layer remained planar with no buckling. The structural relaxation of the disordered hybrid structures preserved the 3-fold symmetry of the pristine MoS2 [56]. The optimized defect hybrid structures exhibit C3v (space group No. 156) symmetry. A defect formation energy of 5.15 (3.33 eV) is obtained for the BM (B2M) vdW hybrid structure. We present in Fig. 5 the electronic properties of the defect hybrid structures. In this current study, we focus on the BM and B2M stacking configurations. Just as in the pristine hybrid structures, the magnetic properties of the disordered hybrid structures are largely nonmagnetic with a negligible local magnetic moment. As could be observed from the electron localization function [Fig. 5(a & b)], the electrons are quasi-localized within the defect site. Observe from the band structure and the corresponding density of states [Fig. 5(c & d)] that the defectinduced states manifest largely as in-gap levels leading to a significant decrease in the bandgap from 0.89 (0.85 eV) in the pristine structure for BM (B2M) to 0.44 (0.42 eV) in the disordered lattice (see also Fig. 6). Our calculations show significant charge transfer of [0.32|e| (2.30 × 1013e/cm2)] per interlayer configuration from BP to MoS2. To further quantify the defect-induced in-gap states, we show in Fig. 6 the partial density of states of the studied hybrid structures. In both the disordered BM [Fig. 6(a)] and B2M [Fig. 6(b)] hybrid structures, the effects of the disorder is the same. Aside from the defectinduced in-gap levels, as expected, there is the broadening of the spectra and the effect of the disorder is predominantly more pronounced around the Fermi level. Due to the missing S2− ion state, a positive Madelung potential is left in the vacant site thereby pushing 389

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

vacant site is extended beyond normal reaching a radius of about 6.40 Å [Fig. 5(a & b)]. In all the disordered hybrid structures studied, the bandgap remained direct at the K− point of the high symmetry zone. However, due to the defect-induced states, the antibonding interaction of Mo-d and S-p states decreased leading to the lowering of their energy. We further calculated the concentration of the in-gap defect states. The defect density of the in-gap levels Dv s is obtained by integrating the product of the number of in-gap impurity states and the Fermi-Dirac probability function f(E) as Dv s = ∫ Dv s (E ) d E to obtain 7.22 × 1011e/ cm2 and 7.19 × 1011e/cm2 for BM and B2M disordered hybrid structures, respectively. In order to quantify the effects of disorder to the transport properties of the vdW heterostructures, we calculated the band effective mass using the same procedure as in the pristine structures. The flat band around CBM (Fig. 5) is a signature of a heavy electron mass, especially along the ky direction. Our fit to the spectra led to electron effective mass of 3.39/3.32 and a hole effective mass of 0.16/0.17 in units of m0 along the kx/ky direction (Table 1). We propose that the defect-induced midgap levels could be explored for applications in photovoltaic devices (PVDs) as a novel material to capture solar energy. An intermediate-band solar cell is a novel photovoltaic device operating at more than 20% higher efficiency than conventional single-gap solar cells [12,13]. Current efforts in PVDs are mainly focused on increasing efficiency while reducing cost. The midgap impurity states in MoS2-based hybrids could be tuned to dramatically increase the efficiency of solar cells via two-phonon processes. In a two-photon process, two sub-bandgap photons are absorbed such that one photon pumps an electron from the valence band to the intermediate band, the other photon will pump from the intermediate band to the conduction band. Such an arrangement provides an extra channel to promote two low-energy photons instead of one from the valence to the conduction band leading to increased photocurrent without decreasing the photovoltage [12,13]. Unlike monolayer MoS2, MoS2-based heterostructure offers the flexibility to tune the bandgap and the midgap impurity levels not just with the disorder but also with the number of stacking layers. As noted above, midgap impurity levels could act as a local scattering center with an increased likelihood of nonradiative recombination of photoexcited electron-hole pairs [11]. One possible solution is to hyperdop the material beyond Mott-Anderson transition [62,63]. The increase in the density of impurity states delocalizes the electrons associated with the midgap levels via insulator-metal-transitions [13] while approaching the radiative limit. This has been demonstrated in several experiments for thin films and quantum-dots (see, e.g., Refs. [64–68]).

Fig. 5. The optimized lattice of (a) BM and (b) B2M hybrid structures with one interfacial sulfur vacancy. The lattices have been plotted with the electron localization function of the defective site. (c) and (d) show the corresponding band structure plotted with the density of states. The large defect-induced ingap states in the density of states, which are flat bands in the k− space (see the band structure) manifest as quasi-localized states around the defect-site.

3. Conclusion In summary, we present detailed first-principles calculations of the electronic and related properties of hybrid structures formed from vertically stacked monolayer of MoS2 and BP. The vibrational spectra of the vdW heterostructure reveals the presence of shear and breathing phonon modes, which enabled superior mechanical and elastic properties. The vdW heterostructures are characterized by large heat capacity and exhibit much better mechanical and elastic properties compared to that of the isolated consistent monolayers. The superior mechanical and elastic properties are shown to be largely preserved even in the presence of disorder due to interfacial sulfur vacancies in the vdW heterostructures. We demonstrate that the vdW heterostructure is a narrow bandgap semiconductor, which could be tuned by both number of stacking layers and defect-engineering. While defectengineering leads to bandgap narrowing, an in-gap impurity levels due mainly to uncompensated Mo-d states are induced, which could be explored for intermediate band solar cell applications. We note that due to the quasi-localized nature of the defect-induced state, electronelectron interactions may be important in understanding fully the practical implications for device applications. Also, the reported results

states in the neighborhood of the sulfide ions below the CBM by ∼ 0.25 eV and above the top of the light-hole valence band. This doubly occupied charged vacancy can trap electrons and form bound states resulting in the occupation of a hydrogenic-like state centered on the vacant site [Fig. 5(a & b)]. Taking the maximum of the density of defect states as a reference, the in-gap state is located at the center approximately 0.35 eV from the conduction band or from the valence band. Our analysis of the defect structures shows that it is comprised mainly of Mo−d states of four Mo atoms around the vacant site. This observation is in agreement with recent experiments and computations investigating native defects in monolayer MoS2 [56,59–61]. The defect states is due to additional states introduced by the uncompensated Mo atoms adjacent to the S vacancy and nearest S site to a sulfur vacancy in the MoS2 layer (Fig. 6). There is also a small contribution from B and P p-states due to charge transfer. Due to the subtle hybridization between Mo, S, B, and P states, the wavefunctions in the neighborhood of the 390

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

Fig. 6. The partial density of states of the Mo site adjacent to S vacancy, nearest S site to a S vacancy, B site, and P sites in (a) BM and (b) B2M vdW heterostructures (solid lines) compared to the partial density of states of the corresponding atom in the pristine crystals (denoted with a gray background). The defect-induced states are shown with an arrow.

One major challenge in designing a heterostructure-based material is the unintended strain from lattice mismatch. In our case, due to the small difference in the lattice constants of monolayer MoS2 and BP, the lattice mismatch between the two structures calculated as 2|(La − Lb)|/ (La + Lb) is < 1.0%, where La and Lb is the in-plane lattice constant in the supercell slabs of MoS2 and BP, respectively. The equilibrium geometry of the heterostructure has an average interlayer separation of ≈ 3.58 Å measured from the bottom of the BP monolayer to the S atom in the S-Mo-S stacking of MoS2. For the pristine structures, we used a 2 × 2 × 1 supercell slab. A larger 4 × 4 × 1 supercell practically led to the same results. The calculated stack normalized formation energy of the hybrid structures are −0.24, −0.31, and −0.40, and −0.42 eV for BM, B2M, B2M2, and B5M5, respectively. To model the defect hybrid structures, we used a 4 × 4 × 1 supercell slab with one missing sulfur atom from the interface of MoS2. The phonon spectra is obtained using density functional perturbation theory [73] as implemented in VASP in conjunction with the Phonopy package [74] using a 3 × 3 ×1 Monkhorst-Pack grid to represent the reciprocal space.

have used a semilocal functional, which may not accurately account for quantum effects in 2DMs. This could be improved by accounting for the many-body effects with first-principles-based many-body approaches. This will be research of future investigation.

4. Method To investigate the structural, electronic, elastic, and mechanical properties of monolayer MoS2 and BP vdW heterostructures, we carried out first-principles calculations using density functional theory (DFT) [69] as implemented in the Vienna Ab Initio Simulation Package (VASP) [70]. The calculations were performed using the Perdew-Burke-Ernzerhof (PBE) [71] exchange-correlation functional, a kinetic energy cutoff for the plane-wave basis set of 550 eV, and a 8 × 8 ×1 Monkhorst–Pack grid to represent the reciprocal space. All calculations included vdW interactions and dipole corrections to avoid spurious interactions between the periodically repeated images of the slab. We used a vacuum of at least 25Å along the out-of-plane direction to further eliminate the artifacts of the periodic boundary condition. All the structures were relaxed until the energy (charge) is converged to within ∼10−4 (10−8) eV and the forces dropped to ∼10−3 eV/ Å. The optimized crystal structure of monolayer MoS2 and BP are shown in Fig. 1(a). The lattice constant of the monolayers is obtained by fitting a polynomial to the energy-lattice constant profile [see Fig. 1(f) for MoS2 and Fig. 1(i) for BP]. The predicted equilibrium lattice constants are 3.183 and 3.213 Å for MoS2 and BP, respectively. The predicted equilibrium lattice constants are in good agreement with the reported values for the free-standing crystals in the range 3.16–3.20 and 3.18–3.22 Å for monolayer MoS2 and BP, respectively [20,72]. The representative crystal structures (plotted with the electron localization function difference) for the various heterostructures studied herein are shown in Fig. 1(b)–(d). We note that the bond length of Mo-S and B-P is 2.41 and 1.84, respectively in the isolated MoS2 and BP monolayer, respectively, which only changed by ± 0.38% after relaxing the heterostructures.

Acknowledgment The research was sponsored by the Army Research Laboratory (ARL) and was accomplished under the Cooperative Agreement Number W911NF-11-2-0030 as an ARL Research [George F. Adams] Fellow. Supercomputer support is provided by the DOD HighPerformance Computing Modernization Program at the Army Engineer Research and Development Center, Vicksburg, MS. References [1] a A.K. Geim, I.V. Grigorieva, Nature 499 (2013) 419; b T. Georgiou, R. Jalil, B.D. Belle, L. Britnell, R.V. Gorbachev, S.V. Morozov, Y.J. Kim, A. Gholinia, S.J. Haigh, O. Makarovsky, L. Eaves, L.A. Ponomarenko, A.K. Geim, K.S. Novoselov, A. Mishchenko, Nat. Nanotechnol. 8 (2012) 100; c A. Pospischil, M.M. Furchi, T. Mueller, Nat. Nanotechnol. 9 (2014) 257;

391

Materials Today Communications 19 (2019) 383–392

C.E. Ekuma, et al.

[2]

[3]

[4] [5]

[6] [7] [8]

[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

[20] [21] [22] [23] [24]

[25] [26]

[27]

125428. [28] J.-W. Jiang, Front. Phys. 10 (2015) 287–302. [29] J.-W. Jiang, B.-S. Wang, J.-S. Wang, H.S. Park, J. Phys.: Condens. Matter 27 (2015) 083001. [30] R.M. Elder, M.R. Neupane, T.L. Chantawansri, Appl. Phys. Lett. 107 (2015) 073101. [31] J.-W. Jiang, Z. Qi, H.S. Park, T. Rabczuk, Nanotechnology 24 (2013) 435705. [32] Y. Huang, J. Wu, K.C. Hwang, Phys. Rev. B 74 (2006) 245413. [33] G. Casillas, U. Santiago, H. Barrón, D. Alducin, A. Ponce, M. Jos-Yacamán, J. Phys. Chem. C 119 (2015) 710–715. [34] a S. Timoshenko, S. Woinowsky-Krieger, Theory of plates and shells, 2nd edn., McGraw-Hill, 1987; b C. Wang, L. Lan, H. Tan, Phys. Chem. Chem. Phys. 15 (2013) 2764–2773. [35] X. Zhang, X.-F. Qiao, W. Shi, J.-B. Wu, D.-S. Jiang, P.-H. Tan, Chem. Soc. Rev. 44 (2015) 2757–2785. [36] Y. Zhao, X. Luo, H. Li, J. Zhang, P.T. Araujo, C.K. Gan, J. Wu, H. Zhang, S.Y. Quek, M.S. Dresselhaus, Q. Xiong, Nano Lett. 13 (2013) 1007–1015. [37] J. Ribeiro-Soares, R.M. Almeida, E.B. Barros, P.T. Araujo, M.S. Dresselhaus, L.G. Cançado, A. Jorio, Phys. Rev. B 90 (2014) 115438. [38] A. Molina-Sánchez, L. Wirtz, Phys. Rev. B 84 (2011) 155413. [39] L. Matthes, P. Gori, O. Pulci, F. Bechstedt, Phys. Rev. B 87 (2013) 035438. [40] B. Peng, H. Zhang, H. Shao, Y. Xu, X. Zhang, H. Zhu, RSC Adv. 6 (2016) 5767–5773. [41] R. Yan, J.R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A.R. Hight Walker, H.G. Xing, ACS Nano 8 (2014) 986–993. [42] X. Gu, R. Yang, Appl. Phys. Lett. 105 (2014) 131903. [43] L. Lindsay, D.A. Broido, T.L. Reinecke, Phys. Rev. Lett. 111 (2013) 025901. [44] a P.H. Tan, W.P. Han, W.J. Zhao, Z.H. Wu, K. Chang, H. Wang, Y.F. Wang, N. Bonini, N. Marzari, N. Pugno, G. Savini, A. Lombardo, A.C. Ferrari, Nat. Mater. 11 (2012) 294; b S. Kitipornchai, X.Q. He, K.M. Liew, Phys. Rev. B 72 (2005) 075443; c K.H. Michel, B. Verberck, Phys. Rev. B 78 (2008) 085424; d N.J. Roome, J.D. Carey, ACS Appl. Mater. Interface 6 (2014) 7743–7750. [45] L. Liang, J. Zhang, B.G. Sumpter, Q.-H. Tan, P.-H. Tan, V. Meunier, ACS Nano 11 (2017) 11777–11802. [46] C.H. Lui, L.M. Malard, S. Kim, G. Lantz, F.E. Laverge, R. Saito, T.F. Heinz, Nano Lett. 12 (2012) 5539–5544. [47] C.H. Lui, T.F. Heinz, Phys. Rev. B 87 (2013) 121404. [48] M.-L. Lin, J.-B. Wu, X.-L. Liu, P.-H. Tan, J. Raman Spectrosc. 49 (2017) 19–30. [49] S. Shafraniuk, Thermoelectricity and Heat Transport in Graphene and Other 2D Nanomaterials, Elsevier Science, 2017. [50] N.B. Le, T.D. Huan, L.M. Woods, ACS Appl. Mater. Interface 8 (2016) 6286–6292. [51] A.I. Chumakov, G. Monaco, A. Fontana, A. Bosak, R.P. Hermann, D. Bessas, B. Wehinger, W.A. Crichton, M. Krisch, R. R“uffer, G. Baldi, G. Carini Jr., G. Carini, G. D’Angelo, E. Gilioli, G. Tripodo, M. Zanatta, B. Winkler, V. Milman, K. Refson, M.T. Dove, N. Dubrovinskaia, L. Dubrovinsky, R. Keding, Y.Z. Yue, Phys. Rev. Lett. 112 (2014) 025502. [52] Q.-F. Li, C.-G. Duan, X.G. Wan, J.-L. Kuo, J. Phys. Chem. C 119 (2015) 8662–8670. [53] X. Shao, K. Wang, R. Pang, X. Shi, J. Phys. Chem. C 119 (2015) 25860–25867. [54] a N. Ashcroft, N. Mermin, Solid State Physics, Holt, Rinehart and Winston, 1976; b T. Nakashima, Y. Umakoshi, Philos. Mag. Lett. 66 (1992) 317–321. [55] D.T. Morelli, J.P. Heremans, Appl. Phys. Lett. 81 (2002) 5126–5128. [56] W. Zhou, X. Zou, S. Najmaei, Z. Liu, Y. Shi, J. Kong, J. Lou, P.M. Ajayan, B.I. Yakobson, J.-C. Idrobo, Nano Lett. 13 (2013) 2615–2622. [57] H. Qiu, T. Xu, Z. Wang, W. Ren, H. Nan, Z. Ni, Q. Chen, S. Yuan, F. Miao, F. Song, G. Long, Y. Shi, L. Sun, J. Wang, X. Wang, Nat. Commun. 4 (2013) 2642. [58] a H.-P. Komsa, J. Kotakoski, S. Kurasch, O. Lehtinen, U. Kaiser, A.V. Krasheninnikov, Phys. Rev. Lett. 109 (2012) 035503; b J.-Y. Noh, H. Kim, Y.-S. Kim, Phys. Rev. B 89 (2014) 205417; c J.D. Fuhr, A. Saúl, J.O. Sofo, Phys. Rev. Lett. 92 (2004) 026802. [59] D. Liu, Y. Guo, L. Fang, J. Robertson, Appl. Phys. Lett. 103 (2013) 183113. [60] J. Hong, Z. Hu, M. Probert, K. Li, D. Lv, X. Yang, L. Gu, N. Mao, Q. Feng, L. Xie, J. Zhang, D. Wu, Z. Zhang, C. Jin, W. Ji, X. Zhang, J. Yuan, Z. Zhang, Nat. Commun. 6 (2015) 6293. [61] D. Le, T.B. Rawal, T.S. Rahman, J. Phys. Chem. C 118 (2014) 5346–5351. [62] N.F. Mott, Rev. Mod. Phys. 40 (1968) 677–683. [63] D. Belitz, T.R. Kirkpatrick, Rev. Mod. Phys. 66 (1994) 261. [64] Y. Zhou, F. Liu, M. Zhu, X. Song, P. Zhang, Appl. Phys. Lett. 102 (2013) 222106. [65] E. Antolín, A. Martí, J. Olea, D. Pastor, G. González-Díaz, I. Mártil, A. Luque, Appl. Phys. Lett. 94 (2009) 042115. [66] M.-J. Sher, E. Mazur, Appl. Phys. Lett. 105 (2014) 032103. [67] E. Antolín, A. Martí, C.D. Farmer, P.G. Linares, E. Hernández, A.M. Sánchez, T. Ben, S.I. Molina, C.R. Stanley, A. Luque, J. Appl. Phys. 108 (2010) 064513. [68] A. Martí, D.F. Marr’øn, A. Luque, J. Appl. Phys. 103 (2008) 073706. [69] a P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864–B871; b W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133–A1138. [70] G. Kresse, J. Furthm”uller, Comput. Mater. Sci. 6 (1996) 15–50. [71] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [72] a M.-H. Chiu, C. Zhang, H.-W. Shiu, C.-P. Chuu, C.-H. Chen, C.-Y.S. Chang, C.H. Chen, M.-Y. Chou, C.-K. Shih, L.-J. Li, Nat. Commun. 6 (2015) 7666; b F.A. Rasmussen, K.S. Thygesen, J. Phys. Chem. C 119 (2015) 13169–13183; c C. Ataca, H. Şahin, S. Ciraci, J. Phys. Chem. C 116 (2012) 8983–8999. [73] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, F. Bechstedt, Phys. Rev. B 73 (2006) 045112. [74] A. Togo, F. Oba, I. Tanaka, Phys. Rev. B 78 (2008) 134106.

d W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M.F. Crommie, A. Zettl, Appl. Phys. Lett. 98 (2011) 242105. a H. Terrones, F. López-Urías, M. Terrones, Sci. Rep. 3 (2013) 1549; b A. Ebnonnasir, B. Narayanan, S. Kodambaka, C.V. Ciobanu, Appl. Phys. Lett. 105 (2014) 031603; c Y.-L. Sun, D. Xie, M.-X. Sun, C.-J. Teng, L. Qian, R.-S. Chen, L. Xiang, T.-L. Ren, Sci. Rep. 8 (2018) 5107; d G. Nazir, M.F. Khan, S. Aftab, A.M. Afzal, G. Dastgeer, M.A. Rehman, Y. Seo, J. Eom, Nanomaterials (2018) 8. a K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666–669; b W.-Y. Shan, J. Zhou, D. Xiao, Phys. Rev. B 91 (2015) 035402; c A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys. 81 (2009) 109; d H. Zeng, C. Zhi, Z. Zhang, X. Wei, X. Wang, W. Guo, Y. Bando, D. Golberg, Nano Lett. 10 (2010) 5049; e L. Song, L. Ci, H. Lu, P.B. Sorokin, C. Jin, J. Ni, A.G. Kvashnin, D.G. Kvashnin, J. Lou, B.I. Yakobson, P.M. Ajayan, Nano Lett. 10 (2010) 3209; f J.C. Meyer, A. Chuvilin, G. Algara-Siller, J. Biskupek, U. Kaiser, Nano Lett. 9 (2009) 2683–2689; g B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, A. Kis, Nat Nano 6 (2011) 147–150; h O. Lopez-Sanchez, D. Lembke, M. Kayci, A. Radenovic, A. Kis, Nat Nano 8 (2013) 497–501; i D. Akinwande, N. Petrone, J. Hone, Nat. Commun. 5 (2014) 5678. a C.E. Ekuma, V. Dobrosavljević, D. Gunlycke, Phys. Rev. Lett. 118 (2017) 106404; b C.E. Ekuma, D. Gunlycke, Phys. Rev. B 97 (2018) 201414. a M. Farmanbar, G. Brocks, Phys. Rev. B 91 (2015) 161304; b B. Liu, Y.-Q. Zhao, Z.-L. Yu, L.-Z. Wang, M.-Q. Cai, J. Colloid Interface Sci. 513 (2018) 677–683. M.-L. Tsai, S.-H. Su, J.-K. Chang, D.-S. Tsai, C.-H. Chen, C.-I. Wu, L.-J. Li, L.-J. Chen, J.-H. He, ACS Nano 8 (2014) 8317–8322. W. Zhang, C.-P. Chuu, J.-K. Huang, C.-H. Chen, M.-L. Tsai, Y.-H. Chang, C.-T. Liang, Y.-Z. Chen, Y.-L. Chueh, J.-H. He, M.-Y. Chou, L.-J. Li, Sci. Rep. 4 (2014) 3826. A. Falin, Q. Cai, E.J.G. Santos, D. Scullion, D. Qian, R. Zhang, Z. Yang, S. Huang, K. Watanabe, T. Taniguchi, M.R. Barnett, Y. Chen, R.S. Ruoff, L.H. Li, Nat. Commun. 8 (2017) 15815. C.R. Dean, A.F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K.L. Shepard, J. Hone, Nat. Nanotechnol. 5 (2010) 722. K. Ren, M. Sun, Y. Luo, S. Wang, J. Yu, W. Tang, Appl. Surf. Sci. 476 (2019) 70. a W. Shockley, W.T. Read, Phys. Rev. 87 (1952) 835; b R.N. Hall, Phys. Rev. 87 (1952) 387. a A. Luque, A. Martí, Phys. Rev. Lett. 78 (1997) 5014–5017; b A. Luque, A. Marti, C. Stanley, Nat Photon 6 (2012) 146–152. A. Luque, A. Martí, E. Antolín, C. Tablero, Physica B 382 (2006) 320–327. X. Wei, Z. Meng, L. Ruiz, W. Xia, C. Lee, J.W. Kysar, J.C. Hone, S. Keten, H.D. Espinosa, ACS Nano 10 (2016) 1820–1828. K. Liu, Q. Yan, M. Chen, W. Fan, Y. Sun, J. Suh, D. Fu, S. Lee, J. Zhou, S. Tongay, J. Ji, J.B. Neaton, J. Wu, Nano Lett. 14 (2014) 5097–5103. N. Lu, H. Guo, L. Li, J. Dai, L. Wang, W.-N. Mei, X. Wu, X.C. Zeng, Nanoscale 6 (2014) 2879–2886. J. Dong, H. Li, L. Li, Npg Asia Materials 5 (2013) e56. D. Çakr, D. Kecik, H. Sahin, E. Durgun, F.M. Peeters, Phys. Chem. Chem. Phys. 17 (2015) 13013–13020. a X. Wei, B. Fragneaud, C.A. Marianetti, J.W. Kysar, Phys. Rev. B 80 (2009) 205407; b J. Zhou, R. Huang, J. Mech. Phys. Solids 56 (2008) 1609–1623. H. Şahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R.T. Senger, S. Ciraci, Phys. Rev. B 80 (2009) 155453. R.C. Cooper, C. Lee, C.A. Marianetti, X. Wei, J. Hone, J.W. Kysar, Phys. Rev. B 87 (2013) 035423. S. Bertolazzi, J. Brivio, A. Kis, ACS Nano 5 (2011) 9703–9709. K.S. Novoselov, V.I. Fal’ko, L. Colombo, P.R. Gellert, M.G. Schwab, K. Kim, Nature 490 (2012) 192–200. A.C. Ferrari, F. Bonaccorso, V. Fal’ko, K.S. Novoselov, S. Roche, P. Bggild, S. Borini, F.H.L. Koppens, V. Palermo, N. Pugno, J.A. Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis, J. Kivioja, C. Marinelli, T. Ryhnen, A. Morpurgo, J.N. Coleman, V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G.F. Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A.N. Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen, A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini, A. Tredicucci, G.M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath, B.J. van Wees, H. van der Zant, L. Occhipinti, A. Di Matteo, I.A. Kinloch, T. Seyller, E. Quesnel, X. Feng, K. Teo, N. Rupesinghe, P. Hakonen, S.R.T. Neil, Q. Tannock, T. Lfwander, J. Kinaret, Nanoscale 7 (2015) 4598–4810. K. Hu, D.D. Kulkarni, I. Choi, V.V. Tsukruk, Prog. Polym. Sci. 39 (2014) 1934–1972. a Q. Pei, Y. Zhang, V. Shenoy, Carbon 48 (2010) 898–904; b H. Zhu, Y. Wang, J. Xiao, M. Liu, S. Xiong, Z.J. Wong, Z. Ye, Y. Ye, X. Yin, X. Zhang, Nat. Nanotechnol. 10 (2014) 151; c S.C. Masmanidis, R.B. Karabalin, I. De Vlaminck, G. Borghs, M.R. Freeman, M.L. Roukes, Science 317 (2007) 780–783; d K.A. Marx, Biomacromolecules 4 (2003) 1099–1120. R.C. Andrew, R.E. Mapasha, A.M. Ukpong, N. Chetty, Phys. Rev. B 85 (2012)

392