It is well known that stacking domains form in moiré superlattices due to the competition between the interlayer van der Waals forces and intralayer elastic forces, which can be recognized as polar domains due to the local spontaneous polarization in bilayers without centrosymmetry. We propose a theoretical model which captures the effect of an applied electric field on the domain structure. The coupling between the spontaneous polarization and field leads to uneven relaxation of the domains, and a net polarization in the superlattice at nonzero fields, which is sensitive to the moiré period. We show that the dielectric response to the field reduces the stacking energy and leads to softer domains in all bilayers. We then discuss the recent observations of ferroelectricity in the context of our model.

Twistronics, the study of layered systems with a relative twist angle or lattice mismatch between the layers, resulting in moiré superlattices, is one of the most exciting new topics in condensed matter physics. It was predicted about a decade ago that introducing a small relative twist in a layered system such as bilayer graphene could lead to flat electronic bands, and strongly correlated behavior^{1,2}. Moiré superlattices have since been shown to exhibit superconductivity^{3,4}, metal–insulator transitions^{5}, as well as magnetic^{6}, topological^{7–10}, and excitonic^{11,12} behavior, facilitated by the tuning of the twist angle or lattice mismatch. Recently, ferroelectricity was observed in bilayer graphene^{13} and hexagonal boron nitride (hBN)^{14}, which is highly unusual, because the constituent materials are nonpolar, and bilayer graphene is normally metallic. The ferroelectricity was found to be sensitive to the twist angle and lattice mismatch, with some samples exhibiting no hysteresis and some exhibiting strong hysteresis. The ferroelectricity is clearly very unconventional, and the physical mechanism is currently not well understood.

Structural phenomena in moiré superlattices are generally well understood. It is known that the interlayer separation ripples in space due to the local misalignment of the atoms, which can influence physical properties^{15,16}. In addition, lattice relaxation occurs due to the competition between the in-plane strains and out-of-plane van der Waals interactions, leading to stacking domains^{17–23}. The elastic energy depends on the twist angle and lattice mismatch quadratically, meaning the domains can be tuned. The domain structures have been shown to have a large influence on the properties of the system^{17–19,24–27}, leading to the opening of band gaps and enhanced Fermi velocity, for example. Polar effects have been given less consideration because the typical materials use to fabricate moiré superlattices, graphene, hBN, and transition metal dichalcogenides (TMDs) such as MoS_{2} (see Fig.

There are two main mechanisms by which polar phenomena can manifest in moiré systems. The first is a local spontaneous out-of-plane polarization^{28}, which occurs in bilayers without centrosymmetry and averages to zero over the moiré period. The second is couplings between strain and polarization, namely piezoelectricity^{29–31} and flexoelectricity^{29,32–37}. The strain gradient is largest across the domain walls, and via flexoelectricity, they have an inherent polarization. The flexoelectric response in 2D materials can be estimated by measuring the potential drop across the wall of a nanotube in the large radius limit^{29,38–40}, and it has been estimated that the flexoelectric coefficients in bilayer graphene are similar in magnitude to the clamped-ion flexoelectric response in oxide perovskites^{29,38,41}. The flexoelectric polarization is localized within the relatively narrow domain walls, however.

If we identify the stacking domains as polar domains via the two aforementioned mechanisms, then the stacking domains may serve as the basis for understanding polar phenomena in moiré materials. Thus, in order to understand the observed ferroelectricity, it is essential to understand how the stacking domains respond to an electric field. It is known that the domain structures in moiré materials can lead to interesting effects such as the opening of band gaps, and topologically protected states or channels when an electric field is applied^{7–9,42,43}. To our knowledge, the influence of an applied electric field on the domains themselves has not been considered. It is known that an electric field can modify the interlayer separation and lead to a breakdown of TMD bilayers, for example^{44,45}. The stacking domains are a result of lattice relaxation, which describes the delicate competition between the interlayer interactions and the intralayer elasticity. Since the interlayer interactions are sensitive to an applied field, it is reasonable to expect that the field would change the delicate balance and affect the resulting domain structure.

In this paper, we introduce a model of lattice relaxation in a moiré superlattice which includes the effect of an applied field on the bilayer. The total energy is an integral of the energy density over a moiré supercell_{sc} is the area of the supercell. For a bilayer system, Eq. (

We can model moiré superlattices at different levels of theory depending on the contributions we include in Eq. (

In order to consider the effect of an electric field on the atomic structure, we also include the electrostatic energy induced by an electric field

_{2} monolayer. The two stacking configurations are shown below: ^{∘}). _{2}. The dashed lines indicate the vertical positions of the atoms, and the vectors show the applied field _{2} along the out-of-plane direction averaged in the in-plane directions.

The stacking energy can be included in a number of ways. The simplest is to use the cohesive energy as a function of space, _{0}(^{15}. When considering the effect of an applied field, it is necessary to include the full van der Waals potential because some phenomena cannot be captured at the harmonic level, such as the breakdown of the bilayer for stronger fields^{44}. A detailed study of the stacking energy of bilayer MoS_{2} in the presence of an electric field, both theoretically and verified by first-principles calculations, is provided in Appendices A and B, respectively.

The first term in ^{28}. Bilayer systems without centrosymmetry, such as 3R MoS_{2} (Fig. _{2} (Fig.

The second term describes the dielectric response of the bilayer to the electric field, where _{0} and _{1} are the first two coefficients in the expansion of the polarizability

The typical lattice relaxation procedure is as follows: the local energy densities in Eq. (^{20}, where all of the local stacking configurations in real space are condensed into a single unit cell of relative translations between the layers (see Section _{2}), and translating one layer over the other. Quantities such as _{0}, etc., can be parameterized in configuration space, and Eq. (

The parameterization of _{0}, _{0}, _{1}, and _{0} for 3R and 2H MoS_{2} was done using ^{46} and is shown in Fig. _{2} has a small nonvanishing local dipole density with zero mean, see Fig. _{2} does not.

_{2}: _{0} (black) and _{1} (black) and _{0} (red); _{0} (black) and _{0} (red). The hollow markers are results from first-principles calculations, and the solid curves show the corresponding Fourier interpolations. Black refers to the leftmost vertical axis and red to the rightmost vertical axis. The configurations AA

The lattice relaxation at finite electric fields was then performed in configuration space, including an additional in-plane displacement ^{21,47,48}, with a single domain wall across the path AB → SP → BA (see Section

The 1D FK model, including only the dielectric response to the field, was solved for a fixed lattice mismatch

_{s}_{stack}(

The 2D relaxation including only the coupling between the field and spontaneous polarization was done for a range of twist angles 0. 1^{∘}≤ ^{∘} and field strengths ^{−1}. The results are summarized in Fig. _{2} in configuration space obtained from first-principles calculations, at electric field strengths of 0 and 2 V Å^{−1}, respectively. The electric field increases the depth of the well at AB and decreases the depth of the well in BA by the same amount, breaking the ^{∘}. At zero field, the relaxation reduces the area of the AA regions and increases the area of the AB/BA regions, leading to a triangular domain structure with sharp domain walls. When an electric field is applied, the AB and BA regions relax unevenly, leading to larger AB regions and smaller BA regions reducing the rotation symmetry to _{M} as a function of twist angle, which was obtained by taking the slope of the polarization about zero field. We can see that the susceptibility increases dramatically as the twist angle decreases.

_{2}, including the linear electrostatic term. The top panels show the stacking energy before lattice relaxation at electric field strengths of (^{−1} and (^{−1}. ^{∘}. The AA regions (yellow) shrink and the AB/BA regions (purple) expand. When an electric field is applied, the AB and BA regions relax unevenly, one increasing in area and the other decreasing. _{M} as a function of twist angle, obtained from the slopes from (

We have introduced a model which illustrates the effect of an applied electric field on lattice relaxation in moiré superlattices. The model contains two electrostatic contributions. The first is a linear coupling between the field and the local spontaneous polarization in bilayers without centrosymmetry, which breaks the degeneracy between the AB- and BA-stacking domains. Under an electric field, the AB and BA regions will relax unevenly with one growing and the other shrinking with respect to the relaxed structure at zero field. This leads to a nonzero average out-of-plane polarization in the superlattice. The second contribution is the dielectric response to the field, which occurs in all bilayers. This term leads to a nonuniform increase the layer separation, which reduces stacking energy, leading to softer domains structures under lattice relaxation.

Finally, as our theory does not predict a ferroelectric response, we offer some thoughts on the recent experimental observations of ferroelectricity in the context of our model.

For a system to be considered ferroelectric, it (i) must exhibit a spontaneous polarization at zero field which (ii) must be switchable with an electric field. However, neither the individual stacking domains nor the moiré superlattice as a whole satisfies both conditions: The stacking domains indeed have a local spontaneous polarization at zero field, and while the average polarization of a domain can change via lattice relaxation under an electric field, the sign of the polarization in each cannot be switched. Therefore, the stacking domains in moiré superlattices are in general not ferroelectric. Conversely, the moiré superlattice itself exhibits an average polarization, the direction of which can be changed by the field, but has zero average polarization at zero field. Therefore, under ideal conditions, moiré superlattices are also not ferroelectric.

This idealized picture may not hold in experimental settings, and defects, mislocations or strain induced by the finite size of samples may lead to uneven domains at zero field. Also, the direct couplings between strain and polarization, piezoelectricity and flexoelectricity, have not been considered, which may make it energetically favorable for the domains to relax unevenly and the superlattice to have a nonzero average polarization at zero field. We leave the consideration of these effects for future work. To summarize, for an ideal system, when considering the spontaneous local polarization and lattice relaxation under an electric field, neither the moiré superlattice nor the stacking domains are ferroelectric, since the former does not have a spontaneous polarization at zero field and the latter does not have a switchable polarization.

There have also been reports of a switching of the polarization in a single stacking domain by a sliding of the atoms by half a monolayer lattice constant when a local field was applied to the domain using biased atomic force microscopy (AFM) tip^{14}. This sliding change the stacking configuration: AB ↔ BA, leading to a first-order switching of the polarization. This is a separate mechanism to the one mediated by lattice relaxation, which results in a second-order change in the polarization. We can understand this sliding in the context of our model. When a field is applied to the domain, the linear coupling between the field and polarization will lower the energy if the two are aligned and result in a large energy penalty if they are anti-aligned. In either case, the dielectric response will reduce the stacking energy by the same amount, making it easier for one layer to slide with respect to the other. When the field and polarization are anti-aligned, the energy can be lowered considerably via a sliding by half a monolayer unit cell, flipping the polarization so that it becomes aligned with the field. However, the field is applied to the domain locally via an AFM tip, and it is not clear whether the sliding occurs locally under the tip, or throughout the entire domain. It is also not clear whether or not the domain will remain flipped once the field is removed, or relax back to its original orientation. Thus, it is not clear whether or not this mechanism for a first-order switching of polarization in a stacking domain is truly ferroelectric either.

The model introduced in this paper illustrates, clearly and intuitively, the effect an electric field can have on lattice relaxation in moiré superlattices. We propose an electric field as a third quantity with which the domain structures in moiré superlattices can be tuned. Unlike the twist angle and lattice mismatch which are fixed for a given sample, an electric field can be applied dynamically to tune a sample in situ. Thus, it may serve as a more practical approach to achieve control in moiré superlattices. We have also discussed how our theoretical model can be used to understand the recent observations of ferroelectricity in moiré superlattices. We believe it is inaccurate to consider moiré materials to be truly ferroelectric via lattice relaxation or sliding under an electric field. However, this motivates further study into polar phenomena in moiré materials.

First-principles density functional theory (DFT) calculations were performed using the ^{46} using PSML^{49} norm-conserving^{50} pseudopotentials, obtained from pseudo-dojo^{51}. ^{46,52}, and double-^{53}.

A mesh cutoff of 1200Ry was used for the real space grid in all calculations. A Monkhorst–Pack ^{54} of 12 × 12 × 1 was used for the initial geometry relaxations, and a mesh of 18 × 18 × 1 was used to calculate polarizabilities. Calculations were converged until the relative changes in the Hamiltonian and density matrix were both less than 10^{−6}. For the geometry relaxations, the atomic positions were fixed in the in-plane directions, and the vertical positions and in-plane stresses were allowed to relax until the force on each atom was less than 0.1 me V Å^{−1}. The layer separation _{2} (see Fig _{0}, was obtained by fixing the relaxed geometry, applying an electric field large enough to overcome internal field effects, and measuring the change in the out-of-plane dipole moment of the bilayer. The polarizability to first order in _{1}, was obtained by changing the layer separation by ± 1% with respect to _{0} and measuring the relative change in the polarizability. A detailed first-principles study is provided in Section

Equation (_{G}_{G}^{iG⋅s}, where _{G}}, and can be minimized numerically with respect to the coefficients:

We can take advantage of the _{2} at zero field, we have _{−G} = −_{G}, and the cosine terms vanish. The optimization was done using the independent ^{∘}, six shells (21 vectors) for 0. 5^{∘} > ^{∘} and seven shells (28 vectors) for 0. 3^{∘} > ^{∘}. The total energy was optimized using the limited memory BFGS algorithm (L-BFGS), until the gradient was below 1 × 10^{−5} eV.

D.B. would like to thank E. Artacho and I. Lebedeva for helpful discussions. D.B. acknowledges funding from the EPSRC Centre for Doctoral Training in Computational Methods for Materials Science under grant number EP/L015552/1. B.R. gratefully acknowledges support from the Cambridge International Trust.

D.B. conceived the project and performed the first-principles calculations. D.B. and B.R. developed the theoretical model and performed the lattice relaxation calculations. D.B. wrote the manuscript, with contributions from B.R.

The data presented were generated from first-principles and lattice relaxation calculations as described in the text.

Code is available upon reasonable request.

The authors declare no competing interests.

Spplementary material

The online version contains supplementary material available at