Ultra-flat twisted superlattices in 2D heterostructures

Moir\'e-superlattices are ubiquitous in 2D heterostructures, strongly influencing their electronic properties. They give rise to new Dirac cones and are also at the origin of the superconductivity observed in magic-angle bilayer graphene. The modulation amplitude (corrugation) is an important yet largely unexplored parameter in defining the properties of 2D superlattices. The generally accepted view is that the corrugation monotonically decreases with increasing twist angle, while its effects on the electronic structure diminish as the layers become progressively decoupled. Here we found by lattice relaxation of around 8000 different Moir\'e-superstructures using high scale Classical Molecular Simulations combined with analytical calculations, that even a small amount of strain can substantially change this picture, giving rise to more complex behavior of superlattice corrugation as a function of twist angle. One of the most surprising findings is the emergence of an ultra-flat phase that can be present for arbitrary small twist angle having a much lower corrugation level than the decoupled phase at large angles. A possible experimental realization of the ultra-flat state is revealed by Scanning Tunneling Microscopy (STM) investigations of the graphene/graphite system.


INTRODUCTION
When 2D materials are layered on top of each other an interference pattern called the Moiré-pattern is formed due to the lattice mismatch and relative rotation between the layers. The Morié-pattern acts as an additional external periodic potential, which perturbs the electronic states of the layers leading to a large variety of physical phenomena such as secondary Dirac cones 1-4 , Hofstadter's butterfly [3][4][5][6] , Brown-Zak oscillations 7 , valley polarized currents 8 , Moiré-excitons [9][10][11][12][13] , and recently superconductive 14 , and Mott-insulating states 15 . The corrugation stemming from the Moiré-pattern is often neglected in theoretical calculations as the periodic potential itself is sufficient to account for the main aspects of these phenomena. However, the way 2D heterostructures relax strain through out-of-plane deformation can highly influence the properties of such systems, enabling us to gain more detailed insights into the origins and rich physics of these phenomena. Experimental evidence for the out-of-plane deformation of a freestanding graphene/h-BN single layer 2D heterostructure has been provided by transmission electron microscopy investigations 16 .
Corrugated Moiré-patterns have nonuniform strains which lead to pseudo-magnetic fields and consequently an energy gap in the graphene in the range of 20-60 meV 17,18 . Heavily corrugated Moiré in gr/Ru(0001) shows localized states close to the Fermilevel 19 . It has been shown that lattice relaxation, (which leads to a corrugated surface) is an important parameter to understand the origin of the gap 20,21 and flat-band formation 22 in gr/h-BN systems. The importance of the relaxation was also pointed out for low angle twisted bilayer graphene (TBLG), where an energy gap of up to 20 meV opens at the superlattice subband edge 18,23 , which is also found in experiments 15,24 and is needed to fully describe the Mott-physics of the system. Small-angle TBLG under interlayer bias also shows topologically protected helical modes in stacking boundaries [25][26][27] . Increasing the corrugation of the stacking boundaries cause quantum-well-like gapped states to abound, which give rise to robust conductance channels along the corrugated boundaries 28 . There is also a quest to synthesize ultraflat graphene which is driven by the fact, that height fluctuations can cause charge inhomogeneities in the graphene which degrade its transport properties 29,30 . On crystalline atomically flat surfaces the main height fluctuation arises from the height modulation of the Moiré-pattern.
Despite the role of the Moiré-corrugation in a wide range of phenomena it is still a largely unexplored field and is lacking a comprehensive theoretical description. The commonly accepted picture is that the corrugation is a monotonically decreasing function of the twist angle. This is generally due to the decreasing Moiré-wavelength with increasing twist angle. Moiré-hills with smaller wavelength would have a higher aspect ratio, and therefore a smaller radius of curvature, which would heavily increase the elastic energy if the corrugation were not reduced to compensate for this effect. Few theoretical calculations have been developed to describe the amplitude modulation [31][32][33][34] , however, the substrate was considered to be rigid and/or perfectly flat in these works. These assumptions hinder some of the complex behavior of the Moiré-corrugation as we will demonstrate in this work.
We used high scale Classical Molecular Mechanics simulations for relaxating approximately 8000 differently strained and twisted Moiré-supercells. Our model systems were the twisted trilayer graphene (T3LG) and twisted five-layer graphene (T5LG), both consisting of Bernal stacked graphene layers with one rotated layer on top. Such few-layer graphene structures turned out to be ideal model systems for studying topographical phases of Moirépatterns. We show that by including strain in the system three typical Moiré-topographies (Moiré-phases) can be present in 2D heterostructures. The existence of various Moiré-phases enriches 1 Institute of Technical Physics and Materials Science (MFA), Centre for Energy Research, Hungarian Academy of Sciences, P.O. Box 49, Budapest 1525, Hungary. ✉ email: szendro. marton@energia.mta.hu the corrugation physics causing that even small strain perturbations lead to complex corrugation behavior.

Moiré-phases and topographical phase transitions
To asses the effect of external strain on the Moiré-corrugation we performed lattice relaxations with differently strained and twisted T3LG and T5LG heterostructures. We used the T3LG structure for the more computation-intensive calculations. Twisted few-layer graphene structures proved to be suitable model systems due to the following reasons: (i) we could use DFT validated force fields for these systems (see Supplementary Fig. 2), (ii) it also captures the asymmetry characteristic of 2D layer/substrate systems as opposed to bilayer graphene, (iii) every Moiré-phase is reachable for reasonable strains in a wide range of twist angles (as we will show later).
We applied homogeneous biaxial strain δ independently to the overlayer (δ o ) i.e., the top layer and to the substrate (δ s ) i.e., the remaining layers. Firstly, we calculated the Moiré-corrugation of relaxed T5LG superstructures as a function of twist angle in the range of α = 2 ∘ -16 ∘ for three distinct strain configurations (see Fig. 1). The most striking finding is that the strain crucially influences the behavior of the corrugation curves. In the first case ( Fig. 1 a) the overlayer was compressed by δ o = −0.65% but the substrate was left strain-free δ s = 0%. The corrugation curve of the overlayer shows the expected monotonically decreasing behavior. In the second case (Fig. 1b) the overlayer was stretched by δ o = +0.45%, still leaving the substrate with no external strain. The overlayer corrugation is no longer monotonically decreasing in the whole α range, rather it tends to vanish at around 4 ∘ , then increasing between 4 ∘ -7 ∘ before reaching again the decreasing tendency characteristic to large twist angles. As for the third case ( Fig. 1c) both the substrate and the overlayer were stretched by a small amount of δ o = +0.05%, δ s = +0.05%, which leads to a curve having a plateau in the interval of 0 ∘ -4 ∘ and then monotonically decreasing in the rest. The significant differences of the corrugation curves triggered by very mild strain levels hint an intriguing underlying mechanism to be revealed. Also for a particular α = 4 ∘ , δ o = +0.45%, δ s = 0% parameter set an ultra-flat state surprisingly appears in the overlayer. In order to get a better insight, we investigated in details the resulting topographies. We found that during rotation, not only the periodicity and corrugation but also the morphology of the Moiré-superlattice can change. Such Moiré-phase transitions are marked by vertical dashed lines in Fig. 1. Sometimes these phase transitions do not manifest themselves in corrugation change ( Fig. 1a), while in other cases the corrugation can drastically alter around them ( Fig. 1 b). The corresponding topographies of the top three layers of the phases are depicted in bubbles in Fig. 1.
In 2D heterostructures, corrugation can develop in both upper and lower layers. The Moiré-amplitude in the substrate can be in phase with the overlayer amplitude (bending modes), or in antiphase (breathing mode). There are two bending modes (convex, nanomesh) and one breathing (chain), which we call together Moiré-phases. In the convex-phase the AA regions of the Moirépattern are forming protrusions, while the remaining two AB/BA stackings are depressions (Fig. 1a). The nanomesh-phase is the inverse of this: the AA regions are depressions and AB/BA stackings are the protrusions (Fig. 1b). In the convex and nanomesh phases, both the overlayer and the substrate has the same type of topography (bending-modes). However, the third phase, the chain-phase ( Fig. 1a, b, c) is qualitatively different as the overlayer and the substrate are closely mirror images of each other having a convex topography in the overlayer and a concave in the substrate (breathing mode).
Moiré-phase maps, criticality, ultra-flat states The results above clearly reveal that for different α, δ o , δ s parameters different Moiré-phase realizations occur, and even ultra-flat states can emerge near some of the phase transitions. To capture all the Moiré-phases emerging in our system, as well as the conditions necessary for the realization of ultra-flat states, the whole (α, δ o , δ s ) phase space has to be studied by examining the resulting topographies and corrugations for each point in this space. Therefore we created phase-maps of the Moiré-phases in the space of (δ o , δ s ) for four different twist angles α = 1.1 ∘ , 2.5 ∘ , 4.0 ∘ , 8.2 ∘ . The sampling of the whole phase-space needs a large number of superstructure relaxations, therefore we used the fewer layer T3LG structure for this task instead of the T5LG. For each phase-map with a definite α, we searched for about 60 Moiré T3LG supercells having different δ o , δ s combinations. Then we stretched each supercell structure in steps of Δδ ≈ 0.085%. This way we ended up with approximately 1800-2000 Moiré-supercells for each phase-map. Every Moiré T3LG superstructure was relaxed with a fixed simulation cell using in-plane periodic boundary conditions and the same setup as in the case of Fig. 1. The type of morphology and corrugation was analyzed after the relaxations. Moiré-superlattice corrugations of the top rotated graphene layer for three differently strained T5LG (one twisted graphene layer above four Bernal stacked graphene layers) heterostructures (red dots) from molecular mechanics simulations. Each point represents a commensurate Moiré-superlattice. a The overlayer is externally compressed with~0.65% during minimization which leads to a convex phase for angles below 8.2 ∘ and a chain phase for angles above. The phase transition, indicated by dashed lines, is smooth, and the whole curve is monotonically decreasing. b The overlayer is stretched with~0.45% strain. This cause a nanomesh phase to emerge for angles below 4 ∘ , and a chain phase for angles above. Near the phase transition (dashed line), the corrugation vanishes and an ultra-flat phase appears. c Both the overlayer and the substrate are stretched with a small 0.05% strain. No phase-transition occurs, each point corresponds to a chain structure. In the range of 0-4 ∘ the corrugation is closely constant (plateau-effect). The topographies of the top three layers of the corresponding Moiré-phases are displayed in bubbles.
The results of the~8000 Moiré-superstructures can be seen in Fig. 2.
Three Moiré-phases are present in the phase-maps for the examined angles (Fig. 2 top row). Also for compressive strains higher than a threshold value, the Moiré-pattern is absent, and the whole structure is buckled due to the buckling instability of graphene while having a much larger corrugation. We denoted the buckled states with black. Incidentally, the buckled states are twist-angle dependent: increased buckling threshold can be seen in the convex phase for higher twist angles. There are three phaseboundaries between the Moiré-phases. Two out of the three are more extended, namely the nanomesh-chain boundary and the convex-chain boundary. The convex-nanomesh boundary is small and is vanishing for larger twist angles as this boundary is completely merging with the buckled states. Furthermore, there is a triple point where the three phases coexist. The fascinating thing is that the position of the strain-free δ o = 0, δ s = 0 state is very close to the triple point for angles between 0 ∘ -4 ∘ . What this means is that very delicate strain fluctuations can lead to very different phase realizations, therefore the δ o = 0, δ s = 0 state is very sensitive from this point of view. Not only it is sensitive for the realization of various Moiré-phases but from a corrugation perspective. In the bottom row of Fig. 2, we show the corrugation value of each relaxed Moiré-supercell. When stretching the whole system gradually the corrugation is decreasing as it is expected; however, for a certain set of finite δ o , δ s , the corrugation drops rapidly close to zero (dark regions in the bottom row of Fig. 2). These are the already discussed ultra-flat states. They run along the convex-chain and the nanomesh-chain phase boundaries. The main difference in the two boundaries is that in the former case the substrate corrugation is zero, while in the latter the overlayer. This is the reason why we could not see ultra-flat states in Fig. 1 a) as these states occur in the substrate. The emergence of ultra-flat states can be linked to the convexity change of the topographies.
The convexity of the overlayer is positive both in the convex and in the chain-phase, therefore the transition between these two phases can be continuous in the corrugation. However, when moving from nanomesh to chain the convexity changes sign progressively in the overlayer implying a particular twist angle in between, where the corrugation has to be zero. Nonetheless, this effect takes place in the substrate inversely: the change of the convexity occurs during the convex-chain transition with consequent ultra-flat states appearing in the substrate, but not in the overlayer (see also Supplementary Fig. 7).
The ultra-flat states go along the boundaries, however, they do not meet in the middle at δ o = 0, δ s = 0. Now it can be clearly seen, that around the δ o = 0, δ s = 0 point, the corrugation strongly fluctuates. It means that in realistic samples even small strain fluctuations can be crucial in determining the topography of the system, which can take shape from being ultra-flat to highly corrugated, or with convexity from positive to negative. This is a very important finding for the practical design of 2D heterostructures, pointing out the crucial role of even mild strains that have to be controlled very accurately for reliably engineering the morphology of Moiré-superlattices in 2D heterostructures.
The phase boundaries (convex-chain, nanomesh-chain) are also moving as the twist angle changes. With increasing twist angle we see the chain phase being more dominant on the phase map. The movement of the boundaries takes place in two distinct stages. For smaller twist angles (0 ∘ -4 ∘ ) only the angle between the phase boundaries are changing and the boundaries are rotating. For larger angles (>4 ∘ ) the phase boundaries are perpendicular to each other, and their movement consists of pure translation. The movement of the phase boundaries makes it possible for a system to go under a phase transition while altering the twist angle as it was shown in Fig. 1.

Energetic interpretation
The interpretations of the results above can be attained through energetic considerations. The Moiré-morphology is defined by the balance between the elastic energy and the vdW adhesion. The elastic energy stemming from the corrugation ξ with a Moiréwavelength λ is a power series of ξ/λ 33 . The vdW contribution is determined by the interlayer separation of the AA stacking d 1 and the AB/BA stackings d 2 . The corrugation in the overlayer ξ o and substrate ξ s are coupled through vdW interaction. In particular, it is easy to show, that ξ o = ξ s + Δd (see Supplementary Fig. 3), where Δd = d 1 − d 2 , which is a strict constraint on the geometry of the phases. It can be seen (see Section A of Supplementary Discussion), that in the convex phase ξ o ≥ Δd and ξ o > ξ s , while ξ s ≥ Δd and ξ o < ξ s in the nanomesh phase. In the chain phase 0 ≤ ξ o ≤ Δd, 0 ≤ ξ s ≤ Δd. These geometric properties have a strong impact on the realization of the phases. For example, a more corrugated overlayer than the substrate can be realized only in the convex or the chain phase. In the case when the overlayer does not want to corrugate as much as the substrate then the nanomesh or the chain phase is favorable. When the corrugation is equally unfavorable for the overlayer and for the substrate, we expect the chain phase to emerge. Also as the elastic energy goes with the powers of ξ/λ, the energy increment with a unit dξ will be much bigger for larger twist angles. Therefore, there is a tendency to divide the total corrugation evenly between the substrate and the overlayer as the twist angle increases. This geometry, however, can only be realized in the chain phase. Consequently, with increasing twist angle we expect to see the chain phase being more dominant on the phase maps. Moreover, for small twist angles, ξ/λ tends to zero and the vdW adhesion becomes prominent, while for large twist angles the elastic energy is the major energy term (see Supplementary Fig. 9). The plateau on Fig.  1c) is the consequence of the special geometry of the chain-phase combined with these two energetic regimes (Supplementary Fig.  9). The boundary movements (rotation, translation) also coincides with these two regimes. The detailed analysis of the results on this energetic basis can be found in the Supplementary  Scanning tunneling microscopy The experimental validation of our theoretical predictions is rather difficult as we cannot fully control the amount of heterostrain emerging in the system. Nevertheless, we have prepared and analyzed graphene layers deposited at various rotation angles on graphite substrates. We were able to measure the rotation angle, Moiré-periodicity, as well as the Moiré-amplitude for various rotation angles using room temperature scanning tunneling microscopy. Although, at finite temperatures, thermal fluctuations are present in two dimensional heterostructures, distorting the Moiré superlattice, such fluctuations are rapidly averaged out on the time scale of scanning tunneling microscopy imaging, therefore, an undistorted Moiré can be observed. On shorter time-scales, the finite temperature time-averaged corrugations of the Moiré superlattices are scattered around the corrugation of the 0K simulation 35 . Albeit, due to the lack of control over the heterostrain, we cannot reliably plot curves similar to the theoretical plots in Fig. 1, the general tendency observed was that the corrugation is high at small rotation angles and low at high angles, which is in accordance with the previous expectations. However, we were also able to find a graphene flake where the corrugation was very low almost undetectable at an intermediate (9 ∘ ) rotation angle. We made sure that this is not an artifact, as the atomic corrugation was as expected (no resolution issue). Furthermore, a very faint superlattice signature could still be observed in larger area images, ensuring that the graphene overlayer is not decoupled by contamination from the graphite substrate. The experimentally determined average Moiréamplitudes for the three rotation angles shown in Fig. 3 are as follows: 0.32 Å ± 0.02 Å(4.2 ∘ ± 0.1 ∘ ); 0.05 Å ± 0.02 Å(9.7 ∘ ± 0.1 ∘ ); 0.12 Å ± 0.02 Å(16 ∘ ± 0.5 ∘ ). Consequently, the graphene flake displayed in Fig. 3b) can be regarded as an experimental realization of an ultra-flat state predicted theoretically, although the exact parameters (heterostrain values) could not be directly inferred. Therefore the direct quantitative comparison of the experimentally found ultra-flat state with the simulations is not straightforward as the unknown strain values essentially influence the twist angle at which the ultra-flat states emerge as it was shown in Fig. 2. The exact shape of the vdW interactions also affects the twist angle value of the ultra-flat state (see Supplementary Fig.  8d). These effects can be at the origin of the quantitative discrepancy between the theoretically predicted and experimentally observed twist angles for the ultra-flat states. Changes in the Moiré-morphologies (i.e., from convex to concave) have also been observed experimentally ( Supplementary Fig. 1) on the same graphene flake (for the same rotation angle), which can be induced by the locally varying strain, in agreement with the predicted high strain-sensitivity of the Moiré-phase realizations. The high strain-sensitivity of Moiré-superlattices is also supported by the experimental observation of the coexistence of the convex Fig. 3 Scanning tunneling microscopy observation of an ultra-flat state. Topographic scanning tunneling microscope images of graphene layers deposited on top of a graphite substrate for various relative rotation angles (scale bars: 2 nm). Although the apparent corrugation has a decreasing tendency from a-c, the STM image in panel b reveals an ultra-flat state, much smoother than observed even for high rotation angles (c). The experimental conditions for image acquisition (I tunnel = 1 nA, U bias = 200 mV), data processing, and graphic display parameters are the same for all panels. The relative rotational angles are indicated under each panel. and nanomesh morphologies of graphene superlattices on Au (111) 36 .

Conclusions
In conclusion, we showed using Molecular Mechanics simulations that strain and twist angle alteration can induce phase transitions in twisted 2D superlattices between three typical topographical Moiré-phases. On phase boundaries, ultra-flat states were identified where the corrugation is almost zero either in the overlayer or in the substrate, but not both. These ultra-flat states were shown to be also present for small twist angles, where conventionally a large corrugation is expected. We also showed that for smaller twist angles (0 ∘ -4 ∘ ) the strain-free equilibrium point is critical and even very small strains can induce significant changes in the corrugation and in the morphology of the superlattices. The results and concepts developed here are expected to be essential for the rational design of twisted 2D superlattices and highlight the critical role played by externals strains in defining the nanoscale morphology of such 2D heterostructures.

Classical molecular geometry optimalization
The classical molecular energy minimizations were carried out using the LAMMPS code 37 . In-plane periodic boundary conditions were used. During the relaxation the size of the simulation cell was fixed, which ensured through the periodic boundary conditions that the pre-applied strain can not be fully eliminated, and that δ o , δ s is an additional constraint on the energy minimum. The long range bond order potential was utilized 38 for carbon-carbon interactions within a single carbon layer. The weak Van der Walls forces between the carbon layers were modeled using the Kolmogorov-Crespi potential 39 . We performed geometry relaxations on commensurate Moiré-supercells using Hessian-free Truncated Newton algorithm implemented in LAMMPS as the hftn algorithm. Each supercell used in this paper was found by our script iterating through highly ordered commensurate cells as described in refs. 40,41 to match the criterion of the cell having a definite α, δ o , δ s property. Due to the fact that for a given arbitrary α, δ o , δ s there not necessarily exists a corresponding supercell, the α, δ o , δ s values may display minor fluctuations that limit the angle and strain resolution of our data; nevertheless do not affect their interpretation.

DATA AVAILABILITY
The data supporting this work are available from the corresponding author upon reasonable request.