Structural relaxation and low energy properties of Twisted Bilayer Graphene

The structural and electronic properties of twisted bilayer graphene are investigated from first principles and tight binding approach as a function of the twist angle (ranging from the first"magic"angle $\theta=1.08^\circ$ to $\theta=3.89^\circ$, with the former corresponding to the largest unit cell, comprising 11164 carbon atoms). By properly taking into account the long-range van der Waals interaction, we provide the patterns for the atomic displacements (with respect to the ideal twisted bilayer). The out-of-plane relaxation shows an oscillating ("buckling") behavior, very evident for the smallest angles, with the atoms around the AA stacking regions interested by the largest displacements. The out-of-plane displacements are accompanied by a significant in-plane relaxation, showing a vortex-like pattern, where the vorticity (intended as curl of the displacement field) is reverted when moving from the top to the bottom plane and viceversa. Overall, the atomic relaxation results in the shrinking of the AA stacking regions in favor of the more energetically favorable AB/BA stacking domains. The measured flat bands emerging at the first magic angle can be accurately described only if the atomic relaxations are taken into account. Quite importantly, the experimental gaps separating the flat band manifold from the higher and lower energy bands cannot be reproduced if only in-plane or only out-of-plane relaxations are considered. The stability of the relaxed bilayer at the first magic angle is estimated to be of the order of 0.5-0.9 meV per atom (or 7-10 K). Our calculations shed light on the importance of an accurate description of the vdW interaction and of the resulting atomic relaxation to envisage the electronic structure of this really peculiar kind of vdW bilayers.


I. INTRODUCTION
After the first experimental findings 1,2 , Twisted Bilayer Graphene (TBG) has been subject of intense investigation from both the experimental and theoretical point of view.
Most of the unconventional transport properties of TBG originate from the almost flat bands (FBs) at the Fermi energy, originally predicted in Ref. 9 , whose bandwidth, of the order of ∼ 10 meV, has been confirmed also from tunnel spectroscopy experiments 8,10-12 . The FBs manifold, which can host up to four electrons above the Fermi energy and four holes below it, is separated by an energy gap of ∼ 50 meV from both higher and lower energy bands, and has been clearly observed in recent nano-ARPES measurements 13 . When an external gate tunes the system chemical potential within these gaps, a clear band insulating phase appears. A second, unexpected, insulating phase shows up at half-filling of the FB manifold, both on the electron and on the hole side (±2 electrons with respect to charge neutrality). The correlated insulating phase is attributed to enhanced electronelectron interaction within the FBs respectively 14,15 , al-though some authors are highlighting the relevance of the electron-phonon interaction [16][17][18][19] . After electrostatic doping, achieved by gating the structure, unconventional superconductivity, with a critical temperature ranging from 1.7 to 3 K appears in a strong pairing regime, with a phase diagram very similar to that of the underdoped cuprates, whose origin is still to be understood 17,20 .
Similar physics is being addressed also in twisted bilayers made out of transition metal dichalcogenides 21 , germanium selenide 19 other two dimensional materials 22,23 .
That reveals how the twist angle can be used as a further degree of freedom 24 for combining two-dimensional (2D) materials to implement desired properties [25][26][27][28] . The twisted lattice geometry gives rise to topological properties of TBG [29][30][31] , unlike conventional topological materials 32 , where topological properties are mostly due to spin-orbit interactions 33,34 and Brillouin zone topology.
In this paper we apply large-scale density functional calculations to better elucidate the origin of the FBs in the single-particle band structure and show the fundamental role played by the atomic relaxation. Relaxation mechanisms have been recently addressed using semiclassical techniques 35,36 , by contrast in this manuscript we resort to a DFT approach already presented in Ref. 37 and find new relaxation patters both at the magic angle and in other low angle twisted structures. The properties of the FBs manifold, at the first magic angle is connected arXiv:2004.14323v1 [cond-mat.mtrl-sci] 29 Apr 2020 with the atomic displacements originating from the interlayer van der Waals interaction. We show that the energy gain, induced by the relaxation, becomes of the order of 10 K, much larger than the typical temperature at which unconventional superconductivity or the correlated insulating phase are observed in TBG ( T <1÷2 K). The smaller is the twist angle, the more pronounced are the atomic displacements with respect to the flat bilayer. In particular, we single out an oscillating displacement pattern of the out-of-plane displacements at smaller angles, that is smoothed at the larger angles, and a vortex-like in-plane displacement pattern, where the atoms "rotate" in opposite directions in the two planes.
Tight binding calculations both at the relaxed and unrelaxed positions are also carried out, to provide a further and less expensive tool to reproduce, especially at the first magic angle, the electronic structure. We also give the effective parameters that best approximate the ab-initio band structure within the low energy continuum theory 9 generalized in the presence of atomic relaxation 38 . Interestingly enough, they are largely independent of the twist angle, which makes the continuum model an excellent tool to describe the low energy physics at small twist angles.
The paper is organized as follows: in Sec. II we outline the technical details of the calculations. In Sec. III we extensively discuss the results on the geometrical relaxation and displacement patterns. In Sec. IV the tight binding approach and the continuum model outcomes are compared with the ab initio band structure. Finally, in Sec. V we summarize our findings and draw our final conclusions.

II. METHODS
DFT calculations have been carried out using the Vienna Ab initio Simulation Package (VASP) 39 . The vdW-DF2 exchange-correlation functional 40 has been adopted to properly take into account the long-range interactions taking place between atoms belonging to different graphene layers. A PAW pseudopotential 41,42 has been employed for carbon with the 2p orbitals in valence, and the 1s orbitals frozen in the core. The single particle Bloch waves were expanded with a plane wave basis set, using a cutoff energy of 400 eV.
TBG systems corresponding to four twist angles have been considered: θ = 1.08 • , θ = 1.61 • , θ = 2.65 • , and θ = 3.698 • . The rotation is carried out starting from two perfectly AA stacked graphene layers and rotating around an axis orthogonal to the layers and passing through two C atoms, one on top of the other, belonging to the two planes (that therefore preserve their initial AA stacking). The respective structures can also be classified, according to the notation commonly used in the literature 43 , using the pair of indexes (n, m): (31,30), (21,20), (13,12) and (9,8) Sampling of the Brillouin zone (BZ) for the selfconsistent (SCF) calculations was restricted at the Γ point for all four systems, no significant changes were observed after increasing the size of the sampling of the BZ for the smaller supercells. Single-particle energies at other points in the BZ were obtained by non-SCF calculations.
For the smallest angle, because of the size of the simulation cell, we could only compute one k-point at a time, and the reported single-particle energies were therefore referred to the Fermi energy computed in the SCF calculation. The size of the supercell in the direction orthogonal to the layers (z-axis) was initially fixed at 10Å, corresponding to about 6.5Å vacuum space, introduced to prevent periodic replicas of the TBG supercell from interacting with each other. Full relaxation of the atomic positions was carried out until the residual forces were smaller than 0.002 eV/Å. Additional calculations were repeated using supercells with z-axis of 12Å and 14Å. A small residual (maximum) relaxation of less than 0.002Å was observed as the z-axis was increased to 12Å, but no further relaxation was detectable with the largest 14Å vacuum space. All symmetries were turned off. Further detail on the calculations can be found in our previous paper 37 .
Tight-binding calculations of the TBG electronic structure at different twist angles and geometries were also carried out using the Slater-Koster tight binding parametrization for p z carbon atoms: Here r 0 = 0.184a is the decay length of the transfer integral, a 0 = a/ √ 3 is the first-neighbor distance in graphene, d 0 = 0.335 nm is the intralayer distance, chosen in agreement with that of graphite. V 0 ppπ = −2.7eV and V 0 ppσ = 0.48eV are the in-plane and-out-of plane nearest-neighbours hopping energy as from Ref. 44 .
We also adopt a continuum model generalizing the model proposed in Ref.s 9,38,44,45 , providing an effective low-energy band structure. Within this approach, the two planes are coupled via two overlap coefficients u, u , that can be expressed as integrals involving the tightbinding hopping term t(R) of Eq.1. The special case u = u corresponds to the unrelaxed graphene bilayer. In order to give a minimal model capable of describing (at least) the low energy properties of the ab-initio band structure, we do not calculate u, u but use them as fitting parameters. In the following we will show that the fitted parameters are relatively close to (but quantitatively different from) those obtained performing the hopping integrals 37 . Remarkably, we will show that u, u can be chosen almost independently of the twist angle.

III. GEOMETRIC RELAXATION
Through a proper inclusion, within the ab initio approach, of the long-range inter-layer vdW interaction, we provide a detailed and accurate description of the atomic relaxations arising from the inter-layer interaction. We start discussing the out-of-plane atomic displacements. We give two complementary representations, in Fig. 1 and Fig. 2. The former shows a color/relief map of the atomic displacements with respect to the ideal, unrelaxed structure. To properly understand the results, we should recall that the TBG is built up starting from an ideal AA stacked bilayer, and then rotating around an axis orthogonal to the graphene planes and passing through to two atoms, each belonging to a different layer. After twisting, these latter atoms preserve the original AA stacking, that The displacements are measured as the difference ∆z = z − zavg of the z coordinate of each atom in a given plane and the average z coordinate zavg in that plane (∆z = 0 for all atoms in the initial, twisted but unrelaxed configuration). In each box the higher (lower) panels correspond to the top (bottom) plane, whereas the left (right) panels correspond to the long (short) diagonal (notice the different scale on the horizontal axis). s represents the coordinate along the two diagonals, with s = 0 corresponding in both cases to the unit cell origin (or a lattice equivalent site), where the AA stacking is preserved. To draw this plot the atoms whose projections in the x − y plane lie onto are are closest to the unit cell diagonals have been considered.
is kept also after relaxation. However, due to the steric repulsion of their p z orbitals, they move far apart from each other. This corresponds to the hills (valleys) in the top (bottom) plane, clearly visible in the left (right) panels of Fig. 1, highlighted in red. This can also be easier inferred from Fig. 2. Here, for the different twist angles, we show the out-of-plane displacements of the atoms positioned onto or closest to the unit cell long and short diagonals. Because the starting configuration is that of two twisted but flat graphene planes, these displacements are visualized, in each graphene plane, as a difference of the z coordinate of each atom and the average z in that plane (highlighted with a thin solid line). We can clearly distinguish how the two atoms (one for each plane) on the cell corners preserve their initial AA stacking, with their final distance estimated to be 3.58Å. Correspondingly, we can estimate the distance d between the two graphene planes, as the difference between the averages of the z coordinates in each plane. This is reported in Table I(a) and shown in Fig. 4, where we can observe a reduced distance at small twist angles. From the figure we can also see that the interplane distance lies between the cal-culate interplane distances for the untwisted bilayer at AA and AB stacking (highlighted bi the red solid lines).
The patterns in Fig. 1 show that AA stacking regions, where there is an enhanced distance between an atom in the top plane and the corresponding (closest) atom in the bottom one, alternate with regions with AB stacking, the latter being predominant. One would expect a smooth change of the displacements when moving from AA to AB regions. This naive prediction is verified only for the larger angles. For the smaller twist angles such smooth behavior is replaced by an oscillating pattern resulting in an atomic corrugation, as clearly visible from the relief maps in Fig.s 1(a-b).
The out-of-plane displacements are accompanied with significant in-plane relaxations. The (x,y) displacement field is shown in Fig. 3. Here the color map and the vector lengths are proportional to the displacement with respect to the unrelaxed twisted bilayer. Since the patterns look quite similar at different twist angles only the result for 1.08 • is shown, whereas the color bars distinguish the different systems. Interestingly enough, a vortex-like displacement field shows up in each plane, with the vorticity The color map and vector field of the in-plane displacements in TGB for the top (left) and bottom (right) plane. Each vector points in the displacement direction and its magnitude is proportional to the displacement norm. The latter is also highlighted by the color. Since the patterns look quite similar at different twist angles only the result for 1.08 • is shown. The color bars report the measure of the displacement in the x − y plane with respect to the unrelaxed structure (inÅ units).
(intended as the curl of the displacement field) changing sign when moving from the top to the bottom plane and viceversa. Such result can be explained by considering that AB stacking regions minimize the total energy of the system. Hence, close to an AA stacking configuration, the atoms of the two layers tend to move in their plane in opposite directions, in order to minimize the overlap of their orbitals that is maximum at AA stacking. Indeed, it is observed that: i) the in-plane displacement is exactly null for the AA stacked pair of atoms (no arrow and blue region at the unit cell corners in the figure); ii) the maximum displacements are observed around the unit cell corners, thus for the atoms that mostly feel the "repulsion" due to a stacking that is quite close to AA; iii) no  Figs. 1 and 3 for TBG at the four considered twist angles. Natoms is the number of atoms in each system, ∆Erelax the difference between the energy of TBG at the optimized geometry and that of the same bilayer in the initial configuration (twisted but unrelaxed graphene bilayer). The optimized distance d between the two graphene planes, computed as the difference between the averages of the z coordinates in the two planes, is also reported. displacement is observed in the AB stacking regions, as it can be evinced by the blue, hexagonal regions. While the displacement pattern looks quite similar for all the twist angles (at variance with the out-of-plane displacements), the magnitude of the displacements decreases by an order of magnitude when moving from the smallest to the largest twist angle (the maximum displacement being of the order of ∼ 0.1Å and ∼ 0.01Å in the two cases, respectively). The relaxation energy, is the energy gained by the structure when it is allowed to relax with respect to the unrelaxed configuration, in which the atoms are arranged in ideal honeycomb lattices in two parallel planes. It is reported in Table I. Such gain is extensive, i.e. it depends on the number of atoms in the unitary cells. In order to allow for a fair comparison between different twist angles, in Table I we also report the energy gain peratom. The relaxation energy normalized to the number of atoms in the unit cell, is plotted in Fig. 4(a) as a function of the twist angle (blue dots). In addition to the four angles in Tab.I, we show also points relative to the angles 4.41 • , 5.09 • , 6.01 • , 7.34 • , 9.43 • , 13.17 • (corresponding to (n, m) = (8, 7), (7, 6), (6,5), (5,4), (4,3), (3,2) respectively), to better represent the limit of large twist angles (small unit cells).
The per-atom relaxation energy increases from 3.7 K to 10.2 K when passing from the largest twist angle (corresponding to the smallest unit cell) to the smallest twist angle (corresponding to the largest unit cell). Such relaxation energy is an estimate of the upper bound of the thermal energy that in a real experiment would induce thermal fluctuations of the atomic positions that in turn would destroy the ground-state geometry pattern shown in Figs. 1 and 3. In the case of the 1.08 • magic angle it has to be compared, for instance, with the critical temperature at which zero-resistance states are observed (T c ≤ 1.7K 2 ) or at which the correlated insulator behavior at half-filling is experimentally observed (T < 1.0K 1 ). Interestingly, the relaxation energy increases upon decreasing the twist angle. This is explained by considering  4. (a) The per-atom relaxation energy as a function of the twist angle (blue circles). The red star (diamond) refers to θ = 1.08 • and show the relaxation energy when only the xy (z) coordinates are allowed to relax. (b) The interplane distance as a function of the twist angle calculated as the difference between the averages of the z coordinates in each plane. The same distance calculated for the untwisted bilayer, in both the AA and AB stacking is also highligted by the red lines.
that in the case of large unit cells the atoms have more freedom to relax to lower total energy. In Fig. 4 b) we also show the interplane distance (averaged over the unit cell). It is almost half way between the distance of the AA and the AB bilayer. However, as we have aready noticed before, relaxation mechanisms tend to enlarge the "effectively" AB regions, hence the average distance is slighly closer to that of the AB stacking.
A. Focus on the first magic angle at θ = 1.08 • We now focus on TBG at the first magic angle, θ = 1.08 • , which shows the most intriguing and marked relaxation pattern. In the following, we investigate the interplay between the out-of-plane and in-plane relaxations. This is done performing two optimizations, the first by allowing all atoms to relax only along the z direction, the second by allowing them to relax only their graphene planes (hence the z coordinates are fixed while xy coor-dinate are free to move). These partial relaxations have nontrivial effects on the electronic properties of the TBG.
First of all, we notice that these partial relaxations lead to a smaller gain in the total energy with respect to a full relaxation. This is shown by the red star and diamond in Fig. 4. This result is not surprising and is due to the reduced freedom of the atoms (constrained along z or xy) to find optimal minimal energy configurations. Hence, neither the only-xy, nor the only-z relaxations catch most of the relaxation energy (which for this system amount to ∼ 10 K), that is due to an interplay between them. Such an interplay will result even more evident in the electronic properties, when we will discuss the energy gaps separating the FB manifold from the lower and higher energy bands.

IV. ELECTRONIC PROPERTIES
In this section we discuss the electronic structure for the four angles analyzed, with and without including relaxation mechanisms. The results are summarized in Fig.s 5, in which we plot the band structure at the angles θ = 3.698 • (panels a), θ = 2.65 • (panels b), θ = 1.61 • (panels c), θ = 1.08 • (panels d-e). The left panels show the band structures including relaxation mechanisms, while in the right panels the atoms are fixed in their lattice positions. In Fig. 5(e) we show a zoom close to the Fermi energy in the case of the magic angle θ = 1.08 • . The blue points correspond to the bands calculated using the DFT approach described in Sec. II. The red (full) lines are obtained within a tight binding calculation where, however, the atomic positions for the relaxed structures are the ones optimized within the DFT approach. In general, we notice that there is an excellent agreement between the red (full) lines and the blue points. The DFT calculation could be well approximated also adopting a continuum model (see Sec. II), where the interplane hoppings are parametrized, for the optimized structure, by coefficients u = 0.078 ± 0.002, u = 0.098 ± 0.004 and for unrelaxed ones by u = u = 0.107 ± 0.004 almost independently of the angle. That can be adopted as a minimal single particle description within more complex many body 11,38 approaches in order to study the superconductivity and the Mott insulating state of this system.
Inspection of Fig. 5 shows that, comparing the TB and the DFT curves, there is a tiny mismatch that can be evidenced, for instance, at Γ. This is due to the fact that, despite sharing the same atomic position, the two approaches do not account for e-e interaction in the same way. In particular, the tight binding approach does not takes into account any e-e interaction while in the DFT approximation the e-e is accounted within the local density approximation. This discrepancy, at the Γ point, could indeed be an estimate of the Hartree energy in TBG.
In all panels, it is clearly visible that relaxation mech-anisms tend to maximize the energy gaps at Γ and, in particular, in the case of the magic angle θ = 1.08 • the gap separating the FBs from the higher (lower) energy bands is about 26 meV (16 meV), consistent with the experiments 1,2,8,10-13 . On the other hand, those gap cannot be reproduced at all if no relaxation is allowed. Aimed at identifying if there is a predominant role of the in-plane or of the out-of-plane displacements, these energy gaps have been estimated also when the system is allowed to relax only either in the x − y plane or along the z axis: in the former case, we obtained gaps of ∼ 2 and ∼ 14 meV, that underestimate these values in the fully relaxed system, especially in the hole side; in the latter case, the two gaps turn out to be both zero, showing that a quite important role is played by the x − y displacements.
As far as the FB dispersion is concerned, we obtain a full bandwidth of ∼ 20 meV in the relaxed system, which is of the order of the one measured in the experiment of Ref. 2.

V. CONCLUSIONS
The FBs in TBG at the first magic angle can be intimately related with the atomic displacements arising as effect of the interlayer vdW interaction. Large scale first principles calculations allow us to conclude that experimental gaps cannot be reproduced if we consider a flat bilayer system. Relaxation effects are thus crucial as already noticed in our previous manuscript 37 . In this manuscript, we have also investigated partial relaxation processes with only in-plane or out-of-plane displacements, however the resulting band structure does not show the expected gaps in both cases. Out-of-plane relaxations are characterized by a strongly oscillating pattern, that is smoothed (until it disappears) at large twist angles. On the other hand, the in-plane displacements show a vortex-like configuration, where the vorticity assumes opposite values in the two planes. However, the magnitude of these displacements decreases upon increasing the twist angle. Overall, the combination of the two patterns allow the atoms to override the steric repulsion felt by the p z orbitals, maximizing the regions with AB stacking at expenses of the regions showing AA stacking.
The energy gain induced by the relaxation is the larger, the smaller is the twist angle, decreases by increasing the twist angle and eventually seems to reach a plateau of ∼6.8 K at large twist angles. The smaller angles correspond to larger unit cells, that can easier accommodate the atomic rearrangement, and correspond to large energy gain, of the order of ∼10 K, much larger than the temperature at which the most exotic phenomena, such as correlated insulating phase and superconductivity are detected. Since the latter are intimately related to the presence of the FBs, that in turn we demonstrate being related with the atomic relaxation, such temperature (10 K) should be considered as an upper limit. Higher tem- The left (right) panels corresponds to the relaxed (initial) structure. Blue dashed lines and filled dots correspond to the DFT calculation, red solid lines to the TB approach. Zero energy corresponds to the Fermi energy. The effective continuum model mentioned in Sec. II can be also used to approximate the DFT calculations. The parameters u, u are independent of the angle and, within a reasonable error, can be approximated as u = u = 0.107 ± 0.004 in the unrelaxed case and u = 0.078 ± 0.002, u = 0.098 ± 0.004 for relaxed structures.
peratures would destroy relaxation effects due to thermal atomic oscillations, and as a consequence FBs effects would be hindered.
Including relaxation effects not only we reproduce band gaps consistent with the ones measured in experiments, but we also give the effective parameters of a low energy continuum model to be adopted for further investigation including correlation effects. Interestingly enough, despite the fact that relaxation patterns have different shapes at various angles, the interplane hopping coefficients u, u are found to be almost independent of the twist angle. That could be a useful hint to apply the continuum model also to smaller angles, where the unit cell would become umpractically large to be attacked with atomistic approaches.