Formal justification of a continuum relaxation model for one-dimensional moiré materials

This work was carried out during JZ’s REU project at University of Minnesota during Summer 2024 supervised by ABW. JZ and ABW would like to thank Michael Hott and Mitchell Luskin for stimulating discussions.

We start by defining a natural energy functional for the bilayer system defined through interatomic potentials. For simplicity, we restrict attention to pair potentials. The approach here is influenced by [1].

Just as before, we assume the energy has the form

$$E_{\text{total}}(u_{1},u_{2}):=E_{\text{intra}}(u_{1},u_{2})+E_{\text{inter}}(u_{1},u_{2}).$$

(20)

For $E_{\text{intra}}$, we take

$$\begin{split}&E_{\text{intra}}(u_{1},u_{2}):=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}W_{\text{intra}}\left(ai+u_{1}(ai)-aj-u_{1}(aj)\right)\\ &+\frac{1}{M+1}\sum_{i=1}^{M+1}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}W_{\text{intra}}\left(\frac{a(1-\theta)i+u_{2}(a(1-\theta)i)-a(1-\theta)j-u_{2}(a(1-\theta)j)}{1-\theta}\right),\end{split}$$

(21)

where $W_{\text{intra}}$ is an interatomic pair potential. Note that we assume periodic boundary conditions, so that we must allow for interactions between atoms within the moiré supercell and with their copies in all other cells. The sum over $j$ converges because of decay of $W_{\text{intra}}$.

For $E_{\text{inter}}$, we take

$$\begin{split}&E_{\text{inter}}(u_{1},u_{2}):=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{j=-\infty}^{\infty}W_{\text{inter}}\left(ai+u_{1}(ai)-a(1-\theta)j-u_{2}(a(1-\theta)j)\right),\end{split}$$

(22)

where $W_{\text{inter}}$ is another interatomic pair potential. Again, since we assume that $W_{\text{inter}}$ decays, the sum over $j$ converges. We assume at this point that $W_{\text{intra}}$ and $W_{\text{inter}}$ are both even so that

$$W_{\text{intra}}(-z)=W_{\text{intra}}(z),\quad W_{\text{inter}}(-z)=W_{\text{inter}}(z).$$

(23)

We now nondimensionalize the model with respect to length. We introduce nondimensionalized displacement functions as

$$U_{i}(X):=a^{-1}u_{i}(a_{\mathcal{M}}X),\quad i=1,2,$$

(24)

and re-scale the interatomic potentials $W_{\text{intra}}$ and $W_{\text{inter}}$ as

$$\tilde{W}_{\text{intra}}(X):=W_{\text{intra}}(aX),\quad\tilde{W}_{\text{inter}}(X):=W_{\text{inter}}(aX).$$

(25)

With these scalings, the displacement functions $U_{i}$ are now $1$-periodic, and the intralayer energy becomes

$$\begin{split}&E_{\text{intra}}(U_{1},U_{2})=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(i-j+U_{1}(\epsilon i)-U_{1}(\epsilon j)\right)\\ &+\frac{1}{M+1}\sum_{i=1}^{M+1}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(i-j+\frac{U_{2}(\epsilon(1-\theta)i)-U_{2}(\epsilon(1-\theta)j)}{1-\theta}\right),\end{split}$$

(26)

where we introduce notation for the dimensionless length ratio

$$\epsilon:=\frac{a}{a_{\mathcal{M}}}.$$

(27)

Note that (3) implies that

$$\epsilon=\frac{\theta}{1-\theta}$$

(28)

so that when we take $\epsilon\downarrow 0$ below we must take $\theta\downarrow 0$ with $\epsilon\sim\theta$ as $\epsilon\downarrow 0$ as well. The interlayer energy becomes

$$\begin{split}&E_{\text{inter}}(U_{1},U_{2})=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{j=-\infty}^{\infty}\tilde{W}_{\text{inter}}\left(i-(1-\theta)j+U_{1}(\epsilon i)-U_{2}(\epsilon(1-\theta)j)\right).\end{split}$$

(29)

We have not non-dimensionalized with respect to energy yet: both $\tilde{W}_{\text{intra}}$ and $\tilde{W}_{\text{inter}}$ still carry units of energy. We postpone this until we have simplified the model further.

We now simplify the energy functional and take its continuum limit, an approximation which should be valid in the limit $\epsilon\downarrow 0$. Recalling (4), this is equivalent to taking the limit $a\downarrow 0$ and $M\rightarrow\infty$ at fixed $a_{\mathcal{M}}$.

Before taking the limit, note that by re-defining $j$ and using (23) the intralayer energy can be re-written as

$$\begin{split}&E_{\text{intra}}(U_{1},U_{2})=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq 0\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(j+U_{1}(\epsilon i+\epsilon j)-U_{1}(\epsilon i)\right)\\ &+\frac{1}{M+1}\sum_{i=1}^{M+1}\sum_{\begin{subarray}{c}j=-\infty\\ j\neq 0\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(j+\frac{U_{2}(\epsilon(1-\theta)i+\epsilon(1-\theta)j)-U_{2}(\epsilon(1-\theta)i)}{1-\theta}\right).\end{split}$$

(30)

Taking the limit, assuming smoothness of the integrands, both sums over $i$ converge to integrals from $0$ to $1$ as

$$\begin{split}&\int_{0}^{1}\!\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(j+U_{1}(X+\epsilon j)-U_{1}(X)\right)\,\textrm{d}X\\ &+\int_{0}^{1}\!\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(j+\frac{U_{2}(X+\epsilon(1-\theta)j)-U_{2}(X)}{1-\theta}\right)\,\textrm{d}X.\end{split}$$

(31)

We can further approximate the energy functional by Taylor-expanding the $U_{i}$s in powers of $\epsilon$ and keeping only the leading terms as

$$\int_{0}^{1}\!\tilde{W}_{\text{CB}}\left(\epsilon\partial_{X}U_{1}(X)\right)+\tilde{W}_{\text{CB}}\left(\epsilon\partial_{X}U_{2}(X)\right)\,\textrm{d}X,$$

(32)

where

$$\tilde{W}_{\text{CB}}(z):=\sum_{\begin{subarray}{c}j=-\infty\\ j\neq i\end{subarray}}^{\infty}\tilde{W}_{\text{intra}}\left(j(1+z)\right)$$

(33)

is the Cauchy-Born energy density. We can then further approximate by Taylor-expanding $\tilde{W}_{\text{CB}}$ in powers of $\epsilon\partial_{X}U_{i}$. Assuming that constant displacement is a non-degenerate minimizer for $\tilde{W}_{\text{CB}}$, we have

$$\tilde{W}_{\text{CB}}(z)=\tilde{W}_{\text{CB}}(0)+\frac{\tilde{\kappa}z^{2}}{2}+o(z^{2})\quad z\downarrow 0,$$

(34)

where $\tilde{\kappa}:=\partial_{z}^{2}\tilde{W}_{\text{CB}}(0)$ again characterizes the stiffness of the lattices. Since its value does not affect the form of minimizers, we at this point set $\tilde{W}_{\text{CB}}(0)=0$ without loss of generality. We arrive at the approximation

$$E_{\text{intra}}(U_{1},U_{2})\approx\frac{\tilde{\kappa}\epsilon^{2}}{2}\int_{0}^{1}\!\left[\left(\frac{\partial U_{1}}{\partial X}\right)^{2}+\left(\frac{\partial U_{2}}{\partial X}\right)^{2}\right]\,\textrm{d}X.$$

(35)

We claim that the derivation leading to (35) provides a microscopic justification of the linear elasticity energy (10) introduced in the GSFE model. To see this, note that upon inverting the transformation (24) and setting¹¹1To see why this makes sense, note that $\kappa$ is defined as an energy per unit length, while $\tilde{\kappa}$ is defined through the derivative of an energy with respect to a nondimensionalized length and hence has units only of energy. $\tilde{\kappa}=a_{\mathcal{M}}\kappa$, we obtain exactly (10).

We now consider the interlayer energy. Again, by re-defining $j$, this energy can be re-written as

$$\begin{split}&E_{\text{inter}}(U_{1},U_{2})=\\ &\frac{1}{M}\sum_{i=1}^{M}\sum_{j=-\infty}^{\infty}\tilde{W}_{\text{inter}}\left(\theta i-(1-\theta)j+U_{1}(\epsilon i)-U_{2}(\epsilon(1-\theta)(i+j))\right).\end{split}$$

(36)

In order to take the continuum limit we write everything in terms of the variable $\epsilon i$, so that

$$\theta i=\frac{\theta}{\epsilon}\epsilon i.$$

(37)

Recalling (3) and (27), we have that

$$\frac{\theta}{\epsilon}\epsilon i=(1-\theta)\epsilon i.$$

(38)

Substituting this and taking the continuum limit we obtain

$$\int_{0}^{1}\!\sum_{j=-\infty}^{\infty}\tilde{W}_{\text{inter}}\left((1-\theta)X-(1-\theta)j+U_{1}(X)-U_{2}((1-\theta)X+\epsilon(1-\theta)j))\right)\,\textrm{d}X.$$

(39)

Just as with the intralayer energy, we expand $U_{2}$ in powers of $\epsilon$ and keep only the leading term. This yields

$$\int_{0}^{1}\!\tilde{V}\left((1-\theta)X+U_{1}(X)-U_{2}((1-\theta)X))\right)\,\textrm{d}X,$$

(40)

where we have introduced the $(1-\theta)$-periodic stacking energy

$$\tilde{V}(z):=\sum_{j=-\infty}^{\infty}\tilde{W}_{\text{inter}}\left(z-(1-\theta)j\right).$$

(41)

Finally, we use smallness of $\theta$ to Taylor-expand $U_{2}$, to arrive at

$$\int_{0}^{1}\!\tilde{V}\left((1-\theta)X+U_{1}(X)-U_{2}(X))\right)\,\textrm{d}X.$$

(42)

We now claim that the derivation of (42) above provides a microscopic justification of (11). First, note that inverting the transformation (24) and setting²²2Here we multiply by $a_{\mathcal{M}}$ for the same reason as when relating $\kappa$ and $\tilde{\kappa}$.

$$\tilde{V}(X)=a_{\mathcal{M}}V(aX),$$

(43)

the energy (42) becomes

$$a_{\mathcal{M}}\int_{0}^{1}\!V\left((1-\theta)aX+u_{1}(X)-u_{2}(X))\right)\,\textrm{d}X.$$

(44)

Changing variables in the integral as $x=a_{\mathcal{M}}X$, after recognizing (3) that $\theta=\frac{a(1-\theta)}{a_{\mathcal{M}}}$, yields

$$\int_{0}^{a_{\mathcal{M}}}\!V\left(\theta x+u_{1}(x)-u_{2}(x))\right)\,\textrm{d}x,$$

(45)

which is exactly (11). Note that $(1-\theta)$-periodicity of $\tilde{V}$ is exactly consistent with the $(1-\theta)a$-periodicity of $V$.

We note further that, at the cost of additional error proportional to $\theta$, we can simplify (42) to

$$\int_{0}^{1}\!\tilde{V}\left(X+U_{1}(X)-U_{2}(X))\right)\,\textrm{d}X,$$

(46)

where $\tilde{V}$ is $1$-periodic. If we also simplify $\tilde{V}$ to the form of a single sinusoid as we did for the GSFE model we obtain the approximation

$$E_{\text{inter}}(U_{1},U_{2})\approx-2\tilde{V}_{0}\int_{0}^{1}\!\cos\left(2\pi\left(X+U_{1}(X)-U_{2}(X))\right)\right)\,\textrm{d}X,$$

(47)

where $\tilde{V}_{0}:=a_{\mathcal{M}}V_{0}$.

To summarize the results of this section, we have obtained an approximation of the atomistic energy functional as the sum of (35) and (47). Just as we did with the GSFE model, it is natural to write the functional as a function of new unknowns

$$U_{+}(x):=\frac{U_{1}(x)+U_{2}(x)}{\sqrt{2}},\quad U_{-}(x):=\frac{U_{1}(x)-U_{2}(x)}{\sqrt{2}},$$

(48)

and set (WLOG) $U_{+}=0$. We thus derive an approximation of the atomistic energy functional depending only on $U_{-}$ as

$$E_{\text{total}}(U_{-})\approx\int_{0}^{1}\!\left[\frac{\epsilon^{2}\tilde{\kappa}}{2}\left(\frac{\partial U_{-}}{\partial X}\right)^{2}-2\tilde{V}_{0}\cos\left(2\pi\left(X+\sqrt{2}U_{-}(X)\right)\right)\right]\,\textrm{d}X,$$

(49)

which is nothing but a partially nondimensionalized version of (18).

After the simplifications made in the previous section, we have that the intralayer part of the functional (49) depends only on the energy scale $\tilde{\kappa}=\kappa a_{\mathcal{M}}$, while the interlayer part depends only on the energy scale $\tilde{V}_{0}=V_{0}a_{\mathcal{M}}$. It is natural then to non-dimensionalize with respect to energy by introducing the dimensionless energy

$$\tilde{E}_{\text{total}}:=\frac{E_{\text{total}}}{\tilde{V}_{0}}.$$

(50)

We thus arrive at the fully nondimensionalized model

$$\tilde{E}_{\text{total}}(U_{-})\approx\int_{0}^{1}\!\left[\frac{\epsilon^{2}}{2\delta}\left(\frac{\partial U_{-}}{\partial X}\right)^{2}-2\cos\left(2\pi\left(X+\sqrt{2}U_{-}(X)\right)\right)\right]\,\textrm{d}X,$$

(51)

where we have introduced notation for the dimensionless energy ratio

$$\delta:=\frac{\tilde{V}_{0}}{\tilde{\kappa}}=\frac{V_{0}}{\kappa}.$$

(52)

It is clear that to have a model where the intralayer and interlayer terms balance in the limit $\epsilon\downarrow 0$, i.e., where neither term is asymptotically larger in this limit, we require that $\delta\downarrow 0$ in such a way that the dimensionless ratio

$$\eta:=\frac{\sqrt{\delta}}{\epsilon}=\sqrt{\frac{V_{0}}{\kappa}}\frac{a_{\mathcal{M}}}{a},$$

(53)

remains fixed. The form of minimizers in the limit, moreover, depends only on the value of $\eta$ as the limit is taken. This justifies Nam and Koshino’s observation [8] that minimizers of the GSFE model depend only on the value of this parameter; recall (19).

It is interesting to again consider physical values of $\epsilon$ and $\delta$ and check whether, for example, twisted bilayer graphene at the magic angle is in this regime. Recall that there we have $a\approx.25$ nm, and $\theta\approx\frac{1}{50}\approx 1.1^{\circ}$. We can then estimate

$$\epsilon=\frac{a}{a_{\mathcal{M}}}\approx\frac{.25}{50(.25)}=\frac{1}{50}=0.02.$$

(54)

As for the energy scales, we have $\kappa\approx 50,000$ meV per unit area, and $V_{0}\approx 20$ meV per unit area [4]. We thus have

$$\delta=\frac{V_{0}}{\kappa}\approx\frac{20}{50,000}=\frac{1}{2500}=0.0004,$$

(55)

so that, remarkably,

$$\eta\approx 1$$

(56)

in twisted bilayer graphene. We conclude that, for twisted bilayer graphene near the magic angle, (18) can indeed be expected to be a good model.