Stoner ferromagnetism in low-angle twisted bilayer graphene at three-quarters filling

The stacking geometry of TBGs is characterized, as standard [20, 47, 48], by the twist of one graphene layer with respect to the other around a given site starting from the AA-stacked bilayer, forming a moiré pattern. We describe a TBG system by twisting the upper layer at an angle $\theta$ as described in Ref. [20, 21]. The primitive lattice vectors of the bottom layer are $\textbf{a}^{b}_{1}=\sqrt{3}a_{0}\hat{\bf e}_{x}$ and $\textbf{a}^{b}_{2}=\sqrt{3}a_{0}/2\left(\hat{\bf e}_{x}+\sqrt{3}\hat{\bf e}_{y}\right)$ with the carbon-carbon bond length $a_{0}=1.42$ Å, and those of top layer as $\textbf{a}^{t}_{i}=R(\theta)\textbf{a}^{b}_{i}$, where $R(\theta)$ is the rotation matrix. In this study, we take the spacing between graphene layers as $d_{0}=3.35$ Å for non-relaxed lattices [19].

The lattice structure of a TBG system is commensurate if the periods of the two graphene layers match, giving a finite unit cell. Hence, the periodicity condition requires the lattice translation vector $m\textbf{a}_{1}^{b}+n\textbf{a}_{2}^{b}$ of the bottom (unrotated) layer and $n\textbf{a}_{1}^{t}+m\textbf{a}_{2}^{t}$ in the top (rotated) layer, with $m$ and $n$ integers, to coincide. The twist angle $\theta$ for a commensurate structure is related to $(m,n)$ by [49, 50, 51]

with a lattice constant

The commensurate unit cell contains $N=4(m^{2}+n^{2}+mn)$ atoms.

Due to the symmetry of the honeycomb lattice, when the translation vector of the top layer is fixed to $\boldsymbol{\delta}^{t}=0$, a twist at an angle $\theta=120^{\circ}$ transforms the AA-stacked bilayer into itself. In turn, a twist by $\theta=60^{\circ}$ transforms the bilayer from AA- to AB-stacking. TBGs with $\theta$ and $-\theta$ are mirror images that share equivalent band structures [49].

Each graphene layer is constituted by two sublattices, ${\rm A}^{b({\rm or}\;t)}$ and ${\rm B}^{b({\rm or}\;t)}$. The commensurate structures occur in two different forms distinguished by their sublattice parities [50, 51]. Figure 1(a) shows an example of an even commensurate structure under sublattice exchange. It is characterized by having 3 three-fold symmetric positions that correspond to the stacking of the ${\rm A}^{b}{\rm A}^{t}$ and ${\rm B}^{b}{\rm B}^{t}$ sites and by hexagonal centers that overlap, ${\rm H}^{b}{\rm H}^{t}$. In contrast, Fig. 1(b) shows an odd commensurate structure. Here, the top and bottom coinciding sites correspond to ${\rm A}^{b}{\rm A}^{t}$ at the origin, and the two remaining three-fold symmetric positions are occupied by ${\rm B}^{b}$(or ${\rm B}^{t}$)-sublattice sites of one layer aligned with the hexagon centers ${\rm H}^{t}$(or ${\rm H}^{b}$) of its neighbor layer. Therefore, the angles $60^{\circ}-\theta$ and $-\theta$ followed by a relative translation of the upper layer by $\boldsymbol{\delta}^{t}=(\textbf{a}_{1}^{t}+\textbf{a}_{2}^{t})/3$ form commensurate partners with unit cells of equal areas but opposite sublattice parities.

A TBG with $\theta=30^{\circ}$ is a special case, since its crystal structure is its own commensurate partner and corresponds to an elementary two-dimensional quasicrystalline lattice [52, 53, 54].

Here, we work with odd commensurate structures. In Fig. 1(b), the carbon atom sites, indicated by the green disks, are equivalent due to the three-fold symmetric positions and two-fold rotation axis (green dashed line). These symmetry properties allow us to reduce the numerical computation by identifying the equivalent sites within the primitive unit cell.

Lattice relaxations have been shown to significantly influence the electronic states of low-angle TBGs [55]. Our calculations, which incorporate both in-plane and out-of-plane lattice relaxations, are implemented as follows: For a given twist-angle $\theta$ we determine the atomic relaxation positions of the lattice structures defined above using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [56]. The energy minimization process involves both non-bonded and bonded interactions to accurately model the system. Non-bonded interlayer interactions are handled using a registry-dependent interlayer potential specifically designed for graphene/hBN heterostructures [57, 58, 59], with a cutoff distance of 16 Å, ensuring that all relevant interlayer forces are considered. The parameters for the registry-dependent potential are given in [60]. Bonded intralayer interactions are modeled using the Tersoff potential [61, 62], which is well-suited for describing the covalent bonding in carbon-based materials like graphene. The lattice structure relaxation calculations are performed with an energy tolerance criterion set to $10^{-10}$ eV per atom [63].

We model the electronic properties of low-angle TBG systems using the effective Hamiltonian

where $H_{\rm TB}$ stands for the TBG tight-binding Hamiltonian and $H_{\rm U}$ for a Hubbard on-site electron-electron interaction term [27, 64].

The tight-binding Hamiltonian reads

where the operators $c_{i\sigma}$ and $c^{\dagger}_{i\sigma}$ annihilate and create an electron with spin $\sigma=\left\{\uparrow,\downarrow\right\}$ at site $i$, respectively. $t_{ij}$ is the transfer integral between the Wannier electronic orbitals centered at the carbon sites $i$ and $j$. The transfer integral $t_{ij}\equiv t({\mathbf{R}}_{ij})$, where ${\mathbf{R}}_{ij}={\bf R}_{i}-{\bf R}_{j}$, depends on the interatomic distance and on the relative orientation between $p_{z}$ orbitals. $t_{ij}$ is parameterized as [19, 49, 65]

with

Here, $V^{0}_{pp\pi}=-2.7$ eV, $V^{0}_{pp\sigma}=0.48$ eV, and $r_{0}=0.319a_{0}$ is the decay length of the transfer integral [65, 19]. The transfer integral for $R>5$ Å is exponentially small and can be safely neglected. The intra- and interlayer hopping matrix elements $t_{ij}$ have been determined by fitting the dispersion relations of graphene monolayer and graphene AB-stacked bilayer obtained from ab initio calculations [19].

The Hubbard term reads

where $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the spin-resolved electron number operator at site $i$. The parameter $U>0$ gives the magnitude of the on-site Coulomb repulsion [64, 27]. The magnitude of $U$ in graphene systems is has been extensively discussed [66, 67] with estimates that vary from $0.5\cdots 2.0V^{0}_{pp\pi}$. In this study we consider, unless otherwise stated, $U=V^{0}_{pp\pi}$.

We solve the TBG Hamiltonian for a given filling factor in the mean-field approximation. As standard, we write $n_{i\sigma}\equiv\langle n_{i\sigma}\rangle+\delta n_{i\sigma}$ and neglect the quadratic terms in $\delta n_{i\sigma}$ that are responsible for electronic correlations. The mean field Hubbard Hamiltonian reads

The last term is a constant for a given local magnetic configuration and can be absorbed in the chemical potential.

Let us now describe how the computation of the ground state charge density and magnetic properties are implemented.

Due to the large number of atoms in a low-angle TBG primitive cell, even mean-field calculations can be computationally very costly. Several effective Hamiltonian models have been developed to reduce the numerical effort [68, 46, 69, 70]. Here, we perform a full real-space calculation using the HHK recursive method, which is very efficient to compute the single-particle LDOS of large systems. The HHK method [43, 44, 45, 42] puts forward a very efficient Lanczos-like $O({\cal N})$ recursive procedure that transforms an arbitrary sparse Hamiltonian matrix in a tridiagonal one. Next, it evaluates the diagonal Green’s function in real space by a continued fraction expansion, which is much more numerically amenable than a full diagonalization.

By a suitable choice of the chemical potential, our approach allows us to consider charge neutral as well as systems with a finite doping, namely, $n_{\rm dop}=N_{e}/A_{\rm PUC}$, where $N_{e}$ is the number of electrons in excess to the CNP and $A_{\rm PUC}=3\sqrt{3}a_{0}^{2}/8\sin^{2}(\theta/2)$ is the area of the TBG moiré cell with $N$ atoms.

We implement the self-consistent mean-field calculation as standard:

($i$) We start with an initial set of occupation numbers $\langle n_{i\sigma}\rangle$, with the constraint $\sum_{i\sigma}\langle n_{i\sigma}\rangle=N+N_{e}$, where the sum runs over all sites of the moiré unit cell, set by the considered doping. The occupations can be chosen randomly or respecting some given symmetry condition.

($ii$) Using the HHK recursion technique [42, 43, 44, 45], we compute the LDOS of the electronic system defined by the Hamiltonian, Eq. (3). The spin-resolved LDOS at a given site $i$ can be written as

In practice, the regularization parameter $\eta$ is considered as finite and its magnitude can be attributed the self-energy correction due to disorder, that is ubiquitous in graphene systems. The HHK method enables us to compute the LDOS with $O(\mathcal{N})$ operations and the DOS with $O(\mathcal{N}^{2})$. Hence, it is much more efficient than direct diagonalization, which demands $O(\mathcal{N}^{3})$ operations.

($iii$) Next, we determine $\langle n_{i\sigma}\rangle$, the average electron occupation number with spin $\sigma$ at the site $i$. In general, at a given temperature $T$, $\langle n_{i\sigma}\rangle$ reads

with $f_{\mu}(\epsilon)=\left\{\exp{\left[\beta(\epsilon-\mu)\right]}+1\right\}^{-1}$, where $\beta=1/k_{\rm B}T$, $k_{\rm B}$ is Boltzmann constant, and $\mu$ is chemical potential. At zero absolute temperature, $\mu$ is equal to the Fermi energy $\varepsilon_{F}$ and the Fermi distribution becomes a step function. For simplicity, here we consider $T=0$. As standard, the Fermi energy $\varepsilon_{F}$ (or the chemical potential $\mu$) is determined by the doping.

($iv$) We examine whether the output occupation numbers $\langle n_{i\sigma}\rangle$ coincide with the input ones, within a tolerance of $10^{-6}$. If the answer is positive, we end the self-consistent loop. If not, we return to ($ii$). The computation of the updated spin densities is then repeated iteratively until all values of $\langle n_{i\sigma}\rangle$ are converged.

The self-consistent solution provides the spin polarization at the $i$th site,

The $p_{z,i}$ is translated in a local magnetization $m_{z,i}=-g\mu_{\rm B}p_{z,i}$, where $\mu_{\rm B}$ is the Bohr magneton and the electron $g$-factor is $g=2$. Hence, the total magnetization per moiré cell reads

We study two cases: $(i)$ $N_{e}=0$, corresponding to the CNP, and $(ii)$ $N_{e}=3$ corresponding to $3/4$-filling of the conduction miniband [24]. For the CNP case, a previous theoretical study working in reciprocal space has predicted that low-angle TBG systems have an antiferromagnetic ground state [46]. In this case, we find computationally convenient to start the self-consistent loop with the configuration: $\langle n^{A}_{i\uparrow}\rangle=1/2+\Delta n^{A}$, $\langle n^{A}_{i\downarrow}\rangle=1/2-\Delta n^{A}$, $\langle n^{B}_{i\uparrow}\rangle=1/2+\Delta n^{B}$ and $\langle n^{B}_{i\downarrow}\rangle=1/2-\Delta n^{B}$, where $\Delta n^{B}=-\Delta n^{A}$. For the $3/4$-filling case, we start with a random set of $\langle n_{i\sigma}\rangle$ with the constraint $\sum_{i\sigma}\langle n_{i\sigma}\rangle=N+N_{e}$, where $N_{e}=3$.

We compute the local electronic properties of TBGs using the HHK recursion method [42, 43, 44, 45]. Being a real space approach, the HHK calculations take advantage of the symmetry properties of the odd commensurate TBG structures discussed above. We have numerically verified that the predicted equivalent sites indeed display the same LDOS. This is used to reduce the numerical effort by a factor of 6. For calculational convenience, the regularization factor is set to be $\eta\approx 5\dots 25$ meV. We return to this point later on.

To highlight the sites with the most prominent enhancement of the LDOS in TBGs, we consider the moiré region where the AA stacking is placed at the center of the TBG primitive unit cell. Figure 2(a) shows the LDOS for a moiré unit cell with $(m,n)=(22,23)$ corresponding to a twist angle $\theta\approx 1.47^{\circ}$ containing 6076 carbon atoms. The AB- (or BA-) stacking region is zoomed in the inset. Close to the CNP, well-localized states are found in the AA-stacking region, leading to an enhanced LDOS on atoms around the AA stacking, with much smaller LDOS at AB- and BA-stacking regions, in agreement with previous studies [18, 19].

Figures 2(b) and 2(c) show the LDOS at the CNP, $\nu_{i}({\rm CNP})$, for rigid and relaxed lattice structures with $\theta\approx 1.47^{\circ}$ within the circular area of radius of $31.3$ Å, centered at the AA dimer site (or ${\rm A}^{b}{\rm A}^{t}$ site). The rigid lattice structure shows a LDOS maximum, $\nu_{\rm AA}({\rm CNP})$, that is twice as large as the one computed for the relaxed lattice. For both rigid and relaxed lattices, $\nu_{i}({\rm CNP})$ decays radially with respect to the AA dimer site. However, the LDOS decays more rapidly for the reconstructed lattices due to in-plane relaxation, which reduces the AA-stacking area. In contrast, there is no significant contrast in the LDOS at the CNP in the AB-stacking region, $\nu_{\rm AB}({\rm CNP})$, for either rigid or relaxed lattices.

Figures 3(a) and (b) show the LDOS at the AA dimer site as a function of the energy $\varepsilon$ around the CNP energy, $\nu_{\rm AA}(\varepsilon)$, for $\theta\approx 1.30^{\circ}$, $1.05^{\circ}$, $0.88^{\circ}$ and $0.60^{\circ}$, that belong to the family of odd commensurate moiré structures with $n-m=1$, for both rigid and relaxed lattices, respectively. Here we use $\eta=5$ meV. The appearance of the central peak can be interpreted in terms of the continuum model [21], which associates the enhancement of the LDOS at the vicinity of the magic angle to the formation of a flat band at the CNP. Figure 3(c) displays the evolution of the AA dimer site LDOS at the CNP as a function of the small twist angles $\theta$ for the $n-m=1$ family of rigid and relaxed lattice structures.

Figures 3(a) to (c) show that when the twist angle decreases, the LDOS of the AA-stacking region increases significantly at the CNP. For the rigid lattice structures, the maximum of the peak appears at the twist angle of $\theta\approx 1.16^{\circ}$ [red vertical line in Fig. 3(c)] corresponding to $(m,n)=(28,29)$ with 9748 carbon atoms within the moiré cell. For the $n-m=1$ structures considered here, this twist angle is the closest to the first magic angle $\theta^{\rm MA}_{1}\approx 1.1^{\circ}$ predicted by the continuum model [21, 23]. We find that the flat-band peak width is larger in the relaxed TBG systems than in the rigid counterpart. The maximum of the LDOS peak appears at $\theta\approx 1.08^{\circ}$ [black vertical line in Fig. 3(c)], corresponding to $(m,n)=(30,31)$ with 11164 carbon atoms within the moiré cell, indicating that lattice relaxation shifts the LDOS peak towards lower twist angles.

Figure 3(c) shows another peak maximum at the twist angle close to $0.53^{\circ}$ (red vertical line), corresponding to $(m,n)=(62,63)$ with $46876$ carbon atoms within the moiré cell for the rigid TBG systems. This twist angle is the closest to the second magic angle $\theta^{\rm MA}_{2}\approx 0.55^{\circ}$ predicted by the continuum model [21]. In Fig. 3(a), as the twist angle decreases, the satellite peaks around the main one, both in the valence and conduction band, become increasingly intense and move closer to the CNP. In distinction, for relaxed TBG systems, neither the satellite peaks around the main one nor the second magic angle are observed in the electronic properties of twist angles below $\theta^{\rm MA}_{1}\approx 1.1^{\circ}$ [71, 72]. In contrast, we find that the van Hove singularities at the AB-stacking region are not depleted, but merely shifted due to lattice relaxation (not shown here), in line with the DOS reported in Ref. [55].

Let us now depart from the CNP and discuss finite doping. Figure 4(a) shows the DOS of the single-particle miniband of a TBG with $\theta\approx 1.16^{\circ}$. Anticipating the appearance of a gap at the CNP when the interaction is switched on, we define a valence and a conduction miniband separated at the CNP. The red area represents the half-band filling of a TBG system where the Fermi energy is at the CNP. This filling corresponds to 4 electrons per moiré cell in the valence miniband. The green area represents the electron doping at the $3/4$ filling of the conduction miniband. This doping corresponds to 3 electrons per moiré cell in excess of the CNP at a Fermi energy $\varepsilon_{F}=872$ meV.

Figure 4(b) shows the carrier distribution $\Delta n_{i}$ in excess to the CNP distribution within the moiré unit cell for $\theta\approx 1.16^{\circ}$. The areas of the green disks are proportional to the carrier occupation number at the atomic site $i$, $\Delta n_{i}$. Owing to the localized LDOS, see Fig. 2, the carrier distribution is also concentrated at the AA-stacking region.

Let us now examine the ferromagnetic phase in rigid and relaxed TBG systems with a twist angle $\theta\approx 1.16^{\circ}$ at $3/4$-filling of the conduction miniband whose origin has been the subject of discussion in the literature.

We consider an odd commensurate lattice with $(m,n)=(28,29)$ corresponding to 9748 carbon atoms. Our results are obtained by following the self-consistent procedure described in Sec. II.3 for a finite doping. Here, we set $U/V^{0}_{pp\pi}=1$. It is worth to note that we find that a small regularization parameter, such as $\eta=5$ meV, is not essential for obtaining converged results in the interacting case, thereby saving computational cost. In our study, we have carefully selected an $\eta$ value that is sufficiently large to ensure accurate integrated DOS and occupation near the Fermi energy, while remaining small enough not to influence the mean-field calculations. The optimal $\eta$ value may vary between 5 and 25 meV depending on the stacking region considered.

Figures 5(a) and (c) display the local spin polarization, $p_{z,i}$, within a circular area with a radius of $39.9$ Å centered at the AA dimer site for rigid and relaxed TBG systems, respectively. These results show that the ground state is ferromagnetic as characterized by the imbalance between $\langle n_{i\uparrow}\rangle$ and $\langle n_{i\downarrow}\rangle$ in the AA stacking region. Since the spin polarization is largest at the regions of enhanced LDOS, we attribute the emergence of ferromagnetism to the Stoner instability mechanism.

Figures 5(b) and (d) show the local spin polarization, $p_{z,i}$, as a function of the distance $D$ from the AA dimer site for rigid and relaxed TBG systems, respectively. The spin polarization is largest at the AA dimer site and becomes smaller as $D$ increases and eventually vanishes as the $i$th site approaches the AB non-dimer position.

Since the ground state is ferromagnetic, the total magnetization is not zero. For the rigid lattice structure, we have a maximum local spin polarization of approximately $1.80\times 10^{-4}$ at the AA dimer site, and we find that the total spin polarization per moiré cell is $\sim 0.19$ and $M_{\rm PUC}=-0.22$ ${\rm eV}\cdot{\rm T}^{-1}$. However, for the relaxed lattice structure, the maximum local spin polarization is $\sim 1.19\times 10^{-4}$ at the AA dimer site with total spin polarization per moiré cell is $\sim 0.13$ and $M_{\rm PUC}=-0.15$ ${\rm eV}\cdot{\rm T}^{-1}$. Furthermore, the local spin polarization decays rapidly within the circular area, see Figs. 5(c) and (d), due to the reduced area of the AA-stacking region.

Figure 6(a) and (c) present the LDOS calculated at the AA dimer site, ${\rm LDOS}({\rm AA},\epsilon)$, within an energy window that contains the flat minibands near the Fermi energy (set to $\varepsilon_{F}=0$ for clarity) for rigid and relaxed structures, respectively. We argue that the LDOS(AA,$\varepsilon$) reproduces the behavior of DOS($\varepsilon$) in the vicinity of $\varepsilon_{F}$, similarly as in the CNP case studied in Ref. [46] (see, App. A). The separation of the two LDOS peaks around the Fermi energy is approximately $\Delta\approx 29$ meV for a rigid lattice structure and $\Delta\approx 36$ meV for a relaxed one. We speculate that this modest enhancement is due to the fact that, despite the slightly smaller DOS observed in the relaxed structure, the reduced size of AA-stacking region amplifies Coulomb repulsion effects. We find that the gap sizes are spatially constant within the AA-stacking region. The existence of two LDOS peaks in the valence miniband suggests the existence of two minibands that overlap in the CNP, as indicated in Ref. [3].

Figure 6(b) depicts the spin up and spin down LDOS at the AA dimer site, denoted as $\nu_{{\rm AA},\uparrow}(\epsilon)$ and $\nu_{{\rm AA},\downarrow}(\epsilon)$, respectively, for rigid lattice structure. A clear spin-imbalance is evident between the spin-up and spin-down occupation numbers in the filled miniband (fmb). Quantitatively, the occupation numbers are $\langle n^{\rm fmb}_{{\rm AA},\uparrow}\rangle=9.7\times 10^{-4}$ and $\langle n^{\rm fmb}_{{\rm AA},\downarrow}\rangle=6.1\times 10^{-4}$. This imbalance translates to a spin polarization $p_{z,{\rm AA}}=1.8\times 10^{-4}$ at the AA dimer site, corresponding to a local magnetic moment $m_{z,{\rm AA}}=-0.21$ ${\rm meV}\cdot{\rm T}^{-1}$. In Fig. 6(d), which corresponds to the relaxed lattice structure, the occupation numbers are $\langle n^{\rm fmb}_{{\rm AA},\uparrow}\rangle=6.4\times 10^{-4}$ and $\langle n^{\rm fmb}_{{\rm AA},\downarrow}\rangle=4.0\times 10^{-4}$, with a spin polarization $p_{z,{\rm AA}}=1.2\times 10^{-4}$ at the AA dimer site.

In summary, the Stoner mechanism is robust against the suppression of the flat-band LDOS peak due to lattice relaxation. The ferromagnetic phase is preserved, but the magnitude of the local magnetic moments becomes smaller than the ones computed for rigid structures.