Electrically tunable stacking domains and ferroelectricity in moir\'e superlattices

It is well known that stacking domains form in moir\'e 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\'e 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.


INTRODUCTION
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 , metalinsulator 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 non-polar, 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 . Additionally, 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][18][19][20][21][22][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][18][19][24][25][26][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 Figs. 1 (a)-(c)), are non-polar.
There are two main mechanisms by which polar phenomena can manifest in moiré systems. The first is a local spon- * db729@cam.ac.uk E P Eint taneous 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][30][31] and flexoelectricity 29,[32][33][34][35][36][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][39][40] , and it has been estimated that the flexoelectric coefficients in bilayer graphene is 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 topo-logically protected states or channels when an electric field is applied [7][8][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.

RESULTS
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 where A sc is the area of the supercell. For a bilayer system Eq. (1) is a discrete sum over atomic sites, but generalizes to a continuum field theory when the moiré period is much larger than the lattice constants of the monolayers. We can model moiré superlattices at different levels of theory depending on the contributions we include in Eq. (1). The stacking energy V stack captures the weak van der Waals interactions between the layers in terms of the layer separation d. The elastic energy V elastic allows for in-plane displacements of the atoms in the layers. Together, the stacking and elastic energies provide a good description of the atomic structure in moiré superlattices, namely lattice relaxation and the formation of stacking domains. Having obtained a realistic description of the structure, one could proceed to obtain the electronic bands from tight-binding theory.
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 E perpendicular to the bilayer, see Fig. 1 (d). The total energy density is then Each contribution is derived and discussed in detail in Section I of the supplementary information. In the elastic energy, summation is assumed, C is the linear elasticity tensor and ε i j = 1 2 (∂ i U j (r) + ∂ j U i (r)) is the strain tensor, written in terms of a relative in-plane displacement U.
The stacking energy can be included in a number of ways. The simplest is to use the cohesive energy as a function of space, V 0 (r), assuming that at each point in the supercell the layer separation takes the value that minimizes the local stacking energy: d(r) = d 0 (r). Other studies have allowed the layer separation to vary by performing a harmonic expansion about the equilibrium layer separation 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 V elec is the coupling between the electric field and the out-of-plane spontaneous dipole moment of the bilayer 28 . Bilayer systems without centrosymmetry, such as 3R MoS 2 ( Fig. 1 (b)), have a local dipole moment throughout the supercell which averages to zero, whereas systems with centrosymmetry have no local dipole moment anywhere in the supercell, such as 2H MoS 2 ( Fig. 1 (c)).
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 α about the equilibrium layer separation. A bilayer system cannot simply be treated as a pair of capacitor plates; due to the overlap of states in the vacuum region between the layers, it is more appropriate to treat the system as a single slab with a nonuniform charge density, see Fig. 1 (e). Thus, changing the layer separation will affect the polarizability of the system, so we perform a Taylor expansion in d. In Section II of the supplementary information we show with first-principles calculations that the polarizability is linear in d. In addition, the polarizability will vary throughout the superlattice due to the different local stacking configurations and equilibrium layer separations. The dielectric response occurs in all layered systems, irrespective of symmetry.
The typical lattice relaxation procedure is as follows: the local energy densities in Eq. (2) are parameterized using firstprinciples calculations. Practically, this is done using the mapping between real space and 'configuration space' 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 II of the supplementary information). Configuration space can be traversed using first-principles calculations with only a single primitive cell (6 atoms for bilayer MoS 2 ), and translating one layer over the other. Quantities such as V 0 , d 0 , etc., can be parameterized in configuration space, and Eq. (1) can then be minimized with respect to the layer separation d and in-plane displacements U in order to obtain the relaxed structure.
The parameterization of V 0 , d 0 , α 0 , α 1 and P 0 for 3R and 2H MoS 2 was done using SIESTA 46 and is shown in Fig. 2. Starting from the metal over metal configuration (AA), one layer was fixed and the other was translated along primitive cell diagonal in small increments, and the aforementioned quantities were measured at each point. By taking advantage of the C 3 rotation symmetry of the moiré superlattices, the data were interpolated by a low order 2D Fourier expansion throughout configuration space, greatly reducing the number of calcula-  tions required. The stacking energy and equilibrium layer separation both vary by about 1Å / 0.1 eV, respectively, which is expected. We also found that both polarizability parameters vary significantly throughout the configuration space. Additionally, 3R MoS 2 has a small nonvanishing local dipole density with zero mean, see Fig. 2 (e), whereas 2H MoS 2 does not.
The lattice relaxation at finite electric fields was then performed in configuration space, including an additional inplane displacement u(s) in the stacking and electrostatic energies, and minimizing the total energy numerically. The two terms in V elec were included separately in order to illustrate their individual effects. The quadratic term does not break any symmetries, and thus we can study its effect on the domain structures by mapping the 2D moiré superlattice to a 1D Frenkel-Kontorova (FK) model 21,47,48 , with a single domain wall across the path AB → SP → BA (see Section III of the supplementary information for more details). The linear term breaks the C 6 rotation symmetry of the 3R moiré superlattice, leading to a splitting between the AB and BA domains, which cannot be captured by a 1D FK model.  The 1D FK model, including only the dielectric response to the field, was solved for a fixed lattice mismatch θ , and for several values of electric field, see Fig. 3, with similar results at zero field for various lattice mismatches shown in Fig S10. At zero field, decreasing θ increases the displacement u(s) in between AB/BA and SP points, leading to a domain structure with wide AB/BA regions and narrow SP regions, separated by a domain wall, with width proportional to θ in configuration space. The domain structure is also evident from the stacking energy and equilibrium layer separation profile in configuration space, shown in Figs. 3 (c) and (d). When an electric field is applied, the equilibrium layer separation increases everywhere in configuration space, which decreases the stacking energy. With the stacking energy reduced, it is not as favorable for the atoms to relax; the displacement decreases, and the domain walls soften.
The 2D relaxation including only the coupling between the field and spontaneous polarization was done for a range of twist angles 0.1 • ≤ θ ≤ 1.0 • and field strengths 0 ≤ E ≤ 2 VÅ −1 . The results are summarized in Fig. 4, and additional plots for all angles and electric field strengths can be found in Section IV of the supplementary information. In Figs 4 (a) and (b) we show the stacking energy in 3R MoS 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 C 6 rotation symmetry. The panels below, Figs 4 (c) and (d), show the corresponding stacking energies after lattice relaxation for a twist angle of θ = 0.1 • . At zero field, the relaxation re-  duces 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 C 3 . When the AB and BA domains are no longer equal in area, the polarization no longer averages to zero. In Fig. 4 (e) we show the average spontaneous polarization P 0 (s) in configuration space as a function of field for different twist angles. We can see that the response to the field is very sensitive to the twist angle. In Fig. 4 (f) we show the susceptibility of the moiré superlattice χ 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.

DISCUSSION
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 satisfy 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 favourable 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 a 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 calculations
First-principles density functional theory (DFT) calculations were performed using the SIESTA code 46 using PSML 49 norm-conserving 50 pseudopotentials, obtained from pseudodojo 51 . SIESTA employs a basis set of numerical atomic orbitals (NAOs) 46,52 , and double-ζ polarized (DZP) orbitals were used for all calculations. The basis sets were optimised by hand, following the methodology in Ref. [53].
A mesh cutoff of 1200 Ry was used for the real space grid in all calculations. A Monkhorst-Pack k-point grid 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 inplane stresses were allowed to relax until the force on each atom was less than 0.1 meVÅ −1 . The layer separation d was taken to be the distance between the carbon atoms in bilayer graphene an the distance between the metals in bilayer MoS 2 (see Fig 1 (e)), and the stacking energy is calculated as V stack = V bilayer −2V mono , where V bilayer and V mono are the total energies of the bilayer and monolayer systems, respectively. The polarizability to zeroth order in d, α 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 d, α 1 , was obtained by changing the layer separation by ±1% with respect to d 0 and measuring the relative change in the polarizability. A detailed firstprinciples study is provided in Section II of the supplementary information.

Lattice Relaxation
Eq. (1) can be minimized by using variational methods and solving the resulting differential equations. This can be relatively demanding in 2D. Instead, we minimized the total energy using numerical optimization methods. Eq. (1) is a functional of the in-plane displacement u in configuration space. If we perform a plane-wave expansion, u(s) = ∑ G u G e iG·s , where G are the reciprocal lattice vectors of the commenusrate bilayer, then the total energy becomes a function of {u G }, and can be minimized numerically with respect to the coefficients: ∇ U G V tot = 0. This was done in JULIA, using the OPTIM package to do the optimization.
We can take advantage of the C 3 symmetry of our model to greatly reduce the number of independent G vectors: When there is a C 6 symmetry, i.e. for 3R MoS 2 at zero field, we have u −G = −u G , and the cosine terms vanish. The optimization was done using the independent G in the first five shells (10 vectors) for θ ≥ 0.5 • , six shells (21 vectors) for 0.5 • > θ ≥ 0.3 • and and seven shells (28 vectors) for 0.3 • > θ ≥ 0.1 • . The total energy was optimized using the limited memory BFGS algorithm (L-BFGS), until the gradient was below 1 × 10 −5 eV.
Data Availability. The data presented were generated from first-principles and lattice relaxation calculations as described in the text. Code Availability. Code is available upon reasonable request. A twisted bilayer system is composed of two layers with a relative twist angle θ between them, which we call the reference (r) and twisted (t) layers, respectively. The lattice vectors of the reference layer and the twisted layer are respectively, where a is the monolayer lattice constant and R θ = cos (θ ) − sin (θ ) sin (θ ) cos (θ ) . For a general θ , the two layers are incommensurate, i.e. they form a supercell which is infinitely large. If the two layers form a commensurate supercell we can define supercell lattice vectors L i as a linear combinations of the lattice vectors of either layer: L 1 = pa r,1 + qa r,2 =pa t,1 +qa t,2 =pR θ a r,1 +qR θ a r,2 L 2 = R π/3 L 1 . ( The second line of Eq. (2) leads to a Diophantine equation in the integers p, q,p andq, and the set of angles which result in a commensurate supercell is given by [1][2][3] θ (p, q) = cos −1 Consider an atom in the reference layer, with in-plane position r 0 . The corresponding atom in the twisted layer has position r = R θ r 0 . The displacement due to twisting is δ (r) = r − r 0 = I − R −1 θ r. If the layers form a commensurate supercell with lattice vectors L M i , say, the atoms will realign when the displacement is equal to a lattice vector, δ (L M i ) = a r,i : The cell spanned by L M i is known as a moiré superlattice, and is known as the the moiré period, which is is not necessarily equal to the supercell period: For p = q, there will be (p − q) 2 moiré periods in the supercell.
Having established a geometric description, we can model structural or electronic phenomena in a twisted bilayer using phenomenological or tight-binding models. These models can be paramaterized using first-principles calculations. However in practice this is difficult to do in real space because the size of and the number of atoms in the supercell becomes prohibitively large at smaller twist angles. Fortunately, we can take advantage of a useful mapping which allows us to parameterize systems at arbitrary twist angles using just a single cell of a commensurate bilayer.
The set of displacements of every atom in a moiré superlattice can be described by δ (r). Alternatively, we can describe the displacements by a set of local translations i.e. at r, the system is locally equivalent to an untwisted bilayer, with a relative in-plane slide s θ (r) between the layers, see Fig. 1.
The set of translations s(r) is contained in a single primitive cell of the reference layer, even in the continuum limit. We call this space of translations configuration space [4][5][6][7] . This space can be traversed by taking an untwisted commensurate bilayer and sliding one layer over the other. Physical properties can be measured in configuration space, and we can use the inverse map, to parametrize a moiré superlattice in real space for arbitrary twist angles 7 . Derivatives in real space and configuration space are related by which is useful for mapping quantities which depend on spatial derivatives, such as strains, to configuration space.

Stacking energy
It is well known layered systems interact via long-range van der Waals forces. For non-polar materials, this is facilitated by induced dipole−induced dipole, or London interactions, where fluctuations in the charge density of one layer lead to a dipolar response in the other layer, and vice versa. This can be described by a van der Waals potential in the layer separation d: where V 0 is the depth of the potential well or cohesive energy per unit cell, d 0 is the equilibrium separation, and the indices (n, m) determine the curvature of the well about the minimum. Eq. (10) can be parametrized using first-principles calculations, which is done in Appendix B. For two slabs which extend infinitely in area, the energy is expected to behave like d −2 at larger separations, since the dipole-dipole interaction is integrated twice over an infinite area 8 . From first-principles calculations we found that this is not the case, and the long-range interactions in bilayer graphene and MoS 2 decay with n ≥ 6. We suspect that this is because a bilayer is not adequately described as a pair of capacitor plates, but rather a single slab with a non-uniform charge density as in Fig. 1 (e). The non-zero overlap of states in the vacuum region may be screening the long-range interactions and could be responsible for larger than expected values of (n, m), although a more detailed study is required to verify this.

Electrostatic energy
The electrostatic energy of a dielectric slab in the presence of a perpendicular electric field E is where P and P 0 are the total and spontaneous polarization, respectively, A is the in-plane area of the bilayer, Ω = Ad is the volume of the bilayer, and χ(d) is the dielectric susceptibility, which in general depends on the geometry of the system. Eq. (11) assumes a linear response to the applied field, P(E, d) = P 0 + ε 0 χ(d)E, which inserting into Eq. (11) gives: where α = Ωχ is the polarizability, which describes the linear response of the dipole moment to the applied field, and p 0 = ΩP 0 is the spontaneous dipole moment. When an electric field is applied, the layers are no longer non-polar, and the long-range interactions can no longer be considered induced dipole−induced dipole interactions. However, we can think of the interactions between the layers as arising due to fluctuations about a finite polarization, rather than zero polarization. The effect of an applied field can be estimated from first-principles by measuring the stacking energy as a function layer separation for different field strengths. In Ref. 9, it was found that the electric field lowers the stacking energy of bilayer MoS 2 for large separations, leading to a breakdown of the bilayer at E ∼ 2 VÅ −1 , although the precise behavior of the equilibrium layer separation as a function of applied field is unclear. In order to clarify this, we performed detailed first-principles calculations in Appendix B. In Fig. 5 (a) we plot the potential energy curves for bilayer MoS 2 as a function of d for different values of applied field. We see that the electric field lowers the energy at larger separations, and the bilayer becomes unstable at E crit ≈ 2.25 VÅ −1 , similar to the results obtained in Ref. 9. By calculating the Mülliken charges on each layer in Fig. 5 (b), we find that there is an interlayer charge transfer when an electric field is applied. The charge transfer is linear in E above E ≈ 0.27 VÅ −1 , below which little or no charge transfer is observed. We attribute this to internal field effects, similar to the depolarizing field observed in ferroelectric thin films. see Fig. 1 (d). The internal field could be accounted for in the electrostatic energy as 9 V int = − 1 2 E int P. In the interest of simplicity, we neglect internal field effects in our model.
In Fig. 6 (a) we plot the dipole moment p(E) for several values of d. For field strengths strong enough to overcome internal field effects, the dipole moment is linear in E, but also in d. Thus, as shown in Fig. 6 (b), the polarizability is approximately linear in d, where α 0 and α 1 are both positive and have units ofÅ 3 . Taking (n, m) = (n, 2n), we can write the energy, exlcuding elastic contributions, as By plotting contours of the energy as a function of d for fixed values of E, see Fig. 2 (a), we can see that the stacking energy is lowered by the electric field, consistent with results from first-principles calculations. At zero field, V stack has a minimum at the equilibrium separation d min = d 0 , and a maximum (for d > d min ) at d max → ∞. When a field is applied, V stack diverges as d → ∞. This makes physical sense: increasing d will increase the total dipole moment and lower the total energy. For a non-zero field, d max has a finite value, at the top of the energy barrier which separates d min and d → ∞. As the field strength increases, d min and d max move closer together, eventually meeting at a critical point (d crit , E crit ) where the energy barrier vanishes, and the bilayer becomes unstable, where d ≡ d d 0 and E ≡ ε 0 α 1 d 0 |V 0 | E. In Fig. 2 (b) we show d min and d max as a function of E for several values of n. The solid red line shows Eq. (15), which separates d min and d max . Interestingly, the behavior is highly sensitive to the value of n. For larger values of n, which are predicted from first-principles calculations, d min does not change by much until E ∼ E crit , where it increases by about 10% before the bilayer becomes unstable. For smaller values of n, the layers separate more easily because the potential well is shallower, but the breakdown occurs at smaller field strengths.

Elastic energy
The elastic energy in Eq. 2 is given as a contraction of the strain tensors with the fourth rank elasticity tensor 7,10 : and is given explicitly as The strain tensor in configuration space is 7 which allows us to write down the elastic energy in configuration space: It is convenient to work in terms of the lattice vectors a r,1 and a r,2 . Under this transformation, the displacement transforms as u → gu and the strain tensor transforms as ε i j → g −1 iα ε αβ g β j , where g = 1 1/2 0 √ 3/2 . The elastic energy is then ∂ a x u a y + ∂ a y u a x 2 + ∂ a x u a x − ∂ a y u a y 2 (20)

SUPPLEMENTARY NOTE 2: FIRST-PRINCIPLES CALCULATIONS
Because functionals within the generalized gradient approximation (GGA) tend to underestimate the cohesive energy of van der Waals systems, the different VDW-corrected functionals in SIESTA were tested by measuring the stacking energy of bilayer graphene as a function of layer separation, see Fig. 3. The rest of the calculations were performed using the C09 functional, as it was found to give good results for both bilayer graphene and MoS 2 .

Parametrizing the stacking energy
We can parameterize the stacking energy by calculating V 0 , d 0 and the two indices (n, m). Using suitable values for (n, m), a fit to the stacking energy as a function of layer separation was obtained using Eq. (10), which is shown in Fig. 4 for bilayer graphene and MoS 2 with 3R stacking. The parameters used to fit Eq. (10) to the first-principles calculations are given in Table I.  To parametrize V elec , an electric field was applied in the out-of-plane direction. A dipole correction [20][21][22][23] was used in the vacuum region to prevent long-range interactions between periodic images. The resulting dipole moment in the out-of-plane direction was then measured.
We first performed geometry relaxations of bilayer MoS 2 for electric fields of increasing strength. We found that the layer separation increases only marginally (a very low force tolerance is required to see this, otherwise the layers do not move) until around E ∼ 2 VÅ −1 , where the layers begin to separate and the bilayer quickly becomes unstable. In order to investigate this peculiar behaviour, we calculated the stacking energy curves for various electric field strengths, see Fig. 5 (a). We can see that  Table I. The dashed lines show the same parameterization, but with the smaller index changed to 2.
the electric field lowers the stacking energy, in agreement with our model. In Fig. 5 (b) we show the Mülliken chares as forces on each layer as a function of electric field, which indicates that there is a transfer of charge across the layers as an electric field is applied, verifying similar results obtained in Ref. 9.
For each value of d, we calculate the dipole moment as a function of electric field, p(E), as shown in Fig. 6 (a). We see that, at smaller values electric field strength, E < 0.27 VÅ −1 , changing the distance has little effect on the polarizability α. As mentioned previously, we attribute this to internal field effects. In order to neglect these effects, we measure the polarizability α(d) at field strengths which are large enough to overcome the internal fields. At stronger values of field α(d) increases linearly with d, see Fig. 6 (b).

Parameterization in configuration space
Calculations were repeated to parameterize V 0 , d 0 , α 0 and α 1 as a function of s. The high symmetry stacking configurations AA (metal over metal), AB/BA (metal over chalcogen), and the saddle point (SP) all lie along the configuration space diagonal. Thus, it was sufficient to perform a series of calculations along the diagonal, and use a 2D Fourier interpolation to parameterize the model everywhere in configuration space. One layer was fixed, and the other layer was translated by s = s(a 1 + a 2 ), s ∈ [0, 1], and the aforementioned quantities were measured on a fine grid of values of s. The interpolation was done following similar approaches in previous studies 7, [24][25][26] : each quantity is written as a Fourier expansion. The results from the first-principles calculations at different values of s are used to fit the Fourier coefficients, and a smooth interpolation of each quantity is obtained everywhere in configuration space. The reciprocal lattice vectors G of the same length are sorted into shells, and the first few shells are sufficient to obtain good parameterizations. The first-principles results and parameterization using the first three shells are shown in Fig. 4 for 3R and 2H stacked bilayer MoS 2 .
After parameterizing the model, we can examine the effect of an electric field on the system in configuration space. Fig. 7  (a) shows the layer separation in configuration space at different electric field values for 3R MoS 2 . Without lattice relaxation, only the quadratic electrostatic term affects the layer separation. We can see that the layer separation increases non-uniformly in configuration space, leading to a corresponding reduction in stacking energy, shown in Fig. 7 (b). We only show up to E = 1.68 VÅ −1 because that is the smallest critical field at which d min → ∞ somewhere in configuration space. We examine this in more detail in Fig. 7 (c) by showing d min as a function of electric field at the AA, AB and SP points. We can see that the critical field values for the AB and SP points are considerably larger. We could reach stronger field strengths by including terms proportional to ∇d in the elastic energy from von Karman plate theory 24 , but this would make the model considerably more difficult to solve, requiring either solutions to fourth order differential equations or a more difficult optimization problem. After lattice relaxation, the AA regions shrink considerably, and the critical fields of interest for the domain wall are E SP crit and E AB crit ,   since the domain wall is across the path AB → SP → BA. However, without including nonlinear terms in the elastic energy, we are limited by the smallest critical field, E AA crit . In Fig. 8 (a) we show the spontaneous polarization P 0 (s) along the diagonal in configuration space. The coupling between the polarization and the field leads to an even splitting of the AB and BA wells, increasing one and decreasing the other, which is shown in Fig. 8 (b).     The Frenkel-Kontorova (FK) model was first introduced in the 1930s and is widely used in condensed matter physics 27,28 . It is a discrete model where a chain of atoms is subject to a rigid periodic potential, from a substrate for example 27,[29][30][31] . The FK model has been used as a 1D representation of twisted bilayer systems 5,25,32 , replacing the rigid substrate with a second chain of atoms which can also deform, see Fig. 9, and has been used to describe commensurate-incommensurate phase transitions 32,33 , lattice relaxation and domain structures 7,25,34 .
Because the dielectric response does not break any symmetries, we study its effect on the domain walls using a 1D FK model. For 3R MoS 2 , the domain walls are along the path AB → SP → BA, so we use those points to parameterize the stacking energy 32,33,35 in a one dimensional version of configuration space, described by a single variable s, the relative displacement between the atoms in the two chains, and similarly for d 0 , α 0 and α 1 . The total energy in configuration space is where θ is the lattice mismatch, B is the bulk modulus of MoS 2 7 and V stack and V elec include the displacement u(s) to allow for relaxations. Eq. 22 can be also be minimized using optimization techniques, but in 1D the differential equations are easy enough to solve. Minimizing with respect to both u and d, we get The second equation can be solved independently to obtain d min (s, E). This is inserted into the first equation, which can be solved numerically.
In Fig. 10 we show results for various values of θ at zero field. For larger θ , the atoms do not move much from their initial stacking configuration. As θ decreases, the elastic energy is reduced, and the atoms relax more. We see from Fig. 10 (a), the total displacement s + u(s), that the atoms move to maximize the area around s = 1 2 , the stacking configuration with lowest energy. Figs of which is proportional to θ . Fig. 10 (b) shows the change in displacement ∂ s u(s), which is proportional to the strain tensor in 1D, from which we see that there is a large strain gradient across the domain wall. At s = 0, we have the analytic approximation 32 : which is shown alongside the numerical solutions to Eq. (23) in Fig. 11 (a). We can approximate the width of the domain wall as w ∼ 1 2 illustrating the dependence on both twist angle and electric field. In Fig. 11