Weak ferroelectric charge transfer in layer-asymmetric bilayers of 2D semiconductors

Abstract

In bilayers of two-dimensional semiconductors with stacking arrangements which lack inversion symmetry charge transfer between the layers due to layer-asymmetric interband hybridisation can generate a potential difference between the layers. We analyse bilayers of transition metal dichalcogenides (TMDs)—in particular, $\hbox {WSe}_2$—for which we find a substantial stacking-dependent charge transfer, and InSe, for which the charge transfer is found to be negligibly small. The information obtained about TMDs is then used to map potentials generated by the interlayer charge transfer across the moiré superlattice in twistronic bilayers.

Introduction

Many two-dimensional (2D) materials¹ lack inversion symmetry in their stoichiometric monolayers. These include hexagonal boron nitride (hBN), all TMDs, such as $\hbox {WSe}_2$, $\hbox {MoSe}_2$, etc., and the post-transition metal chalcogenides (InSe, GaSe). Depending on the orientation of the unit cells, bilayers of these materials can have inversion symmetry restored (like in 2H TMDs) or absent, like in $\gamma$-InSe. The orientation of the unit cells can therefore play an important role in determining electronic properties of twistronic structures of 2D semiconductors, where local stacking of the layers varies following the moiré superlattice structure. It is possible for layer-asymmetric structures to develop interlayer hybridisation between occupied (valence) and unoccupied (conduction) monolayer bands which leads to a net charge transfer between the layers^{2,3,4,5,6,7,8}, resulting in a vertical electric field piercing the 2D crystal. In particular, for heterobilayers and nearly parallel (P-stacked) homobilayers of TMDs the charge transfer and resulting potential difference between the two layers will determine the features of band edge states of twisted bilayers.

Here, we focus on the theoretical modelling aspect of ferroelectric charge transfer in 2D materials. The potential variation across a layer-asymmetric 2D material system in which charge transfer has occurred is a natural subject for first principles density functional theory (DFT) calculations. In modelling a 2D system using a plane-wave DFT code, to satisfy the requirement for periodicity in the out-of-plane direction, it is typical to construct a supercell such that repeated images of the few-layer crystal are separated by a large vacuum, minimising the interaction between them. Charge transfer will give the electrostatic vacuum potentials on either side of the 2D system different values, violating the requirement for out-of-plane periodicity in the DFT calculations. In this study, we compare results which do not correct for the effect of a polar bilayer with calculations using two means by which periodicity is often satisfied—first the construction of a supercell containing two images of the system, one with the layers interchanged, second the approximate method of applying a compensating step potential in the vacuum region—a ‘dipole correction’^9,10.

This article is structured as follows. First we give a full account of the DFT methods used in this work, with $\hbox {WSe}_2$ employed as a test system. Next, we explore charge transfer effects across a range of semiconducting TMDs. Then, we analyse InSe bilayers and show that they exhibit much weaker interlayer charge transfer as compared to $\hbox {WSe}_2$, $\hbox {WS}_2$, $\hbox {MoSe}_2$ and $\hbox {MoS}_2$. Finally, we implement information collected about stacking-dependent interlayer charge transfer in TMD bilayers to discuss its manifestation in the domain structure of twistronic TMDs.

DFT calculations of interlayer charge transfer in semiconductor bilayers

Achieving out-of-plane periodicity

Density functional theory (DFT) calculations of ultrathin films of layered materials using plane-wave based methods require the resolution of a crystal which is only periodic into two dimensions into a three-dimensionally periodic system. This is usually achieved by the construction of a supercell in which images of the two-dimensional (2D) layers are repeated periodically along the third dimension, with a large vacuum between them to ensure the layers are isolated from each other. Here, we compare results for polar $\hbox {WSe}_2$ bilayers using three commonly-used methods of addressing the requirement for periodicity in first-principles calculations of 2D materials: (1) a supercell with a single bilayer⁷, (2) a single bilayer supercell, but with a dipole correction applied¹¹, and (3) a supercell with two mirror-reflected images of the bilayer³.

We consider these methods in XM$'$-stacked $\hbox {WSe}_2$ (the prime symbol indicating that of the vertical metal–chalcogen pair in the bilayer, the metal atom is in the top layer¹²), the structure of which we show in Fig. 1. In the upper panel of Fig. 2, we show its plane-averaged local electrostatic potential (ionic and Hartree contributions) relative to the vacuum level on the Se side of the vertical W–Se pair, which is set to 0 eV. While the detailed variation of the local potential in the vicinity of the atomic planes depends sensitively on how the local and non-local parts of the pseudopotentials are set up, an important meaningful quantity can be extracted from the difference between the vacuum levels on both sides. As modelled in a previous work⁸, layer-asymmetric hybridisation between occupied and unoccupied states gives rise to charge transfer between the layers, giving the bilayer a finite dipole moment which results in the vacuum potentials on either side of the slab having different values. As shown schematically as an inset, the calculation was carried out using a supercell containing two mirror images of the bilayer (the second image thus being MX$'$-stacked), separated by large vacuums. Then, the requirement for out-of-plane periodicity in a plane wave DFT code is met, with the vacuum potentials matching at the supercell boundary.

Figure 1

Upper and lower panels: top and side view, respectively, of XM$'$-stacked $\hbox {WSe}_2$ (the prime indicating that of the vertical metal(M)–chalcogen(X) pair, the metal atom is in the top layer). The solid and dashed vertical lines in the upper panel emphasise the symmetry-breaking of this stacking, with the metal–chalcogen pair vertically opposite in one direction, but not in the other. The italicised labels in the upper panel ($M_1, X_{12}$ etc.) indicate the atomic positions referred to in Table 3.

Figure 2

Upper panel: Out-of-plane dependence of in-plane averaged electron potential energy of XM$'$ bilayer $\hbox {WSe}_2$ (including ionic and Hartree contributions), relative to the vacuum potential on the Se side of the vertical W–Se pair, $U_{vac}(Se)$, which is set to 0 eV. The calculation is for a double-bilayer supercell, with the structure shown as a schematic inset. The charge transfer between the layers gives each bilayer a finite dipole moment, with a consequent difference ($\Delta ^P$) between the vacuum levels on either side of a bilayer. Left-hand lower panel: potentials with isolated monolayer contributions subtracted, to better show the potential drop and other features. Three calculation methods are shown: black line—first 30 Å of the upper panel, red line—using only a single bilayer supercell, showing how the mismatched vacuum levels give a finite displacement field as necessitated by the periodic boundary conditions, blue line—a single bilayer supercell, but with a compensating dipole correction applied. The right-hand lower panel shows the same quantities, calculated using the LDA for comparison. 0 meV is set to the double-supercell vacuum level on the Se side of the vertical W–Se pair.

In the left-hand lower panel of Fig. 2, we compare the results for the first bilayer of the double-bilayer supercell with two methods from a supercell containing a single bilayer. To more easily see the effects of interlayer charge transfer, we subtract the potentials of isolated monolayers from the bilayer potentials. In a supercell with only one bilayer, the periodicity in the out-of-plane direction is broken, and there is a mismatch between the vacuum potentials at the supercell boundary. The numerical result of this in the DFT calculation is an artificial displacement field experienced across the slab, as can be seen in the red line in Fig. 2. The difference between the single- and double-bilayer supercells is greatest in the vacuum regions, as the effect of the artificial field is partially suppressed by the dielectric response within the $\hbox {WSe}_2$ bilayer itself.

We also show a further set of results for a single supercell, but with the well-known dipole correction applied^9,10. This adds a background sawtooth potential at each electronic self-consistency step of the calculation, with its size calculated from the electric dipole moment of the bilayer. As can be seen from the agreement between the black and blue lines in Fig. 2, this is a very good approximation to the truly periodic double-bilayer, while computationally much cheaper, being a supercell of half the size. In Table 1, we show a few key energies from the band structures resulting from the various supercell methods. While the single- and double-bilayer supercell methods give slightly different results, the application of a dipole correction to the single-bilayer supercell returns the values to good agreement with those for the double-bilayer supercell.

Table 1 Quantities (meV) calculated using different supercell methods for XM$'$ bilayer WSe₂—2×(1×).

Choosing DFT functional

For comparison with the PBE GGA results presented above, calculations were also carried out using the local density approximation (LDA), via the exchange correlation functional of Ceperley and Alder¹³ as parametrized by Perdew and Zunger¹⁴, and in Fig. 2 we show the differences between bilayer and monolayer plane-averaged local electrostatic potentials using LDA alongside the PBE results. In Table 1 the LDA-calculated potential differences and some band energies are compared with the PBE results. The bilayer-monolayer potential differences show a drop across the bilayer similar to that seen using PBE, but with a notable peak in between the layers. The potential difference across the bilayer, and the splitting between the upper K-point valence bands show small differences of only a few meV, but the difference between the $\Gamma$- and K-point valence band energies is notably reduced on going from the PBE approximation to LDA. Since the pseudopotential configurations as shown in Table 2 are nearly identical, we ascribe the differences between the calculations to the approximation to the exchange-correlation potential used.

Table 2 Details of projector augmented wave (PAW) pseudopotentials used in VASP calculations of $\hbox {WSe}_2$ and InSe.

Configuration dependence of weak ferroelectric effect in P-stacked bilayers

We now study various TMD bilayers with XM$'$/MX$'$ configuration as well as a number of various stackings which will enable us to describe the variation of charge transfer across the moiré supercell of a twisted bilayer. The strength of hybridisation between the layers is sensitive to the sublattice composition of the band states, and it is stronger for states residing on the chalcogen sublayers. In Table 3 we compare (with calculations using the QE code and a double-bilayer supercell) the wavefunction projections onto the six atomic layers of MX$'$ stacked TMD bilayers. As noted previously for $\hbox {MoSe}_2$ bilayers⁷, but repeated in a manner common to all four TMDs considered, the K-point wavefunctions are nearly entirely layer-polarised, due to a combination of very weak interlayer intraband hybridisation and the effective electric field between the layers arising from the charge transfer effect discussed above. In contrast, the $\Gamma$-point valence band wavefunction has an almost zero out-of-plane dipole moment, due to the strong interlayer hybridisation of the $\Gamma$-point states. The Q-point, which in some cases forms the conduction band minimum, is an intermediate case.

Table 3 Wavefunction projections onto atomic layers for XM$'$ stacked TMD bilayers, for atomic positions shown in Fig. 1, together with the dipole moment $d_z = e\langle {\psi }\vert z\vert {\psi }\rangle$, for valence (VB) and conduction (CB) states at important points in the Brillouin zone ($\hbox {Q}=\hbox {K}/2$).

The variation of the potential drop across a $\hbox {P-MX}_2$ bilayer with in-plane offset, $\varvec{r}_0$ ($\varvec{r}_0=0$ for XX-stacking corresponding to overlapping of chalcogens in two layers), and interlayer distance, d, can be described using the following expression⁸,

$$\begin{aligned} \Delta ^{P}(\varvec{r}_0,d) = \Delta _a e^{-q(d-d_0)}\sum _{j=1,2,3}\sin (\varvec{G}_j\cdot \varvec{r_0}), \end{aligned}$$

(1)

where values of parameters $\Delta _a$, q and $d_0$, $\varvec{G}_{1,2,3}$ are the shortest reciprocal vectors of a monolayer related by $120^{\circ }$-rotation around z. This formula is applicable to all TMDs with a honeycomb lattice structure, with parameters for $\hbox {MoS}_2$, $\hbox {MoSe}_2$, $\hbox {WS}_2$ and $\hbox {WSe}_2$ calculated using the QE code, shown in Table 4. Since in twisted bilayers the interlayer distance d and in-plane offset $\mathbf {r}_0$ vary continuously, we will use the information presented here concerning the dependence of charge transfer on stacking configuration to map the charge transfer and on-layer potential in a moiré superlattice.

Table 4 Vacuum energy difference across MX$'$ TMD bilayers, $\Delta ^P$(MX$'$), together with the parametrisation of the general configuration dependence of the potential drop, $\Delta ^P(\mathbf {r}_0,d)$, Eq. (1), calculated using the QE code.

Indium selenide

As a comparison to the transition metal dichalcogenides, we consider indium selenide, a member of the family of post-transition metal chalcogenides. The two most commonly found polytypes of bulk InSe in experiments are the $\gamma$¹⁵ and $\varepsilon$¹⁶ polytypes—on exfoliation to a bilayer, these will both have the same layer-asymmetric MX$'$/XM$'$ character to their stacking order. In the two panels of Fig. 3 we show first the local electrostatic potentials with isolated monolayer contributions subtracted comparing the supercell and dipole correction methods, and in the second the same differences, but calculated using LDA. For InSe, the peak in the difference between monolayer and bilayer potentials in the interlayer region shown in the LDA results for $\hbox {WSe}_2$ is present for both PBE and LDA approximations—but is much greater in magnitude for the LDA case. The charge transfer for bilayer InSe is negligible, giving a difference of only $\sim$2 meV between the vacuum potentials across the bilayer, with the consequence that differences between supercell and correction methods are negligible.

Figure 3

Left panel: difference between bilayer and isolated monolayer plane-averaged local potentials for 3R-stacked bilayer InSe, comparison between supercell methods. Right panel: same as left panel, but calculated using LDA. 0 meV is set to the double-supercell vacuum level on the Se side of the vertical In-Se pair.

Mapping charge transfer across the moiré superlattice of twistronic bilayers

While the asymmetry of band-edge states in inversion asymmetric bilayers manifests itself in the linear Stark shift of band energies and of the energies of optical transitions⁷ in vertically biased MX$'$ and XM$'$ bilayers, the ferroelectric potentials can be detected by contrast in the potential maps of systems with laterally varying stacking. Such a variation naturally appears in twisted bilayers with a parallel orientation of monolayer unit cells, where the twist angle determines the period of recurrent stacking configurations, known as moiré superlattice. To describe such variation, we employ the description of interlayer potential and the related size of the double layer of charge, $\pm \delta \rho$, on the top/bottom layers.

To demonstrate the latter, in Fig. 4 we show the z-coordinate dependence of the difference between the plane-averaged charge density of XM$'$ stacked bilayer $\hbox {WSe}_2$ and that of two isolated $\hbox {WSe}_2$ monolayers. The greatest differences and charge transfer are to be found in the interlayer region, with a peak and a trough close to the inner Se atomic layers. This peak(trough) could be related to the (de)population of hybridised s and $p_z$ orbitals on the Se atoms. We can calculate a charge transfer density directly from DFT as

$$\begin{aligned} \delta \rho = \frac{1}{2}\left( \int _0^{15}\rho (z)dz - \int _{15}^{30}\rho (z)dz\right) , \end{aligned}$$

(2)

where $\rho (z)$ is the plane-averaged charge density at the z-coordinates as shown in Figs. 2 and 4, and $z=15$ Å  is the mean plane between the two $\hbox {WSe}_2$ layers. This gives $\delta \rho = 1.9\times 10^{12}\,\hbox {cm}^{-2}$.

Figure 4

Difference between plane-averaged charge density of an XM$'$ bilayer and two isolated monolayers, showing greatest charge transfer in interlayer region between planes of Se atoms. Dashed line is a guide to the eye, showing $\Delta \rho = 0$.

We can also roughly estimate the magnitude of charge density transferred between the layers based on the potential drop, $\Delta ^P$ [Eq. (1)], as $\dfrac{\varepsilon _0\Delta ^P}{e^2d_{XX}}$ where $d_{XX}$ is the distance between the inner chalcogen atomic planes (3.14 Å for $\hbox {WSe}_2$). This gives $\sim 10^{12}$ cm$^{-2}$ for MX$'$/XM$'$ stacked bilayer $\hbox {WSe}_2$, with the net transfer of electrons being to the layer in which the selenium atoms are vertically opposite the tungsten atoms in the other layer. The difference in the precise numerical values between that directly calculated from DFT, and that from the rough estimate, arises as the greatest part of the charge transfer occurs in the interlayer region, over a smaller distance than the $d_{XX}$ used to convert from the potential drop.

In a rigid twisted bilayer (which corresponds to $\theta _P > 2.5^{\circ }$¹²) the relative in-plane shift of the layers is $\varvec{r}_0 = \theta \hat{z}\times \varvec{r}$. In a reconstructed twistronic bilayer of a marginally twisted ($\theta _P\ll 1^{\circ }$) parallel(P)-stacked TMD^12,17, a pattern of triangular domains forms, with alternating MX$'$/XM$'$ stacking, with $\varvec{r}_0 = \theta \hat{z}\times \varvec{r} + \varvec{u}(\varvec{r})$, where $\varvec{u}(\varvec{r})$ is the relative displacement field of the two layers, formed on reconstruction. Within this domain pattern, there will be an excess of bonded electron charges in the top layer in the centre of MX$'$ domains, and a corresponding deficiency in XM$'$ domains, with a general $\varvec{r}_0,d$ dependence in the resulting potential distribution given by Eq. (1). An estimate of the magnitude of the charge density transferred between the layers can be roughly approximated as set out above.

The distribution of potential above a twistronic bilayer resulting from its out-of-plane polarisation can be mapped using scanning Kelvin probe microscopy¹⁸ (SKPM) or a single electron transistor^19,20. The SKPM signal would vary between MX$'$ and XM$'$ regions, with the magnitude of variation dependent on the distance from the scanning tip to the surface, and between the bilayer and metallic back plate, which provides the reference for the locally measured potential. For a structure with a thick dielectric substrate separating the bilayer from the back plate by more than the superlattice period, as shown schematically in Fig. 5a, the potential measured by the tip close to the bilayer surface would display a variation with amplitude $\Delta ^P(\mathrm{MX'})$, as shown on the map. For a structure where the back plate is at a distance from the bilayer much less than the moiré superlattice period (Fig. 5b), the potential variation between XM$'$ and MX$'$ domains would be $2\Delta ^P(\mathrm{MX'})$. The corresponding values of $\Delta ^P$ are listed in Table 4. Subject to the requirement that the bilayer remains undoped (so that lateral potential screening inside the bilayer would not kill the effect), the described behaviour should be expected in all twisted TMD bilayers discussed in this work, as well as in heterobilayers⁸.

Figure 5

(a) Left panel: schematic of Kelvin probe setup where TMD bilayer is placed on dielectric substrate with thickness much larger than moiré period $\propto \ell$. Potential differences are then measured w.r.t. the middle plane of the TMD bilayer, with the resulting patterns of potential exemplified in right panel: map of top layer potential of a nearly parallel-stacked $\hbox {WSe}_2$ bilayer with twist angle $\theta =0.6^{\circ }$, from sum of piezoelectric potential¹² contribution and the charge transfer described in this work. Inset: shape of potential drop on going from an XM$'$ domain to MX$'$ (path marked with an arrow in the map), with an in-plane electric field at the domain wall. (b) Left panel: for a twisted bilayer sample placed separated from a metallic plate by only a very thin dielectric of thickness much less than the moiré superlattice period, the potential in the bottom layer will be the same for all domains. This will double the potential difference measured across a domain wall in the top layer, shown in the right panel. The drop of potential occurs on the length $\lambda \approx 8$ nm corresponding to the width of the domain wall¹².

In conclusion, we have described interlayer charge-transfer effects in 2D semiconductor bilayers. Two means of maintaining the requirement of out-of-plane periodicity in DFT calculations have been discussed. We have shown a substantial effect in the TMDs, demonstrating a general formula for finding the size of the effect for general lattice configurations beyond commonly-found high-symmetry stacking orders, showing how charge transfer will be important in understanding the behavior of twistronic 2D material structures, where the local atomic registry and interlayer distance will vary continuously. In contrast to the TMDs, the size of the charge transfer effect in the hexagonal post-transition metal chalcogenide InSe is found to be negligibly small.

Methods

The DFT calculations for $\hbox {WSe}_2$ and InSe in this work were carried out using the plane-wave based VASP code²¹ using the projector augmented wave (PAW) pseudopotentials as distributed with VASP 5.4.4^22,23. We approximated the exchange correlation functional using the generalised gradient approximation (GGA) of Perdew, Burke and Ernzerhof (PBE method)²⁴, while for the local density approximation (LDA) comparisons in Figs. 1 and 2 we use the exchange correlation functional of Ceperley and Alder¹³ as parametrized by Perdew and Zunger¹⁴. In Table 2 we give cutoff radii and valence electron configurations for the VASP pseudopotentials used. The cutoff energy for the plane-waves is set to 600 eV with the in-plane Brillouin zone sampled by a $12 \times 12$ grid. Monolayer crystal structure parameters and interlayer distances are taken from experimental references for bulk crystals^15,25,26.

We also compare results for four of the TMDs using calculations carried out using the Quantum Espresso (QE) package^27,28. A plane-wave cutoff energy of 1090 eV was used for all QE calculations, where the integration over the Brillouin zone was performed using scheme proposed by Monkhorst–Pack with a grid of $12 \times 12\times 1$. All calculations used full relativistic norm-conserving pseudopotentials with spin-orbit interaction included. The exchange correlation functional was approximated using the PBE method.