# Piezoelectric networks and ferroelectric domains in twistronic superlattices
in WS2/MoS2 and WSe2/MoSe2 bilayers

V. V. Enaldiev1,2,3∗ F. Ferreira1,2, S. J. Magorrian1,2
Vladimir I. Fal’ko,1,2,4∗

1School of Physics and Astronomy, University of Manchester, Oxford Road,
Manchester, M13 9PL, UK
2National Graphene Institute, University of Manchester, Oxford Road,
Manchester, M13 9PL, UK
3Kotel’nikov Institute of Radio-engineering and Electronics,
Russian Academy of Sciences, 11-7 Mokhovaya St, Moscow, 125009 Russia
4Henry Royce Institute for Advanced Materials, University of Manchester,
Oxford Road, Manchester, M13 9PL, UK

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
Twistronic van der Waals heterostrutures offer exciting opportunities for
engineering optoelectronic properties of nanomaterials, in particular, due to
the formation of moiré superlattice structures. In twisted bilayers of
transition metal dichalcogenides moiré superlattice effects are additionally
enriched by the lack of inversion symmetry in each monolayer unit cell. Here,
we use multiscale modeling to establish a rich variety of confinement
conditions for electrons, holes and layer-indirect excitons in twistronic
WX2/MoX2 bilayers (X=S,Se). Such trapping of charge carriers and excitons is
caused by ferroelectric (interlayer) polarisation and piezoelectric effects
generated by the reconstruction of twistronic bilayers into preferential
stacking domains separated by domain wall networks. For almost aligned
bilayers with anti-parallel (AP) orientation of WX2 and MoX2 unit cells, we
find that upon lattice relaxation piezoelectric potential modulation traps
holes and electrons in the opposite corners – WMo and XX (tungsten over
molybdenum versus overlaying chalcogens) – of hexagonal-shaped 2H
(simultaneously WX and XMo) stacking domains, swapping their positions at a
twist angle $\gtrsim 0.2^{\circ}$. This crossover happens at such small angles
(set by a very small lattice mismatch between WX2 and MoX2) that would impose
an alignment accuracy and homogeneity better than $0.1^{\circ}$ for achieving
reproducibility of electronic characteristics of such heterostructures. At the
same time, for all angles, XX corners provide $30$ meV deep traps for the
interlayer excitons. In bilayers with parallel (P) orientation of WX2 and MoX2
unit cells, band edges for both electrons and holes appear in triangular
domains, where WX2 chalcogens set over MoX2 molybdenums. We find that, due to
a weak ferroelectric polarisation, these triangular domains act as $130$ meV
deep quantum boxes for interlayer excitons for twist angles $\lesssim
1^{\circ}$, shifting towards XX stacking sites of the domain wall network at
larger twist angles.

Moiré superlattices are generic for van der Waals structures assembled from
two-dimensional materials (2DM): graphene [1, 2], hexagonal boron nitride [3,
4, 5], and monolayers of transition metal dichalcogenides (TMDs) [6, 7, 8, 9,
10]. For pairs of 2DM with a small lattice mismatch and “twistronic” bilayers
with a small misalignment angle, $\theta$, between the crystallographic axes
moiré superlattices have long periods, $\ell$, giving rise to minibands for
charge carriers [1, 2, 3] and excitons [6, 7, 8, 9, 10]. For twistronic
structures with the longest periods, local lattice reconstruction occurs into
an array of the preferential stacking domains corresponding to the equilibrium
structure of bulk van der Waals crystals [11, 12, 13, 14]. While lattice
reconstruction is stronger in marginally twisted homobilayers of 2D materials
(graphene on graphene, or MX2/MX2 [14, 12, 15]), the very close lattice
constants of TMDs with different metals but the same chalcogen (WX2 and MoX2)
also promote the formation of preferential stacking domains in marginally
twisted heterobilayers [12, 11]. Upon lattice reconstruction, moiré pattern in
such structures takes the form of large-area ($\propto\ell^{2}$) domains
separated by a network of dislocation-like domain walls [12, 15, 11, 16, 17,
18]. Due to the inversion asymmetry of TMD monolayers, the emerging domain
structures differ for bilayers with parallel (P, $\theta\equiv\theta_{P}$) and
anti-parallel (AP, $\theta=\pi+\theta_{AP}$) orientations of their unit cells.
While for P-bilayers the reconstructed moiré pattern constitutes triangular
domains of 3R-type stacking (either WX, or XMo), Fig. 1, domains in AP-
bilayers are hexagonal and feature 2H-type stacking (simultaneously, WX and
XMo) [12, 11, 16, 17].

Here, we demonstrate that such TMD heterostructures exhibit a rich variety of
quantum dot-like features for electrons ($e$), holes ($h$) and interlayer
excitons (iX) which positions within moiré superlattices strongly depend on
the misalignment angle. Because MoX2/WX2 heterostructures have a type-II band
alignment, with the conduction band edge residing in MoX2 and the valence band
edge in WX2, the interlayer hybridization of the band edges is non-resonant
and is superseded by electrostatic potentials due intralayer piezoelectric and
interlayer ferroelectric charge transfers. The resulting variety of features
is illustrated in Fig. 1 for WSe2/MoSe2 heterostructures with both P- and AP-
orientations of the unit cells and various small twist angles. For AP-
heterostructures, strain concentrated at the domain wall network generates
piezocharges with the same signs in the two layers [16], forming distributions
shown in Fig. 1A. This leads to an electrostatic modulation of the electron
(in MoX2) and hole (in WX2) band edges ($\varphi^{e}_{\rm Mo}$ and
$\varphi^{h}_{\rm W}$), with quantum dots forming in XX corners for electrons
and in WMo corners for holes for $|\theta_{AP}^{\circ}|<0.1^{\circ}$. This
arrangement swaps to electron/hole band edges at WMo/XX domain wall vortices
for $0.3^{\circ}<|\theta_{AP}^{\circ}|\lesssim 2^{\circ}$. To compare, lattice
reconstruction in heterostructures with the same orientation of monolayer unit
cells (P-bilayers) produces piezoelectric charges with opposite signs in WX2
and MoX2 monolayers [16], leading to an interlayer polarisation of domain
walls, which interplays with the dominant effect – the interlayer
ferroelectric charge transfer oscillating across moiré supercell, Fig. 1E. For
twist angles $|\theta_{P}^{\circ}|\lesssim 1^{\circ}$, these XMo domains act
as quantum boxes for electrons, holes and interlayer excitons, Fig. 1D. For
$|\theta_{P}^{\circ}|\gtrsim 1^{\circ}$, the piezoelectric interlayer
potential (extended into smaller domains from domain walls) shifts the
electron/hole band edges towards XX stacking areas, followed by a gradual
decrease of the band edge modulation across the moiré supercell and to the iX
delocalisation at larger twist angles.

Figure 1:  Intra- and interlayer charge transfers in twisted WSe2/MoSe2
bilayers. A Maps of piezocharge $\rho$ and B out-of-plane electric
polarisation $P_{z}$ densities in AP-WSe2/MoSe2 bilayers with labeled twist
angles. C,D. Twist angle dependences of band edge energies of electrons and
holes C in XX, WMo, 2H regions for AP- and D in XX, XMo, WX regions for
P-WSe2/MoSe2 bilayers. E. Maps of out-of-plane polarisation, $P_{z}$, for
P-WSe2/MoSe2 bilayers for the same twist angles as for AP-bilayers. Panels on
B and E have the same color bar scale. WS2/MoS2 bilayers show a qualitatively
similar behaviour. For all maps horizontal and vertical directions correspond
to zigzag and armchair crystallographic axes of constituent monolayers,
respectively.

The theory proposed in this paper is based on the multiscale analysis [19, 16]
of atomic reconstruction of twisted bilayers using elasticity theory combined
with density functional theory (DFT) modeling of interlayer adhesion energy
and ferroelectric charge transfer (see Methods). In Fig. 2A,B, we gather the
earlier published DFT results on adhesion energy density [16],
$W(\bm{r}_{0},d)$, for WSe2/MoSe2 and WS2/MoS2 heterostructures. Those are
presented for both P- and AP-orientations (circles and triangles,
respectively) and various stacking configurations characterised by lateral
offsets, $\bm{r}_{0}$, of WX2 and MoX2 honeycomb lattices ($\bm{r}_{0}=0$ for
XX stacking) and interlayer distances, $d$. These DFT data are also encoded
into interpolation formulae for $W(\bm{r}_{0},d)$, as described in Methods,
following the recipes of Refs. [16] and [12]. For P-oriented structures the
most energetically favourable are local stackings with tungstens in WX2
overlaying chalcogens in MoX2 (WX), and molybdenums in MoX2 overlaying
chalcogens in WX2 (XMo). Using analogy with the classification of aligned
phases of bulk TMDs we will refer to WX and XMo domains as 3R-type domains.
For AP-orientation the most favourable is stacking with tungstens in WX2
vertically aligned with chalcogens in MoX2 and molybdenums in MoX2 with
chalcogens in WX2 (similar to 2H-stacking in bulk TMDs). Their energies should
be compared to the configuration-averaged adhesion energy, $f(d)$, which has a
minimum at the interlayer distance $d_{0}$ (configuration 6). The adhesion
energy gained due to forming 2H-type and 3R-type domains, as compared to the
minimum $f(d_{0})$, is about $\approx 0.5$ eV/nm2, which is enough to overcome
the costs of the uniform strain, $\left(\lambda_{\rm W}+\lambda_{\rm
Mo}\right)\delta^{2}\approx 10$ meV/nm2, needed for adjusting the lattice
constants of perfectly aligned WX2 and MoX2 monolayers ($\delta=0.2\%$ for
WS2/MoS2 and $\delta=0.3\%$ for WSe2/MoSe2). This is why moiré superlattice of
slightly twisted WX2/MoX2 bilayers transforms into 2H-type domains for AP
stacking, and into 3R-type domains for P stacking. These domains persist up to
angles $\theta_{P}^{\circ}\sim 2.5^{\circ}$ and $\theta_{AP}^{\circ}\sim
1^{\circ}$, whereas for larger twist angles their areas become too small to
compensate for the energy cost of the dislocation-like domain wall network.

Figure 2:  DFT results for adhesion energy, ferroelectric polarisation, and
their analysis. A and B. Triangles/circles show DFT-computed $d$-dependences
of adhesion energies of AP/P-bilayers for six stacking configurations of
WX2/MoX2 shown in top insets; solid and dashed lines represent fit for AP and
P bilayers; C. DFT-computed values (crosses) and their fit (solid lines) with
equation (1) for electron potential energy shift between MoX2 and WX2
monolayers in P-WX2/MoX2 bilayers for three representative interlayer
distances. Insets show the corresponding $d$-dependence of $\Delta_{0,a,s}$
and their linear fit. D. The same as in C, but for AP-WX2/MoX2 bilayers.
Middle inset: DFT-computed difference between unit cell-averaged potential
energy of a WSe2/MoSe2 heterobilayer (with 3R stacking and distance
$d=6.477$Å), and that of isolated monolayers. The values of $\Delta_{i}$ in C
and D shows that $\Delta_{0}^{AP}=\Delta_{0}^{P}$ with a high accuracy and
$\Delta_{s}\ll\Delta_{0}$.

The interlayer hybridisation of valence band states in one layer with
conduction band states in the other (across the entire band structure) gives
rise to the interlayer charge transfer, which can be considered as a weak
ferroelectric effect. Such a possibility was discussed in relation to
homobilayers of TMDs [20] and was recently observed and analysed in twisted
hBN flakes [21, 22, 23]. The formation of a charge double layer at MoX2/WX2
interface produces potential energy jump $\Delta$, quantified in Figs. 2C and
2D for P- and AP-bilayers, respectively. To analyze such a charge transfer by
means of DFT modeling we constructed a supercell containing two mirror-
reflected WX2/MoX2 heterobilayers placed 30 Å apart (to scale down interaction
between them) and computed $\Delta$ for P- and AP-WX2/MoX2 heterostructures in
the same stacking configurations (same $\bm{r}_{0}$) as in Fig. 2A,B across
the range of interlayer distances that covers the minima of $W(\bm{r}_{0},d)$
for these configurations. All these DFT data are accurately described by an
interpolation formula,

$\Delta(\bm{r}_{0},d)=A_{0}\left[1-q_{0}z\right]+A_{a}^{P/AP}\left[1-q_{a}z\right]\sum_{n=1,2,3}\sin\left(\bm{G}_{n}\bm{r}_{0}\right).$
(1)

Here, $z=d-d_{0}$ is a deviation of the interlayer distance from the minimum
of $f$ (configuration average of $W$). The first term in equation (1) is
specific for heterobilayers ($A_{0}=0$ for homobilayers), and its influence on
the indirect (interlayer) band gap comes from the variation of the interlayer
distance across the moiré pattern. The second term is generic for both homo-
and heterobilayers in P-orientation and broken $z\to-z$ symmetry
($\bm{r}_{0}\neq 0$), but $A_{a}^{AP}=0$ (for homobilayers, this is due to the
exact inversion symmetry of their lattice). The values of parameters in this
interpolation formula are listed in Table 1. This interpolation formula
emerged from fitting the DFT data using a more general periodic function,
$\Delta=\Delta_{0}+\Delta_{a}\sum\sin(\bm{G}_{n}\bm{r}_{0})+\Delta_{s}\sum\cos(\bm{G}_{n}\bm{r}_{0})$,
involving the first star ($\pm\bm{G}_{1,2,3}$) of reciprocal lattice vectors
of the TMD monolayer, $\Delta_{i}=A_{i}e^{-q_{i}(d-d_{0})}\approx
A_{i}\left[1-q_{i}(d-d_{0})\right]$. We find that
$\Delta_{s},\Delta^{AP}_{a}\ll\Delta_{0}$ (see insets in Fig. 2C,D).
Furthermore, the DFT analysis of the atomic-scale potential variations, shown
in the inset in Fig. 2D, indicates that ferroelectric charge transfer happens
between the two inner chalcogen sublayers, which we use to compute interlayer
polarisation maps in Fig. 1.

Table 1: Parameters used in the interpolation formulae for ferroelectric charge transfer (equation (1)); for adhesion energy density parameters of WX2/MoX2 bilayers corresponding to DFT data in Fig. 2 (see Methods); elastic moduli [24, 25] of the WX2 and MoX2 monolayers used in lattice relaxation and their piezocoefficients [26, 27]. X |  $A_{0}$, |  $q_{0}$, |  $A_{a}$, | $q_{a}$
---|---|---|---|---
| meV |  Å-1 | meV |  Å-1
S | 29.9 | 1.73 | 12.8 | 2.53
Se | 26.5 | 1.89 | 9.53 | 3.12
X |  $w_{1}$, |  $w_{2}$, |  $Q$, |  $d_{0}$,
|  eV/nm2 |  eV/nm2 |  Å-1 | Å
S | 0.175 | 0.021 | 3.07 | 6.5
Se | 0.128 | 0.02 | 2.93 | 6.9
X |  $\varepsilon$, |  $\lambda_{\rm W}$/$\lambda_{\rm Mo}$, |  $\mu_{\rm W}$/$\mu_{\rm Mo}$, |  $\left|e_{11}^{\rm W}\right|/\left|e_{11}^{\rm Mo}\right|$,
|  eV/nm4 | N/m | N/m |  $10^{-10}$C/m
S | 214 | 52.5/83.2 | 72.5/70.9 | 2.74/2.9
Se | 189 | 29.7/42.3 | 48.4/49.6 | 2.03/2.14

Lacking inversion symmetry, MoX2 and WX2 monolayers are piezoelectric
crystals, so that inhomogeneous strain in each of them produces piezocharges
with density,
$\rho=e_{11}\left[2\partial_{x}u_{xy}+\partial_{y}(u_{xx}-u_{yy})\right]$. It
is important to note that piezoelectric charge transfer leads to the most
pronounced effects specifically in bilayers with AP-orientation of the unit
cells. This is because atomic reconstruction towards 2H-type stacking domains
brings together the lattices of the two layers, so that their strain tensors
are $u^{\rm W}_{ij}\approx-u_{ij}^{\rm Mo}$, whereas for AP-bilayers the signs
of WX2 and MoX2 piezocoefficients are inverted, $e_{11}^{\rm W}\approx-
e_{11}^{\rm Mo}$. As a result, piezocharges in the two layers have the same
sign and add up, creating a potential varying across moiré supercell
simultaneously on the top and bottom layers. The cumulative effect of strain
and ferroelectric charge transfer on the electron and hole band edges in the
AP-bilayers is illustrated in Fig. 3 using WSe2/MoSe2 heterostructures, as an
example. The results shown in Fig. 3A were obtained neglecting interlayer
hybridisation of (conduction-conduction and valence-valence) band edge states,
which are non-resonant for WX2/MoX2 bilayers (in contrast to homobilayers and
some MX2/M′X${}^{\prime}_{2}$ heterostructures with almost resonant band
edges, where it plays an important role [9, 28]). For very small twist angles,
$|\theta_{AP}^{\circ}|<0.1^{\circ}$, electron/hole band edges appear at XX/WMo
vortices of the domain wall network, which creates strongly triangulated
quantum dot confinement for the corresponding charge carriers, see insets in
Fig. 3A. The electron/hole states localised on a periodic array of quantum
dots would produce a narrow band (creating conditions for strongly correlated
phases in n- and p-doped bilayers, like those observed in graphene
superlattices [1, 2]). For $0.3^{\circ}<\theta_{AP}^{\circ}\lesssim
2^{\circ}$, the conduction and valence band edges swap position, with quantum
dots for electrons/holes appearing in WMo/XX stacking areas, with the computed
confinement profiles shown in Fig. 3B insets. Emergence of deep trapping
potentials for charge carriers in AP-WSe2/MoSe2 bilayers induced by lattice
reconstruction has been recently demonstrated in Ref. [29]. The quantitatively
smaller band edge variations caused by the ferroelectric effect additionally
modulate the layer-indirect band gap in these type-II heterostructures. This
is taken into account in the maps shown in Fig. 3, though it does not change
the positions of the conduction band minima and valence band maxima. For
larger misalignment angles $\theta_{AP}^{\circ}\gtrsim 2^{\circ}$, where
lattice reconstruction is weak [16, 12], the piezo- and ferroelectric
contributions become comparable, producing a harmonic variation of electron’s
(in MoX2) and hole’s (in WX2) band edges,

$\displaystyle\varphi^{e}_{\rm
Mo}=-U_{0}\sum_{n=1,2,3}\sin\left(\bm{g}_{n}\bm{r}\right),\quad\bm{g}_{n}=\delta\cdot\bm{G}_{n}-\theta\hat{z}\times\bm{G}_{n};$
$\displaystyle\varphi^{h}_{\rm
W}=U_{0}\sum_{n=1,2,3}\left[\sin\left(\bm{g}_{n}\bm{r}\right)-\frac{1}{2}\cos\left(\bm{g}_{n}\bm{r}\right)\right],\qquad\qquad\quad$
(2)

with $U_{0}\approx 12$ meV for diselenides and $U_{0}\approx 15$ meV for
disulfides.

We note that the above-mentioned lateral piezocharge transfer does not affect
the energy of the bound electron-hole pairs such as interlayer excitons (iX)
[30], whereas the indirect band gap modulation by vertical (MoX2 $\to$ WX2)
charge transfer does. In twistronic AP-bilayers, this sets the iX energy
minimum in the XX corners of the moiré supercell, see Fig. 3C, typically,
$\sim 20$ meV lower than in 2H and WMo areas for WSe2/MoSe2, and $\sim 30$ meV
for WS2/MoS2. At low temperatures, these XX corners of the domain wall network
would act as traps for both bright and dark interlayer excitons, potentially
providing conditions for accumulation of dark iX (traced to the spin-valley
splitting [31, 32]) due to their lower energy, or even create conditions for
the single-photon emission.

Figure 3:  Energy profiles for electrons, holes and interlayer excitons in AP-
bilayers. A,B. Maps of electron (top) and hole (bottom) band edge energies for
twisted AP-WSe2/MoSe2 bilayers with labeled twist angles. Insets illustrate
shape of quantum dot potentials at marginal twist angles. Horizontal and
vertical directions correspond to zigzag and armchair crystallographic
directions of constituent monolayers, respectively. C. Twist-angle-dependence
of interlayer exciton energies in XX, WMo, and 2H regions of AP-WSe2/MoSe2 and
AP-WS2/MoS2 bilayers.

In contrast to the AP orientation, in P-bilayers the contributions of
piezoelectric and ferroelectric charge transfers produce similar effects. This
is because the piezocoefficients in monolayers with the same orientation of
unit cells have the same sign, $e_{11}^{\rm W}\approx e_{11}^{\rm Mo}$, and as
$\bm{u}^{\rm W}\approx-\bm{u}^{\rm Mo}$, leading to opposite signs of the
piezocharges in top/bottom layers [16] and piezopotentials $\phi_{\rm
Mo}\approx-\phi_{\rm W}$. Consequently, the combined effect of piezo- and
ferroelectric contributions (i.e. charge double layer) electrostatically
modulates the layer-indirect band gap in the type-II heterostructure WX2/MoX2.
and the iX energy, $U_{{\it iX}}$. To mention, for the AP-heterostructures the
lowest spin-valley states of iX are bright, according to the optical
transition selection rules [31, 32, 33]. The resulting energy landscape for
the bright ground state of iX is shown in Fig. 4A, indicating that the XMo
regions form boxes (with almost the same depth $\approx 130$ meV for sulphides
and selenides) accumulating iXs in heterostructures with
$\theta_{P}^{\circ}\lesssim 1^{\circ}$. As for P-bilayers we neglected
additions to the iX energy from the off-resonant interlayer hybridisation of
(conduction-conduction and valence-valence) band edge states (see Methods).
The localisation of layer-indirect excitons in XMo domains is determined by
the dominant role of the ferroelectric interlayer polarisation in such
marginally twisted structures, whereas for twist angles
$\theta_{P}^{\circ}\gtrsim 1.2^{\circ}$, the piezoelectric polarisation of
domain walls between 3R stacking domains pushes the iX energy minima towards
XX-stacking sites of triangular domain wall network, Fig. 4B. For bilayers
with even larger angles ($\theta_{P}^{\circ}\gtrsim 3^{\circ}$), where the
lattice reconstruction is weak, the iX energy difference between XMo and XX
areas progressively disappears, Fig. 4C, opening channels for iX propagation
across the supercell. In the latter regime we expect the formation of
dispersive minibands for iXs, as discussed in Ref. [28]. Similarly, the
crossover of iX confinement from XMo domains to XX stacking areas at
$1^{\circ}<|\theta_{P}^{\circ}|\lesssim 1.2^{\circ}$ may also take a form of
delocalised dispersive excitons.

Overall the proposed theory establishes a rich variety of confinement
conditions for electrons, holes and layer-indirect excitons across moiré
superlattice in twistronic type-II heterostructures WX2/MoX2. The richness of
predicted behaviour stems from the interplay between piezoelectric and weak
ferroelectric effects. In contrast to homobilayers [15], the interplay between
lattice reconstruction and charge transfer in slightly twisted same-chalcogen
heterobilayers makes localisation of individual charge carriers and optically
active excitations in AP-bilayers sensitive to a small variation of the twist
angle, in the range of $\theta_{AP}^{\circ}\sim 0.2^{\circ}$. While such a
sensitivity may represent a long-term opportunity for tuning the electronic
transport properties and photoemission dynamics of such systems, in the first
instance, it suggests that a homogeneity of the twist angle with accuracy
better than $0.1^{\circ}$ will be necessary for achieving a systematic
reproducibility in the optoelectronic characteristics of WX2/MoX2
heterostructures.

Figure 4:  Energy landscapes for interlayer excitons (iX) in P-WSe2/MoSe2
(top) and P-WS2/MoS2 (bottom). A. For $\theta_{P}^{\circ}=0^{\circ}$ iX are
confined inside triangular XMo domains. B. For
$\theta_{P}^{\circ}=1.5^{\circ}$ minima of iX energy appear in XX regions. On
maps A and B horizontal and vertical directions correspond to zigzag and
armchair crystallographic directions of constituent monolayers, respectively.
C. Twist-angle dependence of iX energy in XX, XMo, and WX regions of moiré
supercell, with iX energy minima shifting from XMo domains to XX vortices at
$|\theta_{P}^{\circ}|\approx 1^{\circ}$.

## Methods

Multiscale modeling [19, 16] in the reported study exploits the fact that the
moiré superlattice period $\ell\approx a/\sqrt{\theta_{P,AP}^{2}+\delta^{2}}$
is orders of magnitude longer (for $|\theta_{P,AP}^{\circ}|\leq 6^{\circ}$
which corresponds to $|\theta_{P,AP}|<0.1$) than the lattice constant, $a$, of
each monolayer. Because of such length scale separation local properties of
each small area of a moiré supercell are quantified by performing ab initio
DFT computations for ideally aligned pair of monolayers, however, with the
locally defined lateral offset $\bm{r}_{0}$ and interlayer distance, $d$,
which vary across the supercell. For convenience we set $\bm{r}_{0}=0$ for XX
stacking (coinciding lateral positions of all chalcogens) for both P and AP
orientations of the crystals. For the DFT computations we considered
commensurate monolayers and used van der Waals corrected code (optB-88vdW)
neglecting spin-orbit coupling and setting a plane-wave cutoff of 816.34 eV
(60 Ry). The separation of micro- and mesoscopic scales is accounted for by
$\bm{r}_{0}(\bm{r})=\delta\bm{r}+\theta_{P,AP}\hat{z}\times\bm{r}+\bm{u}^{\rm
W}-\bm{u}^{\rm Mo}$, where the first term originate from a small difference
between MoX2 and WX2 lattice constants, $\delta=2\left(a_{\rm Mo}-a_{\rm
W}\right)/\left(a_{\rm Mo}+a_{\rm W}\right)$, the second term takes into
account a fixed twist angle, and the last comes from lateral deformations in
each layer bringing their lattices to the globally energetically favourable
state. Displacements $\bm{u}^{\rm W/Mo}$ are related to strain in each layer
$2u_{ij}^{\rm W/Mo}=\partial_{i}u^{\rm W/Mo}_{j}+\partial_{j}u^{\rm
W/Mo}_{i}$, and local interlayer distance, $d$, varies across the supercell
adjusting to the optimal value for each local stacking configuration.

We find the energetically favourable global state of the superlattice by
minimising the total energy of the deformed crystal which includes the
interlayer adhesion energy, $W_{P/AP}(\bm{r}_{0},d)$, and the energy cost of
lateral deformations. To cross from the microscopic scales (at which the
adhesion energy is computed using DFT) to the mesoscale, we implement an
interpolative analytical description for $W_{P/AP}(\bm{r}_{0},d)$ and use it
to formulate an effective energy density functional [16]. Numerical
minimisation of $\mathcal{E}$ on a fine grid produces a detailed structure of
the strain field and interatomic distances across the moiré supercell. This
information is, then, used to compute piezoelectric charges, generated by
inhomogeneity of strain and the interlayer charge transfer due to a weak
ferroelectric effect. This last step requires additional DFT calculations,
which are performed for the same set of stacking configurations, as in the
analysis of adhesion energy and used to build an analytical interpolation
formula for the potential step $\Delta(\bm{r}_{0},d)$ caused by the interlayer
charge transfer. Finally, the computed $\bm{r}_{0}(\bm{r})$,
$d(\bm{r}_{0}(\bm{r}))$, on-layer piezopotentials $\phi_{\rm W}$ in WX2 and
$\phi_{\rm Mo}$ in MoX2, and $\Delta$ are combined together to compute maps
displayed in Figs. 1,3, and 4.

Interpolation formula for adhesion energy $W_{P/AP}$. First of all, we
implement DFT to compute adhesion energies of two aligned and commensurate WX2
and MoX2 monolayers for several lateral offsets corresponding to stacking
arrangements, shown in the insets in Fig. 2, and various interlayer distances.
The commensurability of the monolayer lattices is attained by in-plane
stretching/shrinking of one of the monolayers to the lattice of the other,
keeping X-X intersublayer distance in each layer fixed to its monolayer value.
Due to a small lattice mismatch between WX2 and MoX2 ($\delta=0.3\%$ for X=Se
and $0.2\%$ for X=S), the DFT results for adhesion energy (Fig. 2A,B) are
independent of which way such commensurability is imposed [16]. An analytical
description of the interlayer adhesion energy as a function of stacking is
achieved by using the Fourier expansion of $W_{P/AP}(\bm{r}_{0},d)$ for a
fixed $d$, keeping track of the constant term and six harmonics from the first
star of reciprocal lattice vectors, $\pm\bm{G}_{1,2,3}$ ($|\bm{G}_{n}|=G$) of
the monolayer. The resulting interpolational formula reads,

$\begin{split}W_{P/AP}(\bm{r}_{0},d)\approx f(d_{0})+\varepsilon
z^{2}+w_{1}\left[1-Qz\right]\sum_{n=1,2,3}\cos\left(\bm{G}_{n}\bm{r}_{0}\right)\\\
+w_{2}\left[1-Gz\right]\sum_{n=1,2,3}\sin\left(\bm{G}_{n}\bm{r}_{0}+\gamma_{P/AP}\right).\end{split}$

Here, $z=d-d_{0}$ is a deviation of the interlayer distance from the minimum
of $f$ (configuration-averaged $W(d)$ coinciding with $W(d)$ for the 6th
configuration in Fig. 2A,B), so that $f\approx f(d_{0})+\varepsilon z^{2}$ is
an expansion of $f$ around its minimum, $d_{0}$; $\gamma_{AP}=0$ and
$\gamma_{P}=\pi/2$, and $\varepsilon$, $w_{1,2}$, $Q$ are fitting parameters
listed in Table 1. This offers a simplified version of the interpolation
formula used in Ref. [16]. To account for out-of-plane relaxation of WX2 and
MoX2 we minimise $W_{P/AP}$ with respect to $z$, which enables us to find a
local relation $z=Z(\bm{r}_{0})$. Note that the interlayer distance variation,
$Z(\bm{r}_{0})$, from the lowest to the highest energy stackings, $\sim
d_{0}\pm 0.3$Å, is much smaller than interlayer distance itself, which
justifies the expansion implemented in $W_{P/AP}$. Also, as it has been shown
earlier [16], the energy cost of bending deformations of WX2 and MoX2
monolayers is negligibly small (compared with that of lateral strain and
variations of adhesion energy) due to the long range of the deformations of
these atomically thin membranes (the same approximation was used in Ref.[19]
when analysing graphene-hBN heterostructures). Overall this determines the
locally adjusted energy, dependent only on the local interlayer offset
$\bm{r}_{0}$ as

$\displaystyle\widetilde{W}_{P/AP}(\bm{r}_{0})=-\varepsilon
Z^{2}(\bm{r}_{0})+\sum_{n=1}^{3}\left[w_{1}\cos\left(\bm{G}_{n}\bm{r}_{0}\right)+w_{2}\sin\left(\bm{G}_{n}\bm{r}_{0}+\gamma_{P/AP}\right)\right],$
$\displaystyle
Z(\bm{r}_{0})=\frac{1}{2\varepsilon}\sum_{n=1}^{3}\left[w_{1}Q\cos\left(\bm{G}_{n}\bm{r}_{0}\right)+w_{2}G\sin\left(\bm{G}_{n}\bm{r}_{0}+\gamma_{P/AP}\right)\right].$

Mesoscale lattice reconstruction in twisted bilayers. Having defined a locally
adjusted adhesion energy function, $\widetilde{W}_{P/AP}$, we minimise the
total energy,

$\mathcal{E}=\int_{\rm supercell}d^{2}\boldsymbol{r}\left\\{\sum_{l={\rm
W,Mo}}\left[\frac{\lambda_{l}}{2}\left(u_{ii}^{(l)}\right)^{2}+\mu_{l}u_{ij}^{(l)}u_{ji}^{(l)}\right]+\widetilde{W}_{AP/P}\left(\bm{r}_{0}(\bm{r})\right)\right\\},$

with respect to displacement fields $\bm{u}^{\rm W}$ and $\bm{u}^{\rm Mo}$
implicit in
$\bm{r}_{0}(\bm{r})=\delta\bm{r}+\theta\hat{z}\times\bm{r}+\bm{u}^{\rm
W}-\bm{u}^{\rm Mo}$. Here $\lambda$, $\mu$ are elastic moduli, which values
are listed in Table 1. To compute the optimal distributions of the
displacement fields $\bm{u}^{\rm W}(\bm{r})$ and $\bm{u}^{\rm Mo}(\bm{r})$,
which we keep periodic, with the moiré superlattice period, we solve a system
of Lagrange-Euler equations using the finite difference method. Here we use a
homogeneous grid aross a moiré supercell to solve a system of non-linear
algebraic equations, supplemented by periodic boundary conditions based on the
interior point method, as implemented in the GEKKO Optimization Suite [34].

DFT analysis of interlayer charge transfer. The analysis of weak
ferroelectricity in this work employs the VASP package [35] with projector
augmented wave (PAW) pseudopotentials. We approximated the exchange
correlation functional using the generalised gradient approximation (GGA) of
Perdew, Burke and Ernzerhof [36]. The cutoff energy for the plane-waves was
set to 600 eV and the in-plane Brillouin zone sampled by a $12\times 12$ grid.
The in-plane lattice constants of the constituent monolayers are strained to
the mean of their experimental values, taken from Refs. [37] and [38] for MoX2
and WX2, respectively, while the intralayer chalcogen-chalcogen distances are
left unchanged for each constituent monolayer. In the inset in Fig. 2D we
illustrate the DFT-computed behaviour of the difference between the plane-
averaged electron potential energy for a WSe2/MoSe2 heterobilayer (for 3R
stacking and interlayer distance $d=6.477$Å), and that from the sum of the
potential energies from isolated individual WSe2 and MoSe2 monolayers. To
avoid having to artificially resolve the mismatch in potential energy which
would arise at the boundary of a supercell containing a single heterobilayer,
results for heterobilayers in this work were computed using a supercell
containing two mirror-reflected images of the structure.

It should be noted that the van der Waals corrected DFT (optB-88vdW) used in
determining the interlayer adhesion energies [16] gives a slightly larger
optimal interlayer distance ($6.56$Å for WSe2/MoSe2, and $6.25$Å for WS2/MoS2)
for 2H/3R stacking configuration in WX2/MoX2 bilayers than the values found in
experiment for the corresponding homobilayers ($6.46$Å, $6.48$Å, $6.15$Å and
$6.17$Å for MoSe2, WSe2, MoS2, and WS2, respectively [37, 38]). However, usage
of a linear fit for $\Delta_{0}$ and $\Delta_{a}$ with relative interlayer
distance dependence $d-d_{0}$ allows us to take the difference into account by
counting the variations of $\Delta$ from a shifted interlayer distance
$d-d_{0}\to Z(\bm{r}_{0}(\bm{r}))-0.1$Å.

Piezo- and ferroelectric potentials for electrons and holes induced by piezo-
and ferroelectric charge transfers. Here we combine the potential jump
$\Delta$ at the double layer induced by the weak ferroelectric effect with the
potential created by piezoelectric charges, accumulated around domain walls in
the reconstructed moiré superlattice. In the analysis of the latter we also
take into account that the electric field generated by combined piezocharges
of AP-oriented WX2 and MoX2 layers is screened by the dielectric environment,
in, e.g., hBN-encapsulated bilayers. In particular we take into account the
in-/out-of-plane anisotropy of hBN dielectric parameters ([39]
$\epsilon_{||}=6.93$ and $\epsilon_{\perp}=3.76$). Additionally, for both P
and AP orientations we take into account the in-plane dielectric
polarisability of WX2 and MoX2, $\alpha_{2D}^{\rm W/Mo}=d_{0}(\epsilon^{\rm
WX/MoX}_{||}-1)/4\pi$ [40], leading to the induced charges

$\rho^{\rm W/Mo}_{ind}=\alpha^{\rm W/Mo}_{2D}\bm{\nabla}^{2}\phi_{\rm W/Mo},$

where $\epsilon^{\rm WX}_{||}$ and $\epsilon^{\rm MoX}_{||}$ are in-plane
dielectric permittivities of bulk WX2 and MoX2 crystals,
$\bm{\nabla}=(\partial_{x},\partial_{y})$, and $\phi_{\rm W/Mo}$ is the
electrostatic piezopotential in the W/Mo layer. We note that piezocharge
densities $\rho^{\rm W/Mo}$ are symmetric with respect to the mirror plane in
each layer. Therefore, the piezocharges in one monolayer are screened only by
the out-of-plane polarised charges in the opposite layer. Under this
assumption the magnitudes of out-of-plane polarised charges are

$\sigma_{\rm W/Mo}=\frac{\left(\epsilon_{\perp}^{\rm
WX/MoX}-1\right)}{2\epsilon_{\perp}^{\rm WX/MoX}}\rho^{\rm Mo/W},$

where $\epsilon_{\perp}^{\rm WX}$ ($\epsilon_{\perp}^{\rm MoX}$) is the out-
of-layer dielectric permittivity of bulk WX2 (MoX2) crystal. Here, we used the
following values of dielectric permittivities $\epsilon_{||}^{\rm
WX}=14.6/16$, $\epsilon_{||}^{\rm MoX}=16.3/17.9$ and $\epsilon_{\perp}^{\rm
WX}=6.7/7.9$, $\epsilon_{\perp}^{\rm MoX}=7.3/8.5$ known for tungsten and
molybdenum disulfide/diselenide bulk crystals [39]. In our electrostatic model
planes hosting charges $\sigma_{\rm W/Mo}$ are distant from the mirror plane
in corresponding layer by distance $d_{0}/2$, with density $+\sigma_{\rm
W/Mo}$ from outer and $-\sigma_{\rm W/Mo}$ from inner sides of the bilayer
structure (here we assume $d_{0}$ is interlayer distance). Finally, we find
$\phi_{\rm W/Mo}$ by solving the Poisson equation with appropriate boundary
conditions on the planes with the charges discussed above.

The total potential energy of charge carriers in MoX2 and WX2 layers is a sum
of the piezo- and ferroelectric contributions, where the ferroelectric part is
determined by
$\frac{1}{2}\left[\Delta(\bm{r}_{0}(\bm{r}))-\Delta_{0}(d_{2H/3R})\right]$.
Thus, for potential energy of electrons and holes shown in Figs. 3 we use
$\varphi_{\rm Mo}^{e}=-e\phi_{\rm
Mo}+\frac{1}{2}\left[\Delta-\Delta_{0}(d_{2H/3R})\right]$ and $\varphi_{\rm
W}^{h}=e\phi_{\rm W}+\frac{1}{2}\left[\Delta-\Delta_{0}(d_{2H/3R})\right]$,
respectively. Then iX energy (Fig. 4) is their sum $U_{iX}=-e(\phi_{\rm
Mo}-\phi_{\rm W})+\Delta$.

Finally, we note that in calculation of energy profiles for electrons/holes
(Fig. 3) and iXs (Fig. 4) we neglect the energy modulation resulting from
interlayer hybridisation of conduction-conduction and valence-valence band
edge states. This is because in type-II heterostructures WX2/MoX2 the band
edges are non-resonant, with [41, 42, 43] $\sim 500$ meV offsets in conduction
and valence bands, which leads only to small corrections in the second order
of perturbation theory. Therefore, its influence on the band edges in WX2/MoX2
bilayers should be much smaller than in TMD homobilyers where states are
resonantly coupled, and where characteristic tunneling matrix element is,
anyway, smaller ($\sim 20$ meV [44, 28, 9]) than the effects discussed above.

Aknowledgements. We thank Mingxin Chen, Andre Geim, Roman Gorbachev, Hongkun
Park, Alex Summerfield, David Ruiz-Tijerina, Wang Yao, and Hongyi Yu for
useful discussions. This work has been supported by EPSRC grants EP/S019367/1,
EP/S030719/1, EP/N010345/1, EP/V007033/1; ERC Synergy Grant Hetero2D; Lloyd’s
Register Foundation Nanotechnology grant; European Graphene Flagship Project,
and EU Quantum Technology Flagship project 2D-SIPC. V.V.E. (piezo- and
ferrocharge model) acknowledges the support of the Russian Science Foundation
(project no. 16-12-10411).

## References

  * [1] Cao, Y., Fatemi, V., Fang, S., Watanabe, K., Taniguchi, T., Kaxiras, E., and Jarillo-Herrero, P.  Nature 556(7699), 43 (2018).
  * [2] Cao, Y., Fatemi, V., Demir, A., Fang, S., Tomarken, S. L., Luo, J. Y., Sanchez-Yamagishi, J. D., Watanabe, K., Taniguchi, T., Kaxiras, E., and et al.  Nature 556(7699), 80 (2018).
  * [3] Ponomarenko, L. A., Gorbachev, R. V., Yu, G. L., Elias, D. C., Jalil, R., Patel, A. A., Mishchenko, A., Mayorov, A. S., Woods, C. R., Wallbank, J. R., and et al.  Nature 497(7451), 594–597 May (2013).
  * [4] Dean, C. R., Wang, L., Maher, P., Forsythe, C., Ghahari, F., Gao, Y., Katoch, J., Ishigami, M., Moon, P., Koshino, M., and et al.  Nature 497(7451), 598 (2013).
  * [5] Lee, M., Wallbank, J. R., Gallagher, P., Watanabe, K., Taniguchi, T., Fal’ko, V. I., and Goldhaber-Gordon, D.  Science 353(6307), 1526–1529 (2016).
  * [6] Kunstmann, J., Mooshammer, F., Nagler, P., Chaves, A., Stein, F., Paradiso, N., Plechinger, G., Strunk, C., Schüller, C., Seifert, G., and et al.  Nature Physics 14(8), 801 (2018).
  * [7] Rivera, P., Yu, H., Seyler, K. L., Wilson, N. P., Yao, W., and Xu, X.  Nature Nanotechnology 13(11), 1004–1015 (2018).
  * [8] Nayak, P. K., Horbatenko, Y., Ahn, S., Kim, G., Lee, J.-U., Ma, K. Y., Jang, A.-R., Lim, H., Kim, D., Ryu, S., and et al.  ACS Nano 11(4), 4041–4050 (2017).
  * [9] Alexeev, E. M., Ruiz-Tijerina, D. A., Danovich, M., Hamer, M. J., Terry, D. J., Nayak, P. K., Ahn, S., Pak, S., Lee, J., Sohn, J. I., and et al.  Nature 567(7746), 81 (2019).
  * [10] Tran, K., Moody, G., Wu, F., Lu, X., Choi, J., Kim, K., Rai, A., Sanchez, D. A., Quan, J., Singh, A., and et al.  Nature 567(7746), 71 (2019).
  * [11] Rosenberger, M. R., Chuang, H.-J., Phillips, M., Oleshko, V. P., McCreary, K. M., Sivaram, S. V., Hellberg, C. S., and Jonker, B. T.  ACS Nano 14(4), 4550–4558 (2020).
  * [12] Weston, A., Zou, Y., Enaldiev, V., Summerfield, A., Clark, N., Zólyomi, V., Graham, A., Yelgel, C., Magorrian, S., Zhou, M., Zultak, J., Hopkinson, D., Barinov, A., Bointon, T. H., Kretinin, A., Wilson, N. R., Beton, P. H., Fal’ko, V. I., Haigh, S. J., and Gorbachev, R.  Nature Nanotechnology 15(7), 592–597 Jul (2020).
  * [13] Alden, J. S., Tsen, A. W., Huang, P. Y., Hovden, R., Brown, L., Park, J., Muller, D. A., and McEuen, P. L.  Proceedings of the National Academy of Sciences 110(28), 11256–11260 (2013).
  * [14] Yoo, H., Engelke, R., Carr, S., Fang, S., Zhang, K., Cazeaux, P., Sung, S. H., Hovden, R., Tsen, A. W., Taniguchi, T., and et al.  Nature Materials 18(5), 448 (2019).
  * [15] Edelberg, D., Kumar, H., Shenoy, V., Ochoa, H., and Pasupathy, A. N.  Nature Physics 16, 1097 Jul (2020).
  * [16] Enaldiev, V. V., Zólyomi, V., Yelgel, C., Magorrian, S. J., and Fal’ko, V. I.  Phys. Rev. Lett. 124, 206101 May (2020).
  * [17] Naik, M. H. and Jain, M.  Phys. Rev. Lett. 121, 266401 Dec (2018).
  * [18] Carr, S., Massatt, D., Torrisi, S. B., Cazeaux, P., Luskin, M., and Kaxiras, E.  Physical Review B 98(22), 224102 (2018).
  * [19] Jung, J., DaSilva, A. M., MacDonald, A. H., and Adam, S.  Nature Communications 6(1), 6308 Feb (2015).
  * [20] Li, L. and Wu, M.  ACS Nano 11(6), 6382–6388 Jun (2017).
  * [21] Stern, M. V., Waschitz, Y., Cao, W., Nevo, I., Watanabe, K., Taniguchi, T., Sela, E., Urbakh, M., Hod, O., and Shalom, M. B.  arXiv:2010.05182 (2020).
  * [22] Woods, C., Ares, P., Nevison-Andrews, H., Holwill, M., Fabregas, R., Guinea, F., Geim, A., Novoselov, K., Walet, N., and Fumagalli, L.  arXiv:2010.06914 (2020).
  * [23] Yasuda, K., Wang, X., Watanabe, K., Taniguchi, T., and Jarillo-Herrero, P.  arXiv:2010.06600 (2020).
  * [24] Iguiñiz, N., Frisenda, R., Bratschitsch, R., and Castellanos-Gomez, A.  Advanced Materials 31(10), 1807150 (2019).
  * [25] Androulidakis, C., Zhang, K., Robertson, M., and Tawfick, S.  2D Materials 5(3), 032005 (2018).
  * [26] Zhu, H., Wang, Y., Xiao, J., Liu, M., Xiong, S., Wong, Z. J., Ye, Z., Ye, Y., Yin, X., and Zhang, X.  Nature Nanotechnology 10(2), 151 (2015).
  * [27] Rostami, H., Guinea, F., Polini, M., and Roldán, R.  npj 2D Materials and Applications 2(1), 1–6 (2018).
  * [28] Ruiz-Tijerina, D. A. and Fal’ko, V. I.  Physical Review B 99(12), 125424 (2019).
  * [29] Shabani, S., Halbertal, D., Wu, W., Chen, M., Liu, S., Hone, J., Yao, W., Basov, D. N., Zhu, X., and Pasupathy, A. N.  arXiv:2008.07696 (2020).
  * [30] Szyniszewski, M., Mostaani, E., Drummond, N. D., and Fal’ko, V. I.  Phys. Rev. B 95, 081301 Feb (2017).
  * [31] Kormányos, A., Burkard, G., Gmitra, M., Fabian, J., Zólyomi, V., Drummond, N. D., and Fal’ko, V.  2D Materials 2(2), 022001 apr (2015).
  * [32] Yu, H., Wang, Y., Tong, Q., Xu, X., and Yao, W.  Phys. Rev. Lett. 115, 187002 Oct (2015).
  * [33] Danovich, M., Ruiz-Tijerina, D. A., Hunt, R. J., Szyniszewski, M., Drummond, N. D., and Fal’ko, V. I.  Phys. Rev. B 97, 195452 May (2018).
  * [34] Beal, L., Hill, D., Martin, R., and Hedengren, J.  Processes 6(8), 106 Jul (2018).
  * [35] Kresse, G. and Furthmüller, J.  Physical Review B 54, 11169–11186 Oct (1996).
  * [36] Perdew, J. P., Burke, K., and Ernzerhof, M.  Phys. Rev. Lett. 77, 3865–3868 Oct (1996).
  * [37] Bronsema, K. D., Boer, J. L. D., and Jellinek, F.  Zeitschrift fur anorganische und allgemeine Chemie 540(9-10), 15–17 September (1986).
  * [38] Schutte, W., Boer, J. D., and Jellinek, F.  Journal of Solid State Chemistry 70(2), 207–209 October (1987).
  * [39] Laturia, A., Van de Put, M. L., and Vandenberghe, W. G.  npj 2D Materials and Applications 2(1), 1–7 (2018).
  * [40] Cudazzo, P., Tokatly, I. V., and Rubio, A.  Phys. Rev. B 84, 085406 Aug (2011).
  * [41] Gong, C., Zhang, H., Wang, W., Colombo, L., Wallace, R. M., and Cho, K.  Applied Physics Letters 103(5), 053513 (2013).
  * [42] Özçelik, V. O., Azadani, J. G., Yang, C., Koester, S. J., and Low, T.  Phys. Rev. B 94, 035125 Jul (2016).
  * [43] Kang, J., Tongay, S., Zhou, J., Li, J., and Wu, J.  Applied Physics Letters 102(1), 012111 (2013).
  * [44] Wang, Y., Wang, Z., Yao, W., Liu, G.-B., and Yu, H.  Phys. Rev. B 95, 115429 Mar (2017).