Lattice Relaxation in Moiré Heterobilayers

Abstract

We develop an analytical theory for lattice relaxation in twisted moiré heterobilayers, accounting for lattice mismatch, twist, external biaxial heterostrain, and different elastic constants. Starting from continuum elasticity, we derive the self-consistent equations for the in-plane displacement fields and obtain simple perturbative expressions for the layer-resolved in-plane displacement fields induced by lattice relaxation. We apply our theory to graphene on hBN and representative 2H transition metal dichalcogenide heterobilayers, including MoTe₂/WSe₂ and WSe₂/WS₂. Our analytical results agree very well with full numerical solutions over experimentally relevant parameters. We further show that heterobilayers can exhibit a buckling instability near alignment, driven by compressive in-plane strain due to moiré relaxation. Our results provide a simple theoretical framework for incorporating lattice relaxation in realistic moiré heterostructures.

I Introduction

Van der Waals heterostructures provide a versatile platform for discovering both new physics [4, 3, 39] and novel technologies [45, 2, 32]. Already more than 10 years ago it was realized both experimentally [60] and theoretically [29, 30] that stacked atomically thin two-dimensional (2D) materials behave more like elastic membranes than rigid crystals. Hence, the structural and electronic properties of moiré superlattices, formed by stacking layers with a relative twist or lattice mismatch, can be significantly modified by lattice relaxation. In particular, this gives rise to important new physics, including the separation of the flat bands from dispersive bands in magic-angle twisted bilayer graphene [43, 10] and the opening of a topological gap in graphene on hexagonal boron nitride (hBN) [29, 1, 28]. This stabilizes an anomalous Hall phase [7] resulting in the observation of the orbital anomalous Hall effect [52, 51].

Atomic reconstruction in moirés is controlled by the competition between intralayer elastic energy and interlayer stacking energy [60, 29, 11], becoming increasingly important as the moiré period increases. Much of the work on lattice relaxation over the past decade has been heavily numerical, either by solving the fully relaxed density functional theory [36, 42, 9, 38] or molecular dynamics simulations [36]. These methods become unfeasible for small twist angles due to the large number of atoms in the moiré unit cell ($\sim 1/\theta^{2}$) and provide limited physical insights into the relaxation process. However, very recently there have emerged analytical treatments of moiré reconstruction, including a large-angle theory for both monoatomic hexagonal lattices [21, 31] (e.g. twisted bilayer graphene) and binary hexagonal lattices [e.g. transition metal dichalcogenides (TMDs)] [21, 61] that agree well with available numerics. Even more recently, we established an approximate analytical treatment [14] of lattice relaxation for marginal twist angles, matching predictions from ab initio methods [62] and numerical solutions of continuum elasticity [43, 11, 6]. In particular, we find an emergent universality where the theory is characterized by a single twist-angle dependent parameter that captures both the small and large angle regimes. It enables accurate predictions of the acoustic displacement fields of any homo moiré interface built from layers with $D_{3h}$ symmetry (twist near $0^{\circ}$) or $D_{3d}$ (twist near $60^{\circ}$) and encompasses a vast number of moiré materials of contemporary interest. Examples of $D_{3h}$ homobilayers are: twisted bilayer graphene, 2H TMDs (MoS₂, WSe₂, MoTe₂, WS₂, NbSe₂, etc.), and 2D magnets like 2H Fe₃GeTe₂. Some examples of $D_{3d}$ homobilayers are given by: twisted double bilayer graphene, 1T TMDs (ZrS₂, HfS₂, HfSe₂, TiS₂, etc.), and 2D magnets e.g. NiI₂, MnI₂, and FeCl₂, and CrX₃ systems such as CrI₃.

However, the previous work focused only on homobilayers, where the two twisted sheets are made from the same material. Heterobilayers, comprising twisted sheets of two different materials, are equally interesting: for example, it is well-known that graphene on hBN features a Hofstadter’s butterfly and the fractional quantum Hall effect [17, 48] and strong moiré potentials [34]. Also, in Ref. [37] a topological phase transition from a Mott insulator to a quantum anomalous Hall insulator was observed in MoTe₂/WSe₂ moiré heterobilayers and Ref. [50] observed Mott and generalized Wigner crystal states in WSe₂/WS₂ moiré superlattices. And theoretically, the formation of quantum dots localized at the highly strained nodes of domain wall networks in MoX₂/WX₂ heterostructures have been proposed [54].

The goal of this paper is to generalize the analytical framework [21, 14] to heterobilayers where lattice relaxation is expected to be important even without twisting when the lattice mismatch between inequivalent layers is small. Here, we successfully develop a one-shot perturbation theory for the in-plane acoustic displacement fields that works for any heterobilayers with the same Bravais lattice allowing for lattice mismatch between layers with different in-plane point group symmetry, e.g. for graphene on hBN: graphene has $C_{6v}$ symmetry, while hBN has $C_{3v}$ symmetry, with approximately $2\%$ mismatch between lattice constants. We apply our analytical theory to several moiré systems of interest such as graphene on hBN and 2H TMDs, including MoTe₂/WSe₂, WS₂/WSe₂, and WS₂/MoSe₂, and compare our work with numerical results. In general, the analytical theory is accurate within $10\%$ for any twist angle. Beyond the in-plane results, we also discuss the role of out-of-plane displacements in free-standing heterobilayers, and identify a buckling transition driven by moiré reconstruction. Our current work extends the analytical framework of lattice relaxation to a much larger class of Van der Waals heterostructures.

Figure 1: First-star analytical theory for (a) the relative in-plane displacement field $\bm{u}=\bm{u}_{1}-\bm{u}_{2}$ and (b) the layer center-of-mass field $\bm{U}=\bm{u}_{1}+\bm{u}_{2}$, for graphene on hBN due to lattice relaxation. Here, the magnitude is shown in units of the average lattice constant $a=(a_{1}+a_{2})/2$ with $a_{1}>a_{2}$. (c) First-star stacking potential with $C_{3v}=\left<\mathcal{C}_{3z},\mathcal{M}_{x}\right>$ symmetry in configuration space from Eq. (22) for $\varphi_{1}=50^{\circ}$. (d) Longitudinal Fourier components of $\bm{u}$, where the dots are numerical solutions of the self-consistent equations given in Eqs. 14 and 15, and the curves give the one-shot result. The legend indicates the star following the notation of Table 1. Shown as a function of lattice mismatch $\epsilon$ for $\theta=0$. (e) Longitudinal and (f) transverse components versus twist angle for $\epsilon=0.0187$. Components that are not shown are negligible. Fourier components of the layer center-of-mass field are obtained from $U_{\bm{b}}^{\parallel}=A^{\parallel}u_{\bm{b}}^{\parallel}$ and $U_{\bm{b}}^{\perp}=A^{\perp}u_{\bm{b}}^{\perp}$ with $A^{\parallel}\approx A^{\perp}\approx 0.08$. Parameters used are shown in Tables 2 and 3.

II Structural relaxation in moiré bilayers

We start from the total structural energy of the bilayer system, which is composed of elastic energy, which favors the rigid moiré without deformations, and adhesion energy, which favors deformations that increase regions with energetically-favorable layer stacking. Treating the layers as isotropic continuous membranes, the elastic energy is given by [35]

$$H_{\text{elas}}=\frac{1}{2}\sum_{l=1,2}\int d^{2}\bm{r}\left[\lambda_{l}u_{ii}^{(l)}u_{ii}^{(l)}+2\mu_{l}u_{ij}^{(l)}u_{ji}^{(l)}\right],$$ (1)

where the sum runs over layers $l=1,2$. Here, $\lambda_{l}$ and $\mu_{l}$ are in-plane Lamé constants for layer $l$. Summation over repeated indices is implied and we introduced the strain tensor (relative to the rigid moiré)

$$u_{ij}^{(l)}=\frac{1}{2}\left(\frac{\partial u_{l,j}}{\partial r_{i}}+\frac{\partial u_{l,i}}{\partial r_{j}}\right),$$ (2)

for $i,j=x,y$. Here, the $\bm{u}_{l}$ are the in-plane displacement fields for the acoustic degrees of freedom. To keep the theory analytically tractable, we do not consider out-of-plane deformations at this stage, but we will include them later. Indeed, for heterobilayers, buckling can become important and the quadratic contribution from out-of-plane displacements to the strain tensor can reduce strain at the cost of bending energy [33], but may further reduce the elastic energy.

The moiré interface between the layers is specified by a rigid displacement gradient [20]

$$M=\left(1+\frac{\mathcal{E}}{2}\right)R(\theta/2)-\left(1-\frac{\mathcal{E}}{2}\right)R(-\theta/2),$$ (3)

where $R$ is a standard counterclockwise rotation matrix and $\mathcal{E}$ is a symmetric $2\times 2$ matrix that encodes both lattice mismatch and external heterostrain. In the long-wavelength limit where the modulation of the interlayer stacking varies slowly compared to the interatomic scale, the moiré lattice can be defined by

$$M\left(\bm{r}+\bm{L}\right)=M\bm{r}+\bm{a},$$ (4)

because a shift by a layer lattice vector $\bm{a}$ yields the same local stacking configuration. One thus finds

$$\bm{L}=M^{-1}\bm{a},\qquad\bm{g}=M^{\top}\bm{b}.$$ (5)

In this work, we only consider rigid matrices $M$ that are invertible. When one eigenvalue of $M$ vanishes, one of the moiré lattice vectors diverges, giving a one-dimensional moiré [20]. When both eigenvalues vanish, there is no periodic moiré structure. For the case of lattice mismatch (or biaxial strain) and twist, we have

$$L=|\bm{L}_{1}|=|\bm{L}_{2}|=\frac{a}{\sqrt{\epsilon^{2}\cos^{2}\tfrac{\theta}{2}+4\sin^{2}\tfrac{\theta}{2}}}.$$ (6)

We further consider adiabatic relaxation such that the displacement fields have the same symmetry as the original rigid moiré. Hence, we do not consider broken-symmetry configurations as in Ref. [53], noting that out-of-plane buckling could influence the stability of these states. Hence, for adiabatic relaxation the translational symmetry of the rigid moré is conserved and the in-plane displacement fields can be written as

$$\bm{u}_{l}(\bm{r})=\sum_{\bm{b}}\bm{u}_{l,\bm{b}}e^{i\bm{b}\cdot M\bm{r}},$$ (7)

where the sum runs over reciprocal vectors of a fictitious layer that we introduced in Eq. (4). In particular, for heterobilayers with lattice mismatch, the lattice constant $a$ of this fictitious layer is given by the mean lattice constant of the physical layers:

$$a_{1,2}=\left(1\pm\frac{\epsilon}{2}\right)a,$$ (8)

with

$$a=\frac{a_{1}+a_{2}}{2},\qquad\epsilon=\frac{a_{1}-a_{2}}{a}.$$ (9)

For example, for graphene on hBN, taking $a_{\text{g}}\approx 2.466$ Å and $a_{\text{hBN}}\approx 2.513$ Å obtained from density-functional theory (see Appendix A for details) we obtain $a\approx 2.489$ Å and $\epsilon\approx 0.0187$. Here, we always take $a_{1}>a_{2}$ so that $\epsilon>0$. The moiré lattice constant for $\theta=0$ (aligned) is then given by $L\approx 14$ nm. The second contribution is the adhesion energy that captures the van der Waals bonding energy between the two layers:

$$H_{\text{adh}}=\int d^{2}\bm{r}\,\mathcal{V}\left[M\bm{r}+\bm{u}_{1}(\bm{r})-\bm{u}_{2}(\bm{r})\right],$$ (10)

with

$$\mathcal{V}(\bm{\phi})=\sum_{\bm{b}}\mathcal{V}_{\bm{b}}e^{i\bm{b}\cdot\bm{\phi}},$$ (11)

the stacking potential, shown in Fig. 1(c) for a $C_{3v}$ stacking energy landscape corresponding to graphene on hBN or heterobilayers of 2H and 1T TMDs, including aligned [27] and anti-aligned [37, 63] MoTe₂/WSe₂ and WS₂/WSe₂ [22, 18].

Next, we make use of the fact that a smooth vector field on a torus can always be expressed in a Helmholtz decomposition:

$$\bm{u}_{l,\bm{b}}=\begin{cases}\bm{0}&g=0,\\ \frac{a}{L}\frac{u_{l,\bm{b}}^{\parallel}\bm{g}+u_{l,\bm{b}}^{\perp}\hat{z}\times\bm{g}}{ig^{2}}&g\neq 0,\end{cases}$$ (12)

with $g=|\bm{g}|$ and where the harmonic part (i.e. the part that is both divergenceless and irrotational) has to be a constant because the torus has no boundaries. The constant is set to zero, as this only amounts to an overall shift of the moiré. Here, the longitudinal ($u^{\parallel}$) and transverse ($u^{\perp}$) components of the displacement field give the divergence and curl, respectively. In this way, the elastic energy per unit area becomes

$$h_{\text{elas}}=\frac{a^{2}}{2L^{2}}\sum_{l=1,2;\bm{b}}\left[\left(\lambda_{l}+2\mu_{l}\right)|u_{l,\bm{b}}^{\parallel}|^{2}+\mu_{l}|u_{l,\bm{b}}^{\perp}|^{2}\right].$$ (13)

$n,m$	$u_{n,m}^{\parallel}$	$u_{n,m}^{\perp}$	$h_{n,m}$
$1,0$	$\mathds{C}$	$0$	$\mathds{C}$
$2,1$	$\mathds{R}$	$i\mathds{R}$	$\mathds{R}$
$2,0$	$\mathds{C}$	$0$	$\mathds{C}$
$3,1$	$\mathds{C}$	$\mathds{C}$	$\mathds{C}$
$3,2$	$(u_{3,1}^{\parallel})^{*}$	$-(u_{3,1}^{\perp})^{*}$	$h_{3,1}^{*}$

Table 1: Constraints on the Fourier components of the adiabatic displacements fields of a heterobilayer moiré with $C_{3v}=\left<\mathcal{C}_{3z},\mathcal{M}_{x}\right>$ symmetry. This corresponds to perfectly aligned or anti-aligned heterobilayers, i.e. without any relative twist angle. The reciprocal stars are labeled by the tuple $(n,m)$ with representative $\bm{b}_{n,m}=n\bm{b}_{1}+m\bm{b}_{2}$ where $|\bm{b}_{1}|=4\pi/(\sqrt{3}a)$ and $\bm{b}_{2}=R(2\pi/3)\bm{b}_{1}$. Here $n=1,2,\ldots$ is the shell index and $m=0,\ldots,n-1$ are stars in a given shell.

II.1 Self-consistent solution

Minimizing the total energy with respect to the longitudinal and transverse components, $u_{l,-\bm{b}}^{\parallel}$ and $u_{l,-\bm{b}}^{\perp}$, respectively, yields [43]

$$u_{l,\bm{b}}^{\parallel} =\frac{L}{a}\frac{(-1)^{l+1}}{\lambda_{l}+2\mu_{l}}\frac{\bm{g}}{ig^{2}}\cdot\left(\frac{\partial\mathcal{V}}{\partial\bm{\phi}}\right)_{\bm{g}},$$ (14)

$$u_{l,\bm{b}}^{\perp} =\frac{L}{a}\frac{(-1)^{l+1}}{\mu_{l}}\frac{\hat{z}\times\bm{g}}{ig^{2}}\cdot\left(\frac{\partial\mathcal{V}}{\partial\bm{\phi}}\right)_{\bm{g}},$$ (15)

for $\bm{\phi}=M\bm{r}+\bm{u}_{1}-\bm{u}_{2}$ and $l=1,2$. We immediately notice

$$u_{2,\bm{b}}^{\parallel}=-\frac{\lambda_{1}+2\mu_{1}}{\lambda_{2}+2\mu_{2}}\,u_{1,\bm{b}}^{\parallel},\qquad u_{2,\bm{b}}^{\perp}=-\frac{\mu_{1}}{\mu_{2}}\,u_{1,\bm{b}}^{\perp},$$ (16)

such that $(\lambda_{1}+2\mu_{1})\nabla\cdot\bm{u}_{1}=-(\lambda_{2}+2\mu_{2})\nabla\cdot\bm{u}_{2}$ and $\mu_{1}\nabla\times\bm{u}_{1}=-\mu_{2}\nabla\times\bm{u}_{2}$ for free-standing heterobilayers, which reduces to the homobilayer result $\bm{u}_{1}=-\bm{u}_{2}$ [11, 21, 14] when $\lambda_{1}=\lambda_{2}$ and $\mu_{1}=\mu_{2}$.

Equations (14) and (15) yield four self-consistent equations for every reciprocal vector $\bm{b}$ that are solved numerically with the DIIS algorithm [49]. Generalization to multiple layers is straightforward as long as one can define a single moiré lattice with adhesion potentials $\mathcal{V}_{l,l+1}$ and stacking configurations $\bm{\phi}_{l,l+1}$.

II.2 One-shot result

To obtain approximate expressions for the relaxed configuration, we solve the self-consistent equations perturbatively in lowest order. This is equivalent to the one-shot result obtained by setting $\bm{u}_{l}=\bm{0}$ on the right-hand side of Eqs. (14) and (15). We find

$$\left.\left(\frac{\partial\mathcal{V}}{\partial\bm{\phi}}\right)_{\bm{g}}\right|_{\text{one shot}}=i\bm{b}\mathcal{V}_{\bm{b}},$$ (17)

with $\bm{g}=M^{\top}\bm{b}$. Here, the small parameter is the ratio of the characteristic adhesion and elastic energy. For example, $(\mathcal{V}_{\text{max}}-\mathcal{V}_{\text{min}})/[(\lambda+2\mu)\epsilon^{2}]$ for $\theta=0$. Hence, the one-shot solution is given by

$$u_{l,\bm{b}}^{\parallel} \simeq\frac{L}{a}\frac{(-1)^{l+1}}{\lambda_{l}+2\mu_{l}}\frac{\bm{b}\cdot M\bm{b}}{g^{2}}\,\mathcal{V}_{\bm{b}},$$ (18)

$$u_{l,\bm{b}}^{\perp} \simeq\frac{L}{a}\frac{(-1)^{l+1}}{\mu_{l}}\frac{\left(\hat{z}\times\bm{b}\right)\cdot M\bm{b}}{g^{2}}\,\mathcal{V}_{\bm{b}},$$ (19)

which is applicable to any moiré bilayer composed of layers with the same Bravais lattice type. For lattice mismatch $\epsilon$ and twist $\theta$, we obtain

$$u_{l,\bm{b}}^{\parallel} \simeq\frac{(-1)^{l+1}}{\lambda_{l}+2\mu_{l}}\frac{\epsilon\mathcal{V}_{\bm{b}}\cos\frac{\theta}{2}}{\left(\epsilon^{2}\cos^{2}\frac{\theta}{2}+4\sin^{2}\frac{\theta}{2}\right)^{3/2}},$$ (20)

$$u_{l,\bm{b}}^{\perp} \simeq\frac{(-1)^{l+1}}{\mu_{l}}\frac{2\mathcal{V}_{\bm{b}}\sin\frac{\theta}{2}}{\left(\epsilon^{2}\cos^{2}\frac{\theta}{2}+4\sin^{2}\frac{\theta}{2}\right)^{3/2}}.$$ (21)

We now take a $C_{3v}$ stacking energy:

$$\mathcal{V}(\bm{\phi})=2\mathcal{V}_{1}\sum_{n=1}^{3}\cos\left(\bm{b}_{n}\cdot\bm{\phi}+\varphi_{1}\right),$$ (22)

where we take $\mathcal{V}_{1}>0$ and which is a first-star approximation for graphene on hBN and 2H TMD heterobilayers with $\mathcal{V}_{\bm{b}_{1}}=\mathcal{V}_{\bm{b}_{2}}=\mathcal{V}_{\bm{b}_{3}}=\mathcal{V}_{1}e^{i\varphi_{1}}$. Here, we take reciprocal vectors $\bm{b}_{1}=4\pi\hat{y}/(\sqrt{3}a)$ and $\bm{b}_{n+1}=R(2\pi/3)\bm{b}_{n}$. For graphene on hBN, we use the parameters listed in Tables 2 and 3 obtained from density-functional theory calculations for untwisted bilayers (see Appendix A).

	$a$	$\lambda$	$\mu$	$\nu$
graphene	$2.466$	$26756$	$44702$	$0.230$
hBN	$2.513$	$22621$	$37791$	$0.230$
2H WS2	$3.188$	$30226$	$31405$	$0.323$
2H WSe2	$3.322$	$24587$	$30099$	$0.290$
2H MoTe2	$3.567$	$22130$	$24161$	$0.314$

Table 2: Lattice constant $a$ in Å, 2D Lamé parameters $\lambda$ and $\mu$ in units of meV per monolayer unit cell, and the 2D Poisson’s ratio $\nu=1/\left(1+2\mu/\lambda\right)$.

In general, the acoustic displacement fields of the adiabatic relaxed ground-state configuration are constrained by the symmetry of the rigid moiré [21]. This results in symmetry constraints for the Fourier components of each reciprocal star, given in Table 1. Here, a star is defined as a set of six reciprocal vectors closed under $\mathcal{C}_{6z}$. The self-consistent solution inherits these symmetries from the adhesion potential. In lowest order,

$$\bm{u}_{l}(\bm{r})\simeq\frac{\sqrt{3}a}{2\pi}\sum_{n=1}^{3} \left(u_{l,1}^{\parallel}\hat{g}_{n}+u_{l,1}^{\perp}\hat{z}\times\hat{g}_{n}\right)$$ (23)

$$\times\sin(\bm{g}_{n}\cdot\bm{r}+\varphi_{1}),$$

with the one-shot result

$$u_{l,m}^{\parallel} =\frac{(-1)^{l+1}}{\lambda_{l}+2\mu_{l}}\frac{\epsilon\mathcal{V}_{m}}{(\epsilon^{2}+\theta^{2})^{3/2}},$$ (24)

$$u_{l,m}^{\perp} =\frac{(-1)^{l+1}}{\mu_{l}}\frac{\theta\mathcal{V}_{m}}{(\epsilon^{2}+\theta^{2})^{3/2}},$$ (25)

for small $\theta$ and where $m$ labels the star. These expressions recover the results from Refs. [21, 12] for homobilayers with $\epsilon\rightarrow 0$ and finite $\theta$.

Figure 2: Acoustic displacement fields induced by atomic reconstruction in moiré TMD heterobilayers. Here, we show the longitudinal and transverse Fourier components of the first reciprocal star for the in-plane displacement fields $\bm{u}=\bm{u}_{1}-\bm{u}_{2}$ where the dots are numerical solutions of the self-consistent equations given in Eqs. 14 and 15, and the curves give the one-shot analytical result. Shown as a function of the twist angle for aligned MoTe₂/WSe₂ in panels (a) and (b), and for anti-aligned WSe₂/WS₂ in panels (d) and (e). We also show the transverse component versus $\epsilon$ for aligned MoTe₂/WSe₂ in panels (c) and for anti-aligned WSe₂/WS₂ in panels (f) for $\theta=0$. The center-of-mass fields are obtained from $U_{\bm{b}}^{\parallel}=A^{\parallel}u_{\bm{b}}^{\parallel}$ and $U_{\bm{b}}^{\perp}=A^{\perp}u_{\bm{b}}^{\perp}$. For aligned MoTe₂/WSe₂, $A^{\parallel}\approx 0.09$ and $A^{\perp}\approx 0.11$. For anti-aligned WSe₂/WS₂, $A^{\parallel}\approx 0.05$ and $A^{\perp}\approx 0.02$. Parameters are shown in Tables 2 and 3.

It is now instructive to write the original layer-resolved fields in terms of the layer center-of-mass ($\bm{U}$) and relative displacements ($\bm{u}$). This is motivated by the fact that the latter is dominant, while the former is expected to be much smaller. We define these auxiliary fields as

$$\bm{u}_{1,2}=\frac{\bm{U}\pm\bm{u}}{2}.$$ (26)

Moreover,

$$U_{\bm{b}}^{\parallel} =\frac{\lambda_{2}+2\mu_{2}-\lambda_{1}-2\mu_{1}}{\lambda_{2}+2\mu_{2}+\lambda_{1}+2\mu_{1}}\,u_{\bm{b}}^{\parallel}\equiv A^{\parallel}u_{\bm{b}}^{\parallel},$$ (27)

$$U_{\bm{b}}^{\perp} =\frac{\mu_{2}-\mu_{1}}{\mu_{2}+\mu_{1}}\,u_{\bm{b}}^{\perp}\equiv A^{\perp}u_{\bm{b}}^{\perp},$$ (28)

such that we only specify $u_{\bm{b}}^{\parallel}$ and $u_{\bm{b}}^{\perp}$. This works because the longitudinal and transverse components are decoupled in the elastic energy. It further implies $\bm{U}=A^{\parallel}\nabla\phi+A^{\perp}\nabla\times\bm{w}$, with $\bm{u}=\nabla\phi+\nabla\times\bm{w}$.

In Appendix B we also show the elastic energy in terms of these fields. We compare these perturbative results to the numerical results in Fig. 1. We find excellent agreement between the analytical theory and the numerics for the case of graphene on hBN. We see that the longitudinal components (related to local compression and dilation) are largest in magnitude for $\theta_{\parallel}=0$ while the transverse components (related to twirling or local rotations) are maximal for $\theta_{\perp}=\epsilon/\sqrt{2}$. For graphene on hBN, we find $\theta_{\perp}\approx 0.76^{\circ}$. In addition, we also apply the analytical theory to several 2H TMD heterobilayers in Fig. 2. Now the stacking energy contains non-negligible contributions from the second and third star. For example, up to the second star, $\bm{u}_{l}\simeq\bm{u}_{l,1}+\bm{u}_{l,2}$ with

$$\bm{u}_{l,2}(\bm{r})=\frac{a}{2\pi}\sum_{n=1}^{3}\left(u_{l,2}^{\parallel}\hat{g}_{n}^{\prime}+u_{l,2}^{\perp}\hat{z}\times\hat{g}_{n}^{\prime}\right)\sin(\bm{g}_{n}^{\prime}\cdot\bm{r}),$$ (29)

where $\bm{g}_{1}^{\prime}=2\bm{g}_{1}+\bm{g}_{2}$ and $\bm{g}_{2,3}^{\prime}$ related by $120^{\circ}$ rotations, are reciprocal vectors of the second star. Thus, the relaxed stacking configuration in experimentally relevant heterobilayers (aligned, anti-aligned, or twisted) can be accurately captured using a perturbative solution. This is because the moiré length in these systems is bounded from above by the lattice mismatch.

layer 1/layer 2	$\epsilon$	$a$	$\mathcal{V}_{1}$	$\mathcal{V}_{2}$	$\mathcal{V}_{3}$	$\mathcal{V}_{4,5}$	$\varphi_{1}$	$\varphi_{3}$	$\varphi_{4}=-\varphi_{5}$
hBN/graphene	$0.0187$	$2.489$	$1.421$	$-0.07514$	$0.01734$	$0.02080$	$0.8867$	$-1.491$	$2.535$
P WSe2/WS2	$0.0412$	$3.255$	$11.09$	$-2.090$	$0.7501$	$0.2748$	$0.06173$	$-2.958$	$-0.07944$
P MoTe2/WSe2	$0.0711$	$3.445$	$13.61$	$-2.675$	$0.9122$	$0.3566$	$0.04351$	$-2.985$	$-0.07550$
AP WSe2/WS2	$0.0412$	$3.255$	$9.665$	$-1.608$	$0.6338$	$0.1922$	$-2.301$	$-1.375$	$2.144$
AP MoTe2/WSe2	$0.0711$	$3.445$	$12.04$	$-2.079$	$0.8117$	$0.2598$	$2.354$	$1.480$	$-2.161$

Table 3: Lattice mismatch $\epsilon$, average lattice constant $a$ in units of Å, and Fourier coefficients of the adhesion potential in units of meV per monolayer unit cell (taking the average lattice constant) for the heterobilayers considered in this work. Angles are given in radians. We always take $a_{1}>a_{2}$ such that $\epsilon>0$ and we consider both aligned (P) and anti-aligned (AP) stacking.

III Buckling

We now consider the out-of-plane displacement fields $h_{l}$ ($l=1,2$). Here, we distinguish between breathing $h_{1}-h_{2}$ and buckling $h_{1}+h_{2}$. Buckling can be driven by in-plane compression to reduce strain [5, 8, 13, 40, 58]. Similarly, substrate engineering of out-of-plane deformations in graphene gives rise to local in-plane compressive deformations [47, 15, 16]. In free-standing twisted homobilayers, compression-driven buckling is absent as $\nabla\cdot\bm{u}_{l}\approx 0$ [21, 14]. Instead, homobilayers may exhibit shear-driven buckling with spontaneously broken moiré translational symmetry [58]. Here, we consider buckling in heterobilayers without symmetry breaking. We thus expect that the buckling amplitude is largest for aligned (or anti-aligned) layers, for which the longitudinal components (and thus local compression) are dominant. With increasing twist angle, the latter and thus buckling are suppressed until a buckling transition occurs at a critical twist angle.

Here, we present a simplified theory of buckling by considering a single layer in the presence of the relaxation-induced in-plane displacement fields that we obtained in the previous section. We then minimize the elastic energy of this layer with respect to the height profile. This theory captures the qualitative features of the buckling transition, while details of the buckling profile and the critical angle are expected to differ in a more accurate theory. Such a theory should include the coupling between buckling and breathing modes through the interlayer interaction, i.e. one needs to include the dependence of the interlayer separation in the Fourier components of the adhesion potential [19].

III.1 Symmetry ansatz

The details of the buckling theory are presented in Appendix C. To simplify the calculation, we use an ansatz for the height profile with $C_{3v}$ symmetry:

$$h(\bm{r}) =2\sum_{n=1}^{3}\Big[h_{1}\cos(\bm{g}_{n}\cdot\bm{r}+\gamma_{1})+h_{2}\cos(\bm{g}_{n}^{\prime}\cdot\bm{r})$$ (30)

$$+h_{3}\cos(2\bm{g}_{n}\cdot\bm{r}+\gamma_{3})+\left(\text{higher stars}\right)\Big],$$

with amplitudes $h_{n}$ and phases $\varphi_{n}$. The elastic energy is then computed analytically and numerically minimized with respect to the amplitudes and phases. We have to go beyond the first star because the height profile enters nonlinear in the strain. This is similar to the reverse problem, where for a fixed first-star buckling profile, the elastic ground state contains up to three stars of the in-plane field [16]. We find that the in-plane rotational components ($u^{\perp}$) do not couple to the height profile at this order of approximation. This is enforced by symmetry as detailed in Appendix C. Moreover, the transverse in-plane components are also decoupled from the longitudinal ones ($u^{\parallel}$) [see Eq. (13)], such that the $u^{\perp}$ do not affect buckling.

The in-plane displacement field due to structural relaxation in moiré heterobilayers gives rise to both regions of in-plane tensile and compressive strain. Because the volumetric part of the strain tensor from the height modulation is given by $(1/2)(\nabla h)^{2}\geq 0$, regions of in-plane tensile strain favor a constant height profile, while compressive regions favor a varying $h(\bm{r})$ to lower the total strain. These gains compete with the bending energy $(\kappa/2)(\nabla^{2}h)^{2}$ which disfavors out-of-plane deformations.

Figure 3: Buckling transition in moiré heterobilayers. Shown for graphene on hBN for $\kappa=1$ eV and the parameters in Tables 2 and 3. (a) Fourier amplitudes $h_{n}$ as a function of twist angle for a $C_{3v}$ buckling profile with up to 7 stars. (b) Relative energy gain due to buckling versus twist angle. (c) Phases of the Fourier components. The 2nd and 7th star are real by symmetry and not shown. (d) Divergence of the in-plane displacement field for graphene aligned to hBN ($\theta=0$) using the one-shot analytical result. (e) Buckling profile corresponding to (d) using the results from (a) and (c).

III.2 Buckling transition

We start by considering only the first star in the ansatz from Eq. (30). While, we numerically find that we need up to seven stars for convergence, we can analytically demonstrate that a buckling transition can occur in principle. In particular, we find

$$h_{1}^{2}=\frac{(-1)^{n}3aLu_{1}\left(\mu-\lambda\right)/\pi^{2}-16\kappa}{72\left(\lambda+2\mu\right)},$$ (31)

with $\gamma_{1}=-\left(\varphi_{1}+n\pi\right)/2$ and $n=0,1$. Here, $a=L\sqrt{\epsilon^{2}+\theta^{2}}$ and $u_{1}$ is the longitudinal first-star component. Since $h_{1}$ is real, this requires

$$(-1)^{n}u_{1}>\frac{\kappa}{\mu-\lambda}\frac{16\pi^{2}}{3aL},$$ (32)

or $h_{1}=0$ otherwise. Because we consider a single layer, the solution is degenerate under $h_{1}\rightarrow-h_{1}$. In a moiré bilayer, there is a preferential orientation due to coupling between the buckling and breathing modes. Using our previous result: $u_{1}=(-1)^{l+1}\mathcal{V}_{1}\epsilon/[(\lambda+2\mu)(\epsilon^{2}+\theta^{2})^{3/2}]$ for layer $l=1,2$, together with Eq. (31), we find a nontrivial buckling solution for

$$\frac{\epsilon}{\left(\epsilon^{2}+\theta^{2}\right)^{2}}>\frac{32\pi^{2}}{1-3\nu}\frac{\kappa}{3a^{2}\mathcal{V}_{1}},$$ (33)

where $n=0$ for the first layer and $n=1$ for the second layer with $a_{1}>a_{2}$. Here, we defined $\nu=1/\left(1+2\mu/\lambda\right)$ the 2D Poisson’s ratio and $a$ the average lattice constant of the two layers. Thus, the first-star approximation predicts a buckling transition for $\epsilon\sim\kappa a^{2}/(\mathcal{V}_{1}L^{4})$ or below a critical twist angle

$$\left.\theta_{c}^{2}\right|_{\text{first star}}=\frac{a}{4\pi}\sqrt{\frac{3\epsilon(1-3\nu)\mathcal{V}_{1}}{2\kappa}}-\epsilon^{2},$$ (34)

determined by the ratio $\epsilon(1-3\nu)a^{2}\mathcal{V}_{1}/\kappa$. This is to be expected, as the energy cost for buckling is proportional to $\kappa$ while the gain is proportional to $|u_{1}|\propto\mathcal{V}_{1}$. While this result provides good intuition, an accurate determination of $\theta_{c}$ is done numerically. We take a $C_{3v}$ buckling profile up to the seventh star, and numerically compute the coefficients $h_{\bm{g}}$ that minimize the elastic energy for a fixed in-plane displacement field from Eq. (23). Here, we only consider the layer with the smaller lattice constant. Note that the in-plane displacement field for the other layer has the opposite sign, and the resulting buckling profiles are not equivalent. Using $\kappa=1$ eV, we find a buckling transition for graphene on hBN below critical twist angle $\theta_{c}\approx 1.1^{\circ}$. This is illustrated in Fig. 3(a) and (b) where we show the $h_{n}$ for each star and the relative energy gain, respectively, as a function of twist angle. We also show the corresponding phases $\gamma_{n}$ in Fig. 3(c). The in-plane volumetric strain and buckling profile for $\theta=0$ are shown in Fig. 3(d) and (e), respectively. We see that compressive in-plane strain ($\nabla\cdot\bm{u}<0$) favors buckling, while tensile strain ($\nabla\cdot\bm{u}>0$) disfavors buckling, leading to a bucket-shaped buckling profile. For aligned (anti-aligned) WSe₂/WS₂ we find $\theta_{c}\approx 2.2^{\circ}(2^{\circ})$, while we find no buckling for MoTe₂/WSe₂ heterobilayers.

IV Conclusion

In this work, we developed an analytical theory of lattice relaxation in moiré heterobilayers. We modeled the intralayer elastic properties with continuum elasticity and the interlayer energy with the local stacking approximation using parameters obtained from density functional theory for untwisted bilayers. We then derived the self-consistent equations for the Fourier components of the acoustic in-plane displacement fields and obtained the one-shot perturbative result. We subsequently applied this theory to graphene on hBN and representative 2H transition metal dichalcogenide heterobilayers, and found good agreement between the analytical results and the full numerical solutions over experimentally relevant parameters. These results show that the relaxed stacking configuration in realistic heterobilayer moiré systems can be captured accurately within a simple analytical framework, even in the presence of both lattice mismatch and twist. In addition, using our analytical result for the in-plane displacement field, we demonstrated that moiré heterobilayers exhibit a buckling instability near alignment, driven by compressive strain due to in-plane relaxation. More broadly, our results provide both a practical framework for treating lattice relaxation in realistic moiré heterostructures and a foundation for integrating structural reconstruction into calculations of their electronic and optical properties.

Appendix A First-principles calculations

First-principles density functional theory (DFT) calculations were performed to simulate graphene, hBN, and TMDs using the abinit [23, 24, 57] code. PAW pseudopotentials were used [55], obtained from Pseudo-Dojo [56], and the PBE functional was used [46] to treat the exchange-correlation interactions. abinit employs a plane-wave basis set, which was determined using a kinetic energy cutoff of 1200 eV, and a cutoff of 1500 eV was used for PAW double grid. A Monkhorst-Pack $k$-point grid [41] of $12\times 12\times 1$ was used to sample the Brillouin zone for hBN and the TMDs, and a sampling of $21\times 21\times 1$ was used for graphene, which requires a finer sampling as it is not gapped.

Geometry relaxations were performed for the graphene, hBN, WS₂, WSe₂ and MoTe₂ monolayers to obtain the equilibrium lattice constants and elastic properties. The bulk and shear moduli were calculated by making small deviations of the lattice constant and the angle between the lattice vectors, respectively, and taking the second derivatives of the change in energy. The equilibrium lattice constants, Lamé factors, and Poisson ratios are given in Table 2.

Commensurate hBN/graphene, WS₂/WSe₂, and MoTe₂/WSe₂ heterostructures were then constructued, using the average of the lattice constants of the individual monolayers in each case. To sample the different local environments, a grid containing $31$ stacking configurations not related by symmetry on a triangular grid was used, which is compatible with the $C_{3v}$ symmetry of the adhesion potential of the heterostructures. For each relative stacking, a geometry relaxation was performed to obtain the equilibrium layer separation, while keeping the in-plane lattice vectors and atomic positions fixed, using a force tolerance of 1 meV/Å. The energy as a function of stacking was then linearly interpolated and numerically Fourier transformed to obtain the coefficients

$$\mathcal{V}_{\bm{b}}=\frac{1}{A_{c}}\int_{\text{cell}}d^{2}\bm{r}\,\mathcal{V}(\bm{r})e^{-i\bm{b}\cdot\bm{r}},$$ (35)

given in Table 2.

Appendix B Self-consistent equations

In this section, we give the full derivation of the self-consistent equations obtained from continuum elasticity, for the in-plane displacement field of a general moiré bilayer. To this end, we first write down the Fourier transform of the strain tensor:

$$u_{ij}=\frac{i}{2}\sum_{\bm{b}}\left(g_{i}u_{l,\bm{b},j}+g_{j}u_{l,\bm{b},i}\right)e^{i\bm{g}\cdot\bm{r}},$$ (36)

where $\bm{g}=M^{\top}\bm{b}$ so we can freely move between the atomic lattice scale $\bm{b}$ and the moiré scale $\bm{g}$. Plugging this into the expression for the elastic energy and integrating, using

$$\frac{1}{A}\int d^{2}\bm{r}\,f(\bm{r})e^{-i\bm{g}\cdot\bm{r}}=f_{\bm{g}},$$ (37)

with $A$ the total area and $f(\bm{r})=f(\bm{r}+\bm{L})$ a periodic function, we find

$$h_{\text{elas}} =\frac{1}{2}\sum_{l=1,2}\sum_{\bm{b}}\left[\lambda_{l}\left|\bm{g}\cdot\bm{u}_{l,\bm{b}}\right|^{2}+\frac{\mu_{l}}{2}\left|g_{i}u_{l,\bm{b},j}+g_{j}u_{l,\bm{b},i}\right|^{2}\right]$$ (38)

$$=\frac{1}{2}\sum_{l=1,1}\sum_{\bm{b}}\left[\left(\lambda_{l}+2\mu_{l}\right)\left|\bm{g}\cdot\bm{u}_{l,\bm{b}}\right|^{2}+\mu_{l}\left|\left(\hat{z}\times\bm{g}\right)\cdot\bm{u}_{l,\bm{b}}\right|^{2}\right]$$ (39)

$$=\frac{a^{2}}{2L^{2}}\sum_{l=1,1}\sum_{\bm{b}}\left[\left(\lambda_{l}+2\mu_{l}\right)|u_{l,\bm{b}}^{\parallel}|^{2}+\mu_{l}|u_{l,\bm{b}}^{\perp}|^{2}\right],$$ (40)

where we used that a smooth vector field on a torus can always be expressed in a Helmholtz decomposition:

$$\bm{u}_{l,\bm{b}}=\begin{cases}\bm{0}&g=0,\\ \frac{a}{L}\frac{u_{l,\bm{b}}^{\parallel}\bm{g}+u_{l,\bm{b}}^{\perp}\hat{z}\times\bm{g}}{ig^{2}}&g\neq 0,\end{cases}$$ (41)

with $g=|\bm{g}|$ and where the harmonic part (i.e. a gradient $\nabla f$ with $\nabla^{2}f=0$) has to be a constant because the torus has no boundaries. This constant is set to zero, as this only amounts to an overall shift of the moiré lattice. The longitudinal and transverse components of the displacement field give the divergence and curl, respectively,

$$\nabla\cdot\bm{u}_{l} =\frac{a}{L}\sum_{\bm{h}}u_{l,\bm{g}}^{\parallel}e^{i\bm{g}\cdot\bm{r}},$$ (42)

$$\nabla\times\bm{u}_{l} =\frac{a}{L}\sum_{\bm{h}}u_{l,\bm{g}}^{\perp}e^{i\bm{g}\cdot\bm{r}}.$$ (43)

We also find

$$h_{\text{adh}} =\frac{1}{A_{m}}\int_{\text{moir\'{e} cell}}d^{2}\bm{r}\,\mathcal{V}\left[M\bm{r}+\bm{u}_{1}(\bm{r})-\bm{u}_{2}(\bm{r})\right]$$ (44)

$$=\frac{1}{A_{m}}\sum_{\bm{b}}\mathcal{V}_{\bm{b}}\int_{\text{moir\'{e} cell}}d^{2}\bm{r}\,\exp\left\{i\bm{b}\cdot\left[M\bm{r}+\bm{u}_{1}(\bm{r})-\bm{u}_{2}(\bm{r})\right]\right\},$$ (45)

with $A_{m}=|\bm{L}_{1}\times\bm{L}_{2}|=(L/a)^{2}|\bm{a}_{1}\times\bm{a}_{2}|$ the moiré unit cell area. The ground-state is obtained by minimizing the total elastic and adhesion energy with respect to $\{u_{l,\bm{b}}^{\parallel},u_{l,\bm{b}}^{\perp}\}$ for $l=1,2$. We find

$$\frac{\partial h_{\text{elas}}}{\partial u_{l,-\bm{b}}^{\parallel}}=\frac{a^{2}}{L^{2}}\left(\lambda_{l}+2\mu_{l}\right)u_{l,\bm{b}}^{\parallel},\qquad\frac{\partial h_{\text{elas}}}{\partial u_{l,-\bm{b}}^{\perp}}=\frac{a^{2}}{L^{2}}\mu_{l}u_{l,\bm{b}}^{\perp},$$ (46)

and

$$\frac{\partial h_{\text{adh}}}{\partial u_{l,-\bm{b}}^{\parallel}} =\frac{1}{A_{m}}\int_{\text{moir\'{e} cell}}d^{2}\bm{r}\,\frac{\partial\bm{u}_{l}}{\partial u_{l,-\bm{b}}^{\parallel}}\cdot\frac{\partial\mathcal{V}}{\partial\bm{u}_{l}}\left[M\bm{r}+\bm{u}_{1}(\bm{r})-\bm{u}_{2}(\bm{r})\right]$$ (47)

$$=-\frac{a}{L}\frac{\bm{g}}{ig^{2}}\cdot\frac{1}{A_{m}}\int_{\text{moir\'{e} cell}}d^{2}\bm{r}\,e^{-i\bm{g}\cdot\bm{r}}\frac{\partial\mathcal{V}}{\partial\bm{u}_{l}}\left[M\bm{r}+\bm{u}_{1}(\bm{r})-\bm{u}_{2}(\bm{r})\right]$$ (48)

$$=(-1)^{l}\frac{a}{L}\frac{\bm{g}}{ig^{2}}\cdot\left(\frac{\partial\mathcal{V}}{\partial\bm{\phi}}\right)_{\bm{g}},$$ (49)

$$\frac{\partial h_{\text{adh}}}{\partial u_{l,-\bm{b}}^{\perp}} =(-1)^{l}\frac{a}{L}\frac{\hat{z}\times\bm{g}}{ig^{2}}\cdot\left(\frac{\partial\mathcal{V}}{\partial\bm{\phi}}\right)_{\bm{g}},$$ (50)

with $\bm{\phi}=M\bm{r}+\bm{u}_{1}-\bm{u}_{2}$ the local stacking. Using these equations, we then obtain the self-consistent equations shown in Eqs. (14) and (15) of the main text. In Fig. 4, we show the local stacking energy for the relaxed configuration for graphene on hBN. Numerical results for the Fourier components for graphene on hBN are shown in Fig. 5 and for 2H TMD heterobilayers in in Fig. 6.

Figure 4: Stacking energy landscape (up to an overall constant) in units of meV per monolayer unit cell (using the average lattice constant). Shown for the relaxed moiré configuration $\bm{\phi}(\bm{r})=\epsilon\bm{r}+\bm{u}(\bm{r})$ for graphene aligned on hBN. Here, we used the analytical result for the relaxed in-plane acoustic displacement fields. Parameters used are listed in Tables 2 and 3 of the main text.

Figure 5: Fourier components of the relative displacement field $\mathbf{u}=\mathbf{u}_{1}-\mathbf{u}_{2}$ for graphene on hBN at $\theta=0^{\circ}$ and $\theta=0.8^{\circ}$, with $\epsilon=0.0187$. The transverse and longitudinal components, $u^{\perp}_{\mathbf{b}}$ and $u^{\parallel}_{\mathbf{b}}$, are shown up to the sixth shell. Dot size denotes $|u_{\mathbf{b}}|$, and color denotes the phase $\arg(u_{\mathbf{b}})$. The indicated scale factors are applied to the dot sizes for visual clarity.

Figure 6: Fourier components of the relative displacement field $\mathbf{u}=\mathbf{u}_{1}-\mathbf{u}_{2}$ for TMD heterobilayers. The top row shows aligned $\mathrm{MoTe}_{2}/\mathrm{WSe}_{2}$ with $\epsilon=0.0711$, while the bottom row shows anti-aligned $\mathrm{WSe}_{2}/\mathrm{WS}_{2}$ with $\epsilon=0.0412$. For each system, $u^{\perp}_{\mathbf{b}}$ and $u^{\parallel}_{\mathbf{b}}$ are shown at $\theta=0^{\circ}$ and at the representative finite twist angle labeled above each panel. Results are shown up to the sixth shell. Dot size denotes $|u_{\mathbf{b}}|$, color denotes $\arg(u_{\mathbf{b}})$, and the indicated scale factors are used for visual clarity.

B.1 Perturbation theory

It is convenient to write the original layer-resolved fields in terms of the layer center-of-mass ($\bm{U}$) and relative displacements ($\bm{u}$). This is motivated by the fact that the latter is dominant, while the former is expected to be much smaller. We define these auxiliary fields as

$$\bm{u}_{1,2}=\frac{\bm{U}\pm\bm{u}}{2}.$$ (51)

We can rewrite the elastic energy in terms of these auxiliary fields. Because we only consider in-plane displacements, the strain tensor is linear,

$$u_{ij}^{(1,2)}=\frac{U_{ij}\pm u_{ij}}{2},$$ (52)

and

$$u_{ii}^{(1,2)}u_{ii}^{(1,2)} =\frac{1}{4}\left(U_{ii}\pm u_{ii}\right)\left(U_{ii}\pm u_{ii}\right)$$ (53)

$$=\frac{1}{4}\left(U_{ii}U_{ii}+u_{ii}u_{ii}\pm 2u_{ii}U_{ii}\right),$$ (54)

and similar for the second term in the elastic energy. We therefore obtain

$$H_{\text{elas}} =\frac{1}{4}\int d^{2}\bm{r}\,\Big(\lambda U_{ii}U_{ii}+2\mu U_{ij}U_{ji}$$ (55)

$$+\lambda u_{ii}u_{ii}+2\mu u_{ij}u_{ji}+\bar{\lambda}u_{ii}U_{ii}+2\bar{\mu}u_{ij}U_{ji}\Big),$$

with effective Lamé parameters

$$\lambda =\frac{\lambda_{1}+\lambda_{2}}{2},\qquad \mu =\frac{\mu_{1}+\mu_{2}}{2},$$ (56)

$$\bar{\lambda} =\lambda_{1}-\lambda_{2},\qquad \bar{\mu} =\mu_{1}-\mu_{2}.$$ (57)

To obtain approximate expressions for the relaxed configuration, we solve the self-consistent equations perturbatively in lowest order. Here, the small parameter is given by the ratio of the characteristic adhesion and elastic energy. We now make use of the fact that the acoustic displacement fields of the adiabatic relaxed ground-state configuration are constrained by the symmetry of the rigid moiré. These symmetry constraints (see Table 1) give rise to the following first-star ansatze:

$$\bm{u}(\bm{r}) =\frac{\sqrt{3}a}{2\pi}\sum_{n=1}^{3}\left[u_{1}^{\parallel}\hat{g}_{n}\sin(\bm{g}_{n}\cdot\bm{r}+\alpha_{1}^{\parallel})+u_{1}^{\perp}\hat{z}\times\hat{g}_{n}\sin(\bm{g}_{n}\cdot\bm{r}+\alpha_{1}^{\perp})\right],$$ (58)

$$\bm{U}(\bm{r}) =\frac{\sqrt{3}a}{2\pi}\sum_{n=1}^{3}\left[U_{1}^{\parallel}\hat{g}_{n}\sin(\bm{g}_{n}\cdot\bm{r}+\beta_{1}^{\parallel})+U_{1}^{\perp}\hat{z}\times\hat{g}_{n}\sin(\bm{g}_{n}\cdot\bm{r}+\beta_{1}^{\perp})\right],$$ (59)

where the subscript refers the first star of reciprocal vectors instead of the layer index. For example, we defined $u_{1,0}^{\parallel}=u_{1}^{\parallel}e^{i\alpha_{1}^{\parallel}}$, $u_{1,0}^{\perp}=u_{1}^{\perp}e^{i\alpha_{1}^{\perp}}$, and similarly for the center-of-mass displacement. If we plug these expression in the elastic and adhesion energy, we find

$$h_{\text{elas}} =\frac{a^{2}}{L^{2}}\Big\{\frac{3\mu}{2}\left[(U_{1}^{\perp})^{2}+(u_{1}^{\perp})^{2}\right]+\frac{3\bar{\mu}}{2}u_{1}^{\perp}U_{1}^{\perp}\cos(\alpha_{1}^{\perp}-\beta_{1}^{\perp})$$ (60)

$$+\frac{3\left(\lambda+2\mu\right)}{2}\left[(U_{1}^{\parallel})^{2}+(u_{1}^{\parallel})^{2}\right]+\frac{3\left(\bar{\lambda}+2\bar{\mu}\right)}{2}u_{1}^{\parallel}U_{1}^{\parallel}\cos(\alpha_{1}^{\parallel}-\beta_{1}^{\parallel})\Big\},$$ (61)

$$h_{\text{adh}} =-\frac{6\mathcal{V}_{1}L}{a}\left[u_{1}^{\parallel}\epsilon\cos(\alpha_{1}^{\parallel}-\varphi_{1})+u_{1}^{\perp}\theta\cos(\alpha_{1}^{\perp}-\varphi_{1})\right],$$ (62)

for small twist angles $M\approx\epsilon\sigma_{0}-i\theta\sigma_{y}$. Minimizing the total energy yields $\alpha_{1}^{\parallel}=\alpha_{1}^{\perp}=\beta_{1}^{\parallel}=\beta_{1}^{\perp}=\varphi_{1}$ and

$$u_{1}^{\parallel} =\left(\frac{1}{\lambda_{1}+2\mu_{1}}+\frac{1}{\lambda_{2}+2\mu_{2}}\right)\frac{\mathcal{V}_{1}\epsilon}{(\epsilon^{2}+\theta^{2})^{3/2}},\qquad U_{1}^{\parallel} =\left(\frac{1}{\lambda_{1}+2\mu_{1}}-\frac{1}{\lambda_{2}+2\mu_{2}}\right)\frac{\mathcal{V}_{1}\epsilon}{(\epsilon^{2}+\theta^{2})^{3/2}},$$ (63)

$$u_{1}^{\perp} =\left(\frac{1}{\mu_{1}}+\frac{1}{\mu_{2}}\right)\frac{\mathcal{V}_{1}\theta}{(\epsilon^{2}+\theta^{2})^{3/2}},\qquad U_{1}^{\perp} =\left(\frac{1}{\mu_{1}}-\frac{1}{\mu_{2}}\right)\frac{\mathcal{V}_{1}\theta}{(\epsilon^{2}+\theta^{2})^{3/2}},$$ (64)

which recovers the results from Refs. [21, 12] for homobilayers with $\epsilon\rightarrow 0$ and finite $\theta$.

Appendix C Buckling

C.1 Elastic energy

Including the out-of-plane bending term, the elastic energy of a single layer can be written as [44]

$$H_{\text{elas}}=\frac{1}{2}\int d^{2}\bm{r}\left[\lambda u_{ii}u_{ii}+2\mu u_{ij}u_{ji}+\kappa\left(\nabla^{2}h\right)^{2}\right],$$ (65)

with 2D Lamé parameters $\lambda$ and $\mu$, and where $\kappa$ is the out-of-plane bending rigidity. The strain tensor now also includes the quadratic out-of-plane part:

$$u_{ij}=\frac{1}{2}\left(\frac{\partial u_{j}}{\partial r_{j}}+\frac{\partial u_{i}}{\partial r_{j}}+\frac{\partial h}{\partial r_{i}}\frac{\partial h}{\partial r_{j}}\right).$$ (66)

Using the Helmholtz decomposition for a periodic in-plane displacement field (now without the prefactor $a/L$), we find that the elastic energy can be written in Fourier space as

$$\frac{1}{2}\sum_{\bm{g}}\left[\left(\lambda+2\mu\right)\left|u_{\bm{g}}^{\parallel}+\frac{f_{xx\bm{g}}+f_{yy\bm{g}}}{2}\right|^{2}+\mu\left(\frac{|u_{\bm{g}}^{\perp}|^{2}}{2}-\mathcal{G}_{-\bm{g}}u_{\bm{g}}^{\perp}-\mathcal{F}_{-\bm{g}}u_{\bm{g}}^{\parallel}+\frac{|f_{xy\bm{g}}|^{2}-f_{xx\bm{g}}f_{yy,-\bm{g}}}{2}+\mathrm{c.c.}\right)+\kappa|h_{\bm{g}}|^{2}g^{4}\right],$$ (67)

with $\mathcal{F}_{\bm{g}}=\left(g_{x}^{2}f_{yy\bm{g}}-2g_{x}g_{y}f_{xy\bm{g}}+g_{y}^{2}f_{xx\bm{g}}\right)/g^{2}$ which transforms as a scalar field, and is related to the Hessian determinant of the surface in Monge gauge. Here, we also defined the symmetric tensor $f_{ij}(\bm{r})=(\partial_{i}h)(\partial_{j}h)$ which is the metric tensor up to the identity matrix, with

$$f_{ij\bm{g}}=-\frac{1}{2}\sum_{\bm{g}^{\prime}}h_{\bm{g}^{\prime}}h_{\bm{g}-{\bm{g}}^{\prime}}\left(g_{i}g_{j}^{\prime}+g_{j}g_{i}^{\prime}-2g_{i}^{\prime}g_{j}^{\prime}\right).$$ (68)

Here, we also defined

$$\mathcal{G}_{\bm{g}}=\frac{\left(f_{xx\bm{g}}-f_{yy\bm{g}}\right)g_{x}g_{y}+f_{xy\bm{g}}\left(g_{y}^{2}-g_{x}^{2}\right)}{g^{2}}.$$ (69)

As a check, we recover the results of Refs. [25, 26, 59, 47, 16] if we minimize with respect to $u_{-\bm{g}}^{\parallel}$ for fixed $h_{\bm{g}}$. This yields

$$u_{\bm{g}}^{\parallel} =\frac{\mu\mathcal{F}_{\bm{g}}}{\lambda+2\mu}-\frac{f_{xx\bm{g}}+f_{yy\bm{g}}}{2},$$ (70)

$$u_{\bm{g}}^{\perp} =\mathcal{G}_{\bm{g}},$$ (71)

which agrees with Ref. [16]. Using the definition of $f_{ij\bm{g}}$, the right-hand side of Eq. (69) can be written as

$$\mathcal{G}_{\bm{g}}=\frac{1}{g^{2}}\sum_{\bm{g}^{\prime}}h_{\bm{g}^{\prime}}h_{\bm{g}-\bm{g}^{\prime}}\left(\bm{g}^{\prime}\cdot\bm{g}\right)\left(\bm{g}^{\prime}\times\bm{g}\right),$$ (72)

where we define the cross product between two in-plane vectors as a pseudoscalar. Here, we used

$$\sum_{\bm{g}^{\prime}}h_{\bm{g}^{\prime}}h_{\bm{g}-\bm{g}^{\prime}}\left(\bm{g}^{\prime}\times\bm{g}\right)=0.$$ (73)

Similarly, we find

$$\mathcal{F}_{\bm{g}}=\frac{1}{g^{2}}\sum_{\bm{g}^{\prime}}h_{\bm{g}^{\prime}}h_{\bm{g}-\bm{g}^{\prime}}\left(\bm{g}^{\prime}\times\bm{g}\right)^{2},$$ (74)

which is related to the Hessian determinant:

$$(\partial_{x}^{2}h)(\partial_{y}^{2}h)-(\partial_{x}\partial_{y}h)^{2}=\frac{1}{2}\sum_{\bm{g}}g^{2}\mathcal{F}_{\bm{g}}.$$ (75)

Unlike $\mathcal{F}_{\bm{g}}$ we see that $\mathcal{G}_{\bm{g}}$ transforms as a pseudoscalar under a symmetry $\mathcal{S}$ of the out-of-plane displacement field, i.e., for $h_{\bm{g}}=h_{\mathcal{S}\bm{g}}$. This is consistent with the transformation laws

$$u_{\mathcal{S}\bm{g}}^{\parallel}=u_{\bm{g}}^{\parallel},\qquad u_{\mathcal{S}\bm{g}}^{\perp}=\det(\mathcal{S})u_{\bm{g}}^{\perp}.$$ (76)

Hence, we find that $\mathcal{F}_{\bm{g}}$ and $\mathcal{G}_{\bm{g}}$ obey the same symmetry constraints as $u_{\bm{g}}^{\parallel}$ and $u_{\bm{g}}^{\perp}$, respectively. For a height profile with $\mathcal{C}_{3v}$ symmetry, we note that $\mathcal{G}_{\bm{g}}$ vanishes for the first-star because it transforms as a pseudoscalar:

$$\mathcal{G}_{\mathcal{M}_{x}\bm{g}}=-\mathcal{G}_{\bm{g}},$$ (77)

which implies that $\mathcal{G}_{\bm{g}}=0$ for $\mathcal{M}_{x}\bm{g}=\bm{g}$. This is the case for any multiple of the first-star reciprocal vector $\bm{g}_{1}$. Hence, together with $\mathcal{C}_{3}$ this implies that for a first-star in-plane displacement field, and a buckling profile with $C_{3v}$ symmetry, there is no coupling to the rotational piece. However, this is only strictly true for the aligned (or anti-aligned) case, as the symmetry is reduced to $C_{3}$ when the twist angle is finite. We assume that these corrections, which are proportional to imaginary part of the second star ($\sin\gamma_{2}$) in lowest order, can be neglected.

C.2 Buckling equation

In order to find the out-of-plane displacement field under periodic boundary conditions for a given in-plane displacement field (specified by $u_{\bm{g}}^{\parallel}$ and $u_{\bm{g}}^{\perp}$), we minimize the elastic energy of a single layer with respect to $h_{-\bm{g}}$. To this end, we first compute

$$\frac{\partial f_{ij,\bm{g}^{\prime}}}{\partial h_{-\bm{g}}} =\left[g_{i}(g_{j}+g_{j}^{\prime})+g_{j}(g_{i}+g_{i}^{\prime})\right]h_{\bm{g}+\bm{g}^{\prime}},$$ (78)

$$\frac{\partial\mathcal{F}_{\bm{g}^{\prime}}}{\partial h_{-\bm{g}}} =\frac{2\left(\bm{g}\times\bm{g}^{\prime}\right)^{2}}{g^{\prime 2}}\,h_{\bm{g}+\bm{g}^{\prime}},$$ (79)

$$\frac{\partial\mathcal{G}_{\bm{g}^{\prime}}}{\partial h_{-\bm{g}}} =\frac{\left(\bm{g}\times\bm{g}^{\prime}\right)\bm{g}^{\prime}\cdot\left(2\bm{g}+\bm{g}^{\prime}\right)}{g^{\prime 2}}\,h_{\bm{g}+\bm{g}^{\prime}}.$$ (80)

Setting the partial derivative of the elastic energy of a single laver with respect to $h_{-\bm{g}}$ equal to zero then gives

$$\kappa g^{4}h_{\bm{g}}+\sum_{\bm{g}^{\prime}} \Bigg\{-2\mu\Bigg[\frac{2\left(\bm{g}\times\bm{g}^{\prime}\right)^{2}}{g^{\prime 2}}\,u_{\bm{g}^{\prime}}^{\parallel}+\frac{\left(\bm{g}\times\bm{g}^{\prime}\right)\bm{g}^{\prime}\cdot\left(2\bm{g}-\bm{g}^{\prime}\right)}{g^{\prime 2}}\,u_{\bm{g}^{\prime}}^{\perp}$$ (81)

$$+g_{x}(g_{x}-g_{x}^{\prime})f_{yy\bm{g}^{\prime}}+g_{y}(g_{y}-g_{y}^{\prime})f_{xx\bm{g}^{\prime}}+\left(g_{x}g_{y}^{\prime}+g_{y}g_{x}^{\prime}-2g_{x}g_{y}\right)f_{xy\bm{g}^{\prime}}\Bigg]$$ (82)

$$+\left(\lambda+2\mu\right)\bm{g}\cdot\left(\bm{g}-\bm{g}^{\prime}\right)\left(u_{\bm{g}^{\prime}}^{\parallel}+\frac{f_{xx\bm{g}^{\prime}}+f_{yy\bm{g}^{\prime}}}{2}\right)\Bigg\}h_{\bm{g}-\bm{g}^{\prime}}=0,$$ (83)

which has to be solved numerically in general.