First-principles calculation of gate-tunable ferromagnetism in magic-angle twisted bilayer graphene under pressure

The TBG system we choose to study is labeled by integers (12,11) with a corresponding twisting angle of θ = 2.88^∘. Both for simplicity and because electronic properties of TBG are mostly smooth functions of θ but not of (n, m), we will refer to the system by its twisting angle. The number of atoms per moiré unit cell is 1588, which is the local minimum among those commensurate TBGs that have twisting angles close to 2.88^∘ and pushes to the limit of our computational capability. According to Carr et al [16], the pressure that is required to induce flat bands near the Fermi level in the θ = 2.88^∘ TBG is experimentally accessible. In the following, we decrease the interlayer distance of this system from its value at zero pressure and study the evolution of the band width at the Fermi level. All the calculations in this subsection are non-spin-polarized.

Figure 1(a) shows the band structure of the θ = 2.88^∘ TBG with interlayer distance $d\,=\,3.42$ Å, which corresponds to essentially zero pressure. In this case, all the bands are dispersive, and a clear Dirac point is present at the Fermi level at K. The band that ranges from 0 to $0.28\,\textrm{eV}$ has four-fold degeneracy along the k-point path $M\textrm{–}K\textrm{–}\Gamma$, two-fold from spin degeneracy and two-fold from valley degeneracy. This is also the case for the band that ranges from $-0.24\,\textrm{eV}$ to 0. Therefore, full filling of these two sets of bands corresponds to ±4 electrons doped per moiré unit cell. Figure 1(b) shows the band structure at interlayer distance $d\,=\,2.85$ Å, which from our calculation corresponds to a pressure of about $14\,\textrm{ GPa}$. To calculate this pressure, we sum up the out-of-plane component of force on each atom in one graphene layer and divide the result by the area of the moiré unit cell. At this separation, flat bands appear near the Fermi level. In this case, the single particle gap between the flat bands and the dispersive ones above remains, while the flat bands below become almost degenerate with the dispersive bands at the Γ point. ${\Delta}E_{\Gamma}$ versus the interlayer distance d is plotted in figure 1(c), where $\Delta E_{\Gamma}$ is the total width of the bands at the Fermi energy as measured at the Γ point (see figure 1(a)); for instance, $\Delta E_{\Gamma}\,=\,525\textrm{ meV}$ at $d\,=\,3.42$ Å. The width ${\Delta}E_{\Gamma}$ decreases as d is reduced until 2.85 Å, where it vanishes; $\Delta E_{\Gamma}$ becomes finite again for even smaller d, which is consistent with previous study [16]. Displayed in figure 1(d) is the modulus of a flat-band wavefunction at the K-point of the θ = 2.88^∘ system with interlayer distance $d\,=\,2.85$ Å. The wavefunctions of flat bands exhibit strong localization at AA stacking sites, which makes sense since from a tight-binding point of view localized states have reduced hoping integrals, and this usually indicates band flatness. Thus we have found an MATBG system of θ = 2.88^∘ and $d\,=\,2.85$ Å under pressure that will be our focus in the remaining part of the paper.

Zoom In Zoom Out Reset image size

Next, we investigate single-gate field effects on the θ = 2.88^∘ TBG with interlayer distance $d\,=\,2.85$ Å, which possesses flat bands around the Fermi energy. A single-gate setup allows electrostatic charge doping of the TBG.

Figures 2(a)–(c) display the band structures of the NonM solution of the system under doping levels of two electrons, charge neutral and two holes per moiré unit cell, respectively. Upon finite doping, the valley degeneracy of the flat bands along the k-point path $M\textrm{–}K\textrm{–}\Gamma$ is lifted due to a gate-induced electrostatic potential difference between the two graphene layers. At most doping levels between four electrons per moiré unit cell and four holes per moiré unit cell, both the flat bands and the dispersive bands around the Fermi energy are doped and thus partially occupied. The shift of the dispersive bands in energy relative to $E_\mathrm F$ due to doping is much larger than that of the flat bands. As a result, the dispersive bands cross the flat bands near the Γ point, where a single particle gap opens due to their hybridization.

Zoom In Zoom Out Reset image size

Figures 2(d)–(f) show the band structures of the FM solution of the system under the same doping levels as figures 2(a)–(c), respectively. Here both the valley and spin degeneracy of the flat bands are broken at finite doping levels, resulting in an unbalanced occupation between the spin-up and spin-down channels and thus a net magnetic moment. The inset of figure 2(e) is the spin density in real space of the FM solution of the neutral system, displaying strong localization at AA stacking sites. This is because spin splitting mainly occurs in the flat bands, which are by themselves localized at these sites. Thus the spin density, defined as the difference between spin-up and spin-down electron densities, can only be localized also at AA stacking sites. This localization persists through all doping ranges from four electrons to four holes doped per moiré unit cell.

Finally, in figure 2(g), we show the energy difference between the FM and NonM solutions and the magnetization per moiré unit cell of the FM solution at different doping levels. The total energy of the FM solution is lower than that of the corresponding NonM one at doping levels ranging from two electrons to four holes doped per moiré unit cell. In the case of zero doping (see figures 2(b) and (e)), this can be understood by examining the band energy: the spin-up (spin-down) flat bands are lowered (raised) in energy and become fully occupied (nearly empty) after spin splitting. The FM state is most stabilized at the doping level of two holes per moiré unit cell. In contrast, the NonM solution has lower energy at the doping level of three electrons per moiré unit cell. If we compare the magnetic moment of the FM solution at the doping of three electrons with that at the doping of three holes, we find that the former (0.27 $\mu_\mathrm B$) is significantly lower than the latter (1.27 $\mu_\mathrm B$). A low magnetic moment implies a small spin splitting, as shown by the band structures in the inset of figure 2(g), in dramatic contrast with the significant spin splittings in figures 2(d)–(f). This is correlated with the fact that the nonmagnetic state has lower energy at the doping level of three electrons per moiré unit cell. The magnetization is highest near zero doping and decreases monotonically as the doping level increases in the directions of both electron and hole doping. It vanishes at both four electrons and four holes doped per moiré unit cell. Both electron and hole doping reduce charge imbalance between the two spin channels as well as the splitting between spin-up and spin-down flat bands. For TBG under zero pressure, the FM state was found at a doping level of three electrons per moiré unit cell in experiments [8]. This doping level is different from our calculated doping level at which the FM state is stabilized the most. Such a difference is likely due to a change in pressure. However, the physics of the modulation of FM-state-inducing doping level by pressure is still unclear.

We remark here that we have also calculated θ = 2.88^∘ TBGs with interlayer distances $d\,=\,2.83$ Å and 2.87 Å at doping levels of ±2 charges per moiré unit cell. An FM solution is not found for the $d\,=\,2.87$ Å systems, while in the $d\,=\,2.83$ Å cases, an FM solution with magnetization M per moiré unit cell of only $0.5 \,{\mu}_\textrm{B}$ ($0.3\,{\mu}_\textrm{B}$) exists at a doping level of two holes (electrons) per moiré unit cell. Comparing these to the magnetizations of the $d\,=\,2.85$ Å systems at the same two hole and two electron doping, which both have $M\,=\,2.2\,\;{\mu}_\textrm{B}$, we conclude that the FM solution only exists in a narrow range of interlayer distance around the flat-band condition.

Yndurain [33] suggested creating magnetic structures in compressed TBG with a twisting angle larger than 1.1^∘. In Yndurain's work, logarithmic van Hove singularities shift toward the Dirac point and eventually merge into a single peak as two graphene layers were brought closer to each other. He reported that a TBG with θ = 5.08^∘ becomes magnetic with a magnetic moment of nearly $4\;\,\mu_\mathrm B$ per moiré unit cell at a pressure of $2.19\,\textrm{ GPa}$. Our results are consistent with Yndurain's findings in this regard. Lopez-Bezanilla [34] studied substitutional chemical doping of θ = 5.09^∘ TBG under pressure and found that N-doping can turn FM order into antiferromagnetic order. Compared to chemical doping, electrostatic doping provides continuous and reversible control over the doping level without destroying intrinsic properties of the system under study. However, limited by the large size of our θ = 2.88^∘ system, we include only one moiré unit cell in our simulations, and thus it remains an open question whether antiferromagnetic order occurs upon electrostatic doping.

Atomic structure relaxation in ab initio calculations is known to open and widen a single particle gap below as well as above the flat bands [16, 35]. This gap in the case of zero pressure MATBG is known to be less than ∼50 meV, both from ab initio calculations [16, 35] and from experiments [2, 36], and its size has been shown to decrease if the twisting angle is increased from 1.1^∘. However, the relative shift in energy between the dispersive and flat bands due to charge doping varies over a much larger scale in our calculation of the θ = 2.88^∘ and $d\,=\,2.85$ Å system. This can seen from how the energy, relative to the Fermi level, of the band maximum (minimum) of the dispersive bands, which are entirely below (above) the flat bands at zero doping, evolves with the doping level, as shown in figure 3 for the NonM solutions. Since the energy of the flat bands is not sensitive to the doping level within the range of ±4 charges doped per moiré unit cell (see figures 2(a)–(c)), the band edges plotted here are a good measure of the hybridization between the dispersive bands and the flat bands. The maximum (minimum) of the dispersive bands below (above) the flat bands at zero doping moves up (down) in energy with respect to the flat bands by ∼150 (∼200) meV, when the doping level varies from four electrons (holes) to four holes (electrons) doped per moiré unit cell. Although a full atomic structure relaxation under pressure is beyond the scope of this work, the large size of the relative shift in energy between the dispersive and flat bands shown here indicates that the common assumption that the flat bands are perfectly isolated from the surrounding dispersive ones may not hold true for MATBG at large twisting angles at certain finite doping levels.

Zoom In Zoom Out Reset image size

A dual-gate setup permits an out-of-plane electric field between two gate electrodes as well as charge doping. In this subsection, we examine the effect of out-of-plane electric field on the charge neutral θ = 2.88^∘ TBG with interlayer distance $d\,=\,2.85$ Å. Since our dual-gate setup is symmetric, a negative out-of-plane electric field is equivalent to an oppositely directed electric field with the same magnitude, and it is sufficient to consider only positive electric fields applied to the TBG. Figure 4 shows the magnetization per moiré unit cell M of the FM solution versus the strength of the out-of-plane electric field. M decreases with the electric field from $4.0\,{\mu}_\textrm{B}$ at zero field. This decrease is slow and only $0.17\,{\mu}_\textrm{B}$ ($4.2\%$) for an electric field up to 0.6 V Å⁻¹, which is already a large field in experiments. Figure 4 also shows the energy difference between the FM state and the NonM state as a function of electric field. Notably, a phase transition from the FM state to the nonmagnetic state takes place between 0.4 V Å⁻¹ and 0.6 V Å⁻¹.

Zoom In Zoom Out Reset image size

The reason for such a phase transition can be understood in the band picture in terms of a competition between the spin splitting energy of the FM state and the energy splitting in the valley degree of freedom. Figure 5(d) shows the band structure of the FM solution at zero electric field, where the spin-up flat bands are lowered in energy to be fully occupied while the spin-down flat bands become nearly empty. After such a spin splitting, the occupied flat-band states have lower energy than those of the nonmagnetic solution at zero electric field as shown in figure 5(a). This is the major reason the FM state is energetically favorable at small out-of-plane electric fields. As the strength of the out-of-plane field increases, the valley splitting of flat bands (along the k-point path $M\textrm{–}K\textrm{–}\Gamma$) of the nonmagnetic solution increases monotonically, which can be seen from figures 5(a)–(c) for the out-of-plane electric fields of 0, 0.4 and 0.8 V Å⁻¹ respectively. Consequently, the energies relative to the Fermi level of the occupied states are mostly lowered, which contributes to a decrease in the total energy of the NonM state $E_\textrm{NonM}$. In contrast, the spin splitting of flat bands of the FM solution is insensitive to the out-of-plane electric field, although the valley splitting for each spin channel is comparable to that of the corresponding nonmagnetic solution (see figures 5(d)–(f)). It is also noteworthy that the spin splitting of flat bands of the FM solution is not sensitive to either the valley band index or to the crystal momentum. Overall, the energy of the FM state $E_\textrm{FM}$ does not vary as much as that of the NonM state with out-of-plane electric field, as long as the spin-up (spin-down) flat bands remain fully occupied (empty). Therefore, the energy difference $E_\textrm{FM}\,-\,E_\textrm{NonM}$ increases with the electric field below 0.4 V Å⁻¹. At an out-of-plane electric field $\geqslant$0.4 V Å⁻¹, the valley splitting is so large that the otherwise-empty spin-down flat bands become partially occupied near the Fermi level along $M\textrm{–}K$ cuts in the BZ, raising the energy of the FM state. Consequently, the energy difference $E_\textrm{FM}\,-\,E_\textrm{NonM}$ increases greatly between 0.4 V Å⁻¹ and 0.6 V Å⁻¹, accompanied by an FM-to-NonM phase transition.

Zoom In Zoom Out Reset image size

Next, we examine how impurities affect electronic properties of the TBG with θ = 2.88^∘ and $d\,=\,2.85$ Å. We consider two impurities separately, namely an $\textrm{O}_2$ molecule and a Stone–Wales defect. Both impurities are located at the AA site in our calculations. The oxygen molecule is chosen because it frequently occurs in 2D material experiments [37]. The Stone–Wales defect is chosen since it is the lowest-energy defect in graphene [38].

We determine the adsorption configuration of oxygen under pressure in the following way. Mimicking the experimental setup, we sandwich the oxygen molecule between bilayer graphene and bilayer boron nitride (BN) [2]. The computational load would be prohibitively heavy if we adopt the moiré unit cell. For computational feasibility, we choose a smaller cell which is as large as a $5\,\times\,5$ graphene supercell. During structural relaxation, the bilayer graphene is fixed with an interlayer distance of 2.85 Å. Both BN layers are also kept flat, with the outer BN layer fixed and the inner BN layer, which is in contact with the oxygen molecule, free to move in the out-of-plane direction. We tune the vertical position of the outer BN layer until the pressure on it equals the pressure on the outer graphene layer. We find the most stable adsorption configurations by trying several randomized initial configurations.

The relaxed $\textrm{O}_2$ molecule is then used on TBG for further calculations. Interestingly, we find that ferromagnetism is completely suppressed by the molecule with all different configurations.

Figure 6(a) shows in red the typical band structure of TBG with an oxygen molecule per moiré unit cell located at the locally AA stacking site. With the oxygen molecule, the ground state is now non-magnetic and no FM solution is found. Each band in red shown in figure 6(a) has a spin degeneracy. The energy bands of TBG in the nonmagnetic state without oxygen are also shown for comparison in the same figure. Comparing the two sets of bands, an energy gap of $29\textrm{ meV}$ at the Γ point and larger elsewhere opens around the Fermi level with the adsorption of $\textrm{O}_2$. The spin degenerate flat bands, now two above the Fermi level and two below are still filled with four electrons per moiré unit cell after the absorption of $\textrm{O}_2$. This indicates that no electrons are transferred from an oxygen molecule to the TBG. Compared with the band structure of the ground state of the pure system at zero doping (see figure 2(e)), where spin majority channel of the flat bands has nearly four electrons filled while the spin minority channel has nearly zero, the absorption of the $\textrm{O}_2$ molecule causes a transfer of two electrons from the spin-majority channel to the spin-minority channel, such that the whole system becomes non-magnetic. Figures 6(b) and (c) shows the PDOS of all C (O) atoms, where the PDOS is calculated at the Γ point in reciprocal space. The energy bands near the Fermi level, which is set to zero, originate mainly from TBG. Minor contributions from oxygen atoms are seen around $\pm 7\textrm{ meV}$ and $\pm 10\textrm{ meV}$, indicating hybridization between TBG and $\textrm{O}_2$.

Zoom In Zoom Out Reset image size

We obtain the atomic structure of a Stone–Wales defect based on AA-stacking bilayer graphene in a $7\,\times\,7$ graphene supercell. All atoms are free to relax in the plane, and the interlayer distance is fixed at 2.85 Å. The relaxed Stone–Wales defect is used in TBG for further calculations. Figure 7 shows the spin-polarized band structure. Again, we find that ferromagnetism is suppressed due to charge transfer from the spin-up channel to the spin-down channel.

Zoom In Zoom Out Reset image size

In this subsection, we examine the bonding character in TBG under pressure via Wannier functions, which have been applied to obtain accurate tight-binding models for both monolayer graphene [39] and AB stacked bilayer graphene [40]. Since spin-polarization makes little difference in the shape of the flat energy bands, we performed the Wannier analysis for non-spin-polarized TBG. The left panels of figure 8(a) (figures 8(b) and (c)) show both the DFT band structure and the interpolated band structure based on a 4-band (20-band) model Hamiltonian in the Wannier function representation. As can be seen from the figure, they match well for both the 4-band model and the 20-band model. The right panels of figure 8 show the PDOS projected onto each Wannier function. In the 4-band model, all four Wannier functions contribute almost equally to each flat band. In the 20-band model, the twenty Wannier functions are of two types; the first (type 1) consists of Wannier functions of indices ranging from 1 to 12; and the second (type 2) consists of those of indices ranging from 13 to 20. Figure 8(b) shows that both Wannier function types contribute significantly to the dispersive bands, with more contribution from type 2 than that from type 1. Figure 8(c) shows that only Wannier functions of type 1 contribute to the flat bands; the contribution from any Wannier function of type 2 is essentially zero. In figure 9(a), we show Wannier Hamiltonian matrix elements $h_{mn}(\boldsymbol{R})\,=\,\langle \boldsymbol{0} \, m \vert H \vert \boldsymbol{R} \, n \rangle$ for the 20-band model, where 0 is the origin of the home moiré unit cell, R is the origin of any moiré unit cell, and m and n are indices of Wannier functions in the moiré unit cells at 0 and R respectively. The Hamiltonian matrix elements can be significant even if two Wannier functions are separated by three lattice constants of the moiré unit cell. Beyond three lattice constants, the corresponding matrix elements become less important. We verify this point by setting those matrix elements to zero and still obtain flat bands that are close to the DFT results. The $20\,\times\,20$ matrix $h_{mn}(\boldsymbol{R})$ for any R can be rearranged into a (nearly) block diagonal form (see figure 9(b)) by reindexing the Wannier functions as follows: $1~\rightarrow~1$, $4~\rightarrow~2$, $5~\rightarrow~3$, $8~\rightarrow~4$, $9~\rightarrow~5$, $12~\rightarrow~6$, $13~\rightarrow~7$, $14~\rightarrow~8$, $19~\rightarrow~9$, $20~\rightarrow~10$, $2~\rightarrow~11$, $3~\rightarrow~12$, $6~\rightarrow~13$, $7~\rightarrow~14$, $10~\rightarrow~15$, $11~\rightarrow~16$, $15~\rightarrow~17$, $16~\rightarrow~18$, $17~\rightarrow~19$, and $18~\rightarrow~20$. Major features of both flat and dispersive energy bands are retained when we set the off-block-diagonal matrix elements to zero. If we keep only the Hamiltonian matrix elements due to the first ten (in the new index) Wannier functions, we obtain two flat bands, four dispersive bands below the Fermi level, and four dispersive bands above the Fermi level. Similarly, we can obtain the remaining ten bands using the remaining ten Wannier functions only. The block diagonal form allows us to classify the Wannier functions into two separate sets that are responsible for two disjoint sets of energy bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10(a) (figure 10(b)) shows the Wannier charge centers (WCCs) for the 4-band (20-band) model. In the 4-band model, from a top view all four WCCs are located almost at AA stacking sites. The flat-band states are localized at AA sites, which is consistent with previous DFT calculations [18] and STM measurements [41–43]. From a side view, all four WCCs are in the middle between the two graphene layers. This indicates that each flat band mediates a certain bonding between the two graphene layers (via the AA site). In the 20-band model, from a side view all twenty WCCs are in the middle of the two graphene layers as well. However, the in-plane positions of the WCCs are quite different from those in the 4-band model: type 1 (type 2) WCCs are closer to AA sites (AB sites). Specifically, each type 1 (type 2) WCC is about $0.181\,a_0$ ($0.131\,a_0$) away from an AA (AB) site, where, a₀ is the period of the moiré pattern. That the WCCs are away from AA sites in the 20-band model is likely due to the flat bands being mixed with dispersive bands. We infer that the dispersive bands mediate bonding between the two graphene layers via sites other than AA sites.

Zoom In Zoom Out Reset image size

Figure 11(a) shows the real part of the first Wannier function for the 4-band model. The isosurface level is about $2.3\%$ of the maximum value of the real part in all of space. The spread of this Wannier function is about $18.1\%$ of the area of the moiré unit cell. The spread Ω is defined as

$$\begin{equation} \Omega\,=\,\left\langle w_{n 0} \left|r^{2}\right| w_{n 0} \right\rangle\,-\,\left|\left\langle w_{n 0}|\boldsymbol{r}| w_{n 0})\right\rangle\right|^{2}, \end{equation} \tag{ 2 }$$

where $w_{n 0}$ denotes the $n\textrm{th}$ Wannier function that belongs to the home unit cell. Figure 11(a) also shows three cross sections that pass the WCC as indicated by the gray dot over the isosurface or by the black arrow near each cross section. It can be seen from the first (second) cross section that π-like bonds between neighboring carbon atoms occur near the WCC in the bottom (top) graphene layer. As such, interlayer π–π hybridization seems to be present near the WCC. If we examine a position that is far away from the WCC, we find that the intralayer π bonding between neighboring carbon atoms becomes weaker, while the signature of individual p_z orbitals becomes stronger. Figures 11(b) and (c) show the first Wannier function (type 1) and the thirteenth Wannier function (type 2) respectively for the 20-band model. The spread of the first (thirteenth) Wannier function is about $13.2\%$ ($27.4\%$) of the area of the moiré unit cell. It is reasonable that the first Wannier function has a much smaller spread than the thirteenth Wannier function, since the former originates partially from the flat bands while the latter originates mostly from the dispersive bands. The cross sections in figure 11(b) show that some π-like orbital of the top graphene layer hybridizes with some p_z -like orbital of the bottom graphene layer near the WCC. Even more pronounced p_z -like orbital character can be observed in figure 11(c). We may understand this as follows: the 20-band model involves dispersive energy bands, which possess mainly p_z orbital character, as in the case of monolayer graphene [39]. It is still true for the 20-band model that p_z orbital character becomes obvious at positions far away from the WCC.

Zoom In Zoom Out Reset image size

Lu et al [7] reported states with non-zero Chern numbers in MATBG with one electron or one hole doped per moiré unit cell and under certain magnetic field. Electrically tunable Chern number was also found in an ABC-trilayer graphene/hexagonal boron nitride moiré superlattice [44]. Inspired by these studies, we examine the band topology of the θ = 2.88^∘ TBG under pressure with interlayer distance $d\,=\,2.85$ Å. The characteristic topological invariant Chern number is defined by [45, 46]

$$\begin{equation} C\,=\,\frac{1}{2~\pi} \int_{S} \sum_{n \in B} \nabla_{\boldsymbol{k}}\,\times\,i \left\langle u_{n \mathbf{k}}\left|\nabla_{\boldsymbol{k}}\right| u_{n \boldsymbol{k}}\right\rangle \cdot \textrm{d} \boldsymbol{S}, \end{equation} \tag{ 3 }$$

where n is the band index, B is an isolated set of energy bands, $u_{n \boldsymbol{k}}$ is the periodic part of a Bloch wave, and S is a closed orientable 2D surface (the whole BZ in our case). For the sake of the flat bands being isolated from the rest energy bands, we chose the case of one electron doping per moiré unit cell (see figure 12(a)). It is noteworthy that the TBG exhibits half-metallicity in this case. In practice, the Chern number can be found by examining how the sum of hybrid Wannier charge centers (HWCCs) changes with crystal momentum [47]. Figure 12(b) shows the sum of HWCCs versus the crystal momentum k_y , which was calculated using the Z2Pack package based on a 20-band model Hamiltonian in the representation of Wannier functions [48]. Since the winding number of the curve for each spin in figure 12(b) is zero, the Chern numbers $C_\uparrow$ (for spin up) and $C_\downarrow$ (for spin down) are both zero, and the flat bands in each spin channel of the (θ = 2.88^∘, $d\,=\,2.85$ Å) TBG with one electron doped per moiré unit cell have trivial band topology according to our calculations. It is still possible that non-trivial band topology occurs for the TBG at a different doping level. Previous studies suggest that non-trivial band topology may occur in TBG if strong correlation is taken into account [7, 49].

Zoom In Zoom Out Reset image size