Analytical solution for the relaxed atomic configuration of twisted bilayer graphene including heterostrain

TL;DR

The paper analyzes lattice relaxation in twisted bilayer graphene by comparing two continuum models against Bragg interferometry data, identifying the KaxirasPRB18 model as more accurate. It introduces a dimensionless parameter $\lambda = c_1/(\mathcal{G}\theta^2)$ to quantify adhesion–elastic energy balance and derives a closed-form analytical solution for the in-plane relaxation that remains valid up to the first magic angle $\theta \approx 1.05^{\circ}$, while also incorporating external heterostrain through a strain matrix $S^{\epsilon}$. A Taylor-series expansion in $\lambda$ is shown to have a finite radius of convergence, explaining why perturbative approaches fail near the magic angle; the authors then develop a nonperturbative closed-form for the relaxation, with higher-shell corrections $\zeta_2$ and $\zeta_3$ determined from energy minimization. The methodology is extended to heterostrain, yielding a linearized, tractable expression for the strain-induced corrections $\delta \boldsymbol{U}$ and corresponding Fourier amplitudes, accurate to about 1% for strains up to 1%, and applicable to uniaxial/biaxial cases. These results provide a practical input for electronic continuum models and enable strain engineering analyses in moiré materials.

Abstract

Continuum atomic relaxation models for twisted bilayer graphene involve minimization of the sum of intralayer elastic energy and interlayer adhesion energy. The elastic energy favors a rigid twist i.e. no distortion in the twisted honeycomb lattices, while the adhesion energy favors Bernal stacking and breaking the relaxation into triangular AB and BA stacked domains. We compare the results of two relaxation models with the published Bragg interferometry data, finding good agreement with one of the models. We then provide a method for finding a highly accurate approximation to the solution of this model which holds above the twist angle of $\approx0.7^\circ$ and thus covers the first magic angle. We find closed form expressions in the absence, as well as in the presence, of external heterostrain. These expressions are not written as a Taylor series in the ratio of adhesion and elastic energy, because, as we show, the radius of convergence of such a series is too small to access the first magic angle.

Figures (11)

• Figure 1: Lattice relaxation for the twist angle $\theta = 1.05^{\circ}$ with (a) no external heterostrain and (b) $\epsilon_1 = 0.5\times 10^{-2}$, $\epsilon_2 = -0.8 \epsilon_1$, and $\phi = \pi/12$. The moire lattice vectors as well as various stacking points are also marked.
• Figure 2: $\theta_R = \frac{1}{2} \left( \partial_x \delta U_y - \partial_y \delta U_x \right)$ at two positions: SP stacking ("+") at $\boldsymbol{x} = \frac{1}{2} \boldsymbol{L}_1$ and AB stacking (empty square) at $\boldsymbol{x} = \frac{1}{3}(\boldsymbol{L}_1 + \boldsymbol{L}_2)$, when the heterostrain is absent. The three pairs of plots show the experimental measurements (red), the results of the model in Ref. KoshinoPRB17*KoshinoPRB17Erratum (green) KangLatticeRelax2025, and the results of the model in Ref. KaxirasPRB18 (blue) KangLatticeRelax2025. The experimental data is more consistent with the model in Ref. KaxirasPRB18. More distortion parameters are plotted and compared in Fig. \ref{['FigS:RelaxComp']}.
• Figure 3: (a,b) $\boldsymbol{G}_i$ are the reciprocal lattice vectors of the undistorted monolayer graphene, with $\boldsymbol{G}_1$, $\boldsymbol{G}_2$ defined by the Eq. \ref{['Eqn:GVectors']} and $\boldsymbol{G}_3 = -(\boldsymbol{G}_1 + \boldsymbol{G}_2)$. (c) $\boldsymbol{g}_i$ ($i =1$, $2$, $3$) are the reciprocal vectors of the moire lattice given by Eq. \ref{['Eqn:gdef']}; in the (a) and (c) for we did not include heterostrain, only the twist. (b) the first three $\boldsymbol{G}$-shells marked by red (the first shell), black (the second shell), and blue (the third shell) colored arrows. Note that in our definition the $\boldsymbol{G}$ vectors are unaffected by strain, unlike the $\boldsymbol{g}$ vectors. (c) the first three $\boldsymbol{g}$-shells marked by the same colored arrows as for $\boldsymbol{G}$-shells. Unlike in (a), in the schematics (b) and (c) we do not distinguish between different sizes of $\boldsymbol{G}$ and $\boldsymbol{g}$ vectors.
• Figure 4: Numerical solution KangLatticeRelax2025 to the Eq. \ref{['Eqn:LatticeUqRelax']}, minimizing the sum of the elastic energy $U_E$ and the interlayer adhesion energy $U_B$, is obtained first; it includes many $\boldsymbol{g}$-shells and achieves convergence. (a) Red dots show the amplitude of this solution at the first shell $\pm \boldsymbol{g}_{1,2,3}$ for the parameters of Ref. KaxirasPRB18 in Table \ref{['Tab:RelaxPara']}, expressed in terms of the dimensionless variable $\zeta_1$ (see Eq.\ref{['Eqn:LatRelaxFormula']}). The blue curve is the approximate formula given in the Eq.\ref{['Eqn:Zeta1Appr']} using our approach, upon which we can systematically improve as explained in the main text. The solid black curve and the dashed green curve are the Taylor expansions of the solution to the first and the second order in the dimensionless parameter $\lambda=c_1/(\mathcal{G}\theta^2)$ that quantifies the relative magnitudes of $U_B$ and $U_E$ at a given twist angle $\theta$. The Taylor expansions are seen to be accurate above the twist angle $\sim 1.6^{\circ}$, but to increasingly deviate below. In Sec.\ref{['Sec:radius of convergence']} we explain the poor performance of the Taylor series expansion by finding its radius of convergence which translates to the twist angle $1.14^{\circ}$ and placing the magic angle beyond the reach of the Taylor expansion approach. (b) Dots show the analogous numerical results at the second and third shells, demonstrating that in this range of twist angles, the dimensionless amplitudes $\zeta_{2,3}$ are at least an order of magnitude smaller than $\zeta_1$. The solid lines are our approximate formulas; next order improvements can be found in the Fig. \ref{['FigS:LatticeRelax']}.
• Figure 5: Schematic plot for the branch points of the function $\zeta_1(\lambda)$ in the complex $\lambda$-plane, whose magnitude $|\lambda_{BP}|$ dictates the radius of convergence, $\lambda_c$, of the power series expansion of $\zeta_1(\lambda)$ about $\lambda=0$. Note that the two branch points, marked by the blue crosses, are complex conjugates of each other.