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

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, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {WSe}_2$$\end{document}WSe2—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.

sion, with a large vacuum between them to ensure the layers are isolated from each other. Here, we compare results for polar 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 7 , (2) a single bilayer supercell, but with a dipole correction applied 11 , and (3) a supercell with two mirror-reflected images of the bilayer 3 .
We consider these methods in XM ′ -stacked WSe 2 (the prime symbol indicating that of the vertical metal-chalcogen pair in the bilayer, the metal atom is in the top layer 12 ), 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 8 , 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.
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 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,  Table 3 www.nature.com/scientificreports/ 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.
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 13 as parametrized by Perdew and Zunger 14 , 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 bilayermonolayer 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 Ŵ -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.

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 MoSe 2 bilayers 7 , 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 Ŵ-point valence band wavefunction has an almost zero out-of-plane www.nature.com/scientificreports/ dipole moment, due to the strong interlayer hybridisation of the Ŵ-point states. The Q-point, which in some cases forms the conduction band minimum, is an intermediate case.
The variation of the potential drop across a P-MX 2 bilayer with in-plane offset, r 0 ( 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 8 , where values of parameters a , q and d 0 , G 1,2,3 are the shortest reciprocal vectors of a monolayer related by 120 • -rotation around z. This formula is applicable to all TMDs with a honeycomb lattice structure, with parameters for MoS 2 , MoSe 2 , WS 2 and WSe 2 calculated using the QE code, shown in Table 4. Since in twisted bilayers the interlayer distance d and in-plane offset 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.

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 γ 15 and ε 16 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 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 ∼2 meV between the vacuum potentials across the bilayer, with the consequence that differences between supercell and correction methods are negligible.

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 7 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, ±δρ , on the top/bottom layers.
(1) � P (r 0 , d) = � a e −q(d−d 0 ) j=1,2,3 sin(G j · r 0 ), 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�ψ|z|ψ� , for valence (VB) and conduction (CB) states at important points in the Brillouin zone ( Q = K/2).  (1), calculated using the QE code. The value in parantheses for WSe 2 is that found above using the VASP code, with only ∼10% variation between different codes in the prediction for P . www.nature.com/scientificreports/ To demonstrate the latter, in Fig. 4 we show the z-coordinate dependence of the difference between the planeaveraged charge density of XM ′ stacked bilayer WSe 2 and that of two isolated 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  www.nature.com/scientificreports/ 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 where ρ(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 WSe 2 layers. This gives δρ = 1.9 × 10 12 cm −2 . We can also roughly estimate the magnitude of charge density transferred between the layers based on the potential drop, P [Eq.
(1)], as ε 0 � P e 2 d XX where d XX is the distance between the inner chalcogen atomic planes (3.14 Å for WSe 2 ). This gives ∼ 10 12 cm −2 for MX ′ /XM ′ stacked bilayer 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 θ P > 2.5 •12 ) the relative in-plane shift of the layers is r 0 = θẑ × r . In a reconstructed twistronic bilayer of a marginally twisted ( θ P ≪ 1 • ) parallel(P)-stacked TMD 12,17 , a pattern of triangular domains forms, with alternating MX ′ /XM ′ stacking, with r 0 = θẑ × r + u(r) , where u(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 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 18 (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 � P (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� P (MX ′ ) . The corresponding values of 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 8 .
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 WSe 2 and InSe in this work were carried out using the plane-wave based VASP code 21 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) 24 , while for the local density approximation (LDA) comparisons in Figs. 1 and 2 we use the exchange correlation functional of Ceperley and Alder 13 as parametrized by Perdew and Zunger 14 .
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 × 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 × 12 × 1 . All calculations used full relativistic norm-conserving pseudopotentials with spin-orbit interaction included. The exchange correlation functional was approximated using the PBE method. with thickness much larger than moiré period ∝ ℓ . 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 WSe 2 bilayer with twist angle θ = 0.6 • , from sum of piezoelectric potential 12 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 ≈ 8 nm corresponding to the width of the domain wall 12 . www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.