1 Introduction↩︎

Two-dimensional (2D) materials have become a central topic in solid-state physics due to their distinct properties [1]. Compared with three-dimensional systems, they offer advantages such as tunability by electric fields and proximity coupling to other van der Waals materials, which can reshape their band structure and interactions [2], [3]. Few-layer graphene provides perhaps the simplest examples, including twisted bilayers [4] and rhombohedral, ABC stacked multilayers [5]–[9], where emergent phases such as superconductivity and unconventional magnetism appear. Lattice defects strongly influence the electronic properties of 2D materials [10]–[16]. In multilayer systems, an additional category of defects emerges in the form of stacking faults, which arise from deviations in the regular layer-by-layer arrangement. The significance of stacking order is particularly pronounced in van der Waals layered materials, where it dramatically alters electronic properties [17], [18]. Such changes in the interlayer stacking can be harnessed to tune the properties of heterostructures, a prime example being twisted structures [3], [4]. Multilayer graphene serves as one of the simplest examples of this phenomenon. When a graphene layer is placed on top of graphite surface with at least two layers, it can adopt two thermodynamically stable stacking configurations: the hexagonal (ABA) and the rhombohedral (ABC) positions [19], [20]. This pattern repeats with each additional graphene layer, leading to an exponentially growing number of possible stackings as the layer count increases [21]. Graphite in the bulk typically exhibits ABA (Bernal) stacking, while a small portion adopts ABC (rhombohedral) stacking [22]. Following this pattern, four-layer graphene should display either ABAB or ABCA sequences. However, a third type: ABCB has been experimentally observed, showing ferroelectric properties not present in the defect free structures [21].

Graphene multilayers with ABC stacking sequence, have recently gained prominence as a platform for probing emergent, strongly-correlated electronic phenomena [6]–[9]. Both the lattice defects and electronic properties of this material have started to be investigated [19]–[21], [23]–[25]. Various theoretical approaches have been employed, ranging from simple continuum [26], [27] and tight-binding (TB) models with parameters fitted to either experiments [28]–[30] or first-principle calculations [31], [32]. The most advanced techniques, including the GW approximation [33], [34], density functional theory (DFT) using hybrid functionals [35], [36], and random phase approximation (RPA) [34], have also been applied to few-layer systems. As theoretical methods become more accurate and precise, computational costs increase significantly. Calculations employing hybrid functionals demand substantial high-performance computing (HPC) resources and considerable time even on advanced systems. Conversely, TB calculations offer analytically tractable or computationally inexpensive approaches to investigate electronic structure, making them capable of modelling various lattice defects in large-scale systems containing tens of thousands of atoms [37]–[39]. While TB calculations offer computational efficiency, they traditionally neglect electron-electron interactions, which can be crucial for understanding many-body phenomena in graphene systems. To address this limitation, we apply the methods suggested by Hubbard to include electron-electron interactions in the TB Hamiltonian, enabling the description of conductor-insulator transitions [40]. This model incorporates repulsive contact interactions parametrized by the so-called Hubbard-\(U\) term, with the resulting Hamiltonian often called the TB+U Hamiltonian. This \(U\) is a free parameter in the model and its value should be fitted to experimental results or ab initio calculations.

In this work, we propose a parametrization of the TB Hamiltonian of multilayer graphene systems, using the distance-dependent Slater-Koster parametrization and extend it by adding the Hubbard term. We show that a universal Hubbard-\(U\) value can be used to describe the magnetic properties of multilayer graphene systems with arbitrary stacking configurations. Although ABC graphite is predicted to host several competing many-body ground states [5], [41], [42], we focus on the zero-doping case without external electric or magnetic fields. In this regime, experiments consistently report an insulating ground state at charge neutrality in pristine ABC-stacked samples [43]–[46]. This insulating phase, known as the layer antiferromagnet (LAF) [43], is characterized by antiparallel spin alignment between the outermost layers of the ABC ladder. The LAF state matches ab initio predictions of a magnetic insulator [35] and has been observed in tetralayer ABC [44], [45] and pentalayer samples [46] among others, with possible beyond-mean-field corrections in thicker systems [47], [48]. Here we investigate how the magnetic moments of the LAF phase are modified as the stacking changes from ABC to ABA. More generally, our framework can also describe complex stacking geometries, including twisted multilayers [49] and lateral domain walls between ABC and ABA regions, which have been detected in scanning probe experiments [50], [51].

2 Model and methods↩︎

To investigate the magnetic properties of multilayer graphene systems, we first focus on the two most common stacking arrangements: ABA (Bernal) and ABC (rhombohedral) stacking. These structures were modeled using the following lattice vectors: \[\mathbf{a}_1=\frac{a_{\mathrm{CC}}}{2} \begin{pmatrix} 3\\ \sqrt{3} \end{pmatrix},\;\;\;\;\;\;\;\mathbf{a}_2=\frac{a_{\mathrm{CC}}}{2} \begin{pmatrix} 3 \\ -\sqrt{3} \end{pmatrix} \label{eq:lat95vecs}\tag{1}\] Where \(a_{\mathrm{CC}}=a_0/\sqrt{3}\) is the carbon-carbon distance used in our work, and \(a_0=2.461\) Åis the in-plane lattice constant. The interlayer distance was set to \(d_0=3.347\) Å.

These lattice parameters were chosen to be consistent with the values used in previous HSE calculations [35], [36]. The two stacking arrangements differ in their interlayer registry. In ABA graphene, the second layer is shifted by the vector \(a_{\mathrm{CC}}(1,0)\) relative to the first layer, while the third layer is positioned directly above the first, and the fourth above the second. In ABC stacking, the first and second layers maintain the same relative positioning as in ABA, but the third layer is shifted by \(2a_{\mathrm{CC}}(1,0)\) relative to the first layer. The unit cell structures of both systems are illustrated in the right column of Fig. 1, where even (odd) numbered sites indicate the A (B) sublattices within each layer.

To investigate the robustness of the magnetic structure in ABC graphene, we studied two types of systems consisting of both ABC and ABA parts, as simple examples for mixed stackings. The first type consists of an \(M\)-layer ABC region stacked on top of an \(N\)-layer ABA region. The second type features an \(N\)-layer ABA region sandwiched between two ABC regions of thickness \(M_1\)and \(M_2\) layers, respectively. Since these mixed stackings are composed of alternating rhombohedral (ABC) and Bernal (ABA) multilayer segments, we adopt the naming convention R\(M_1\)-B\(N\)-R\(M_2\), where R denotes rhombohedral stacking, B denotes Bernal stacking, and the subscripts indicate the number of layers in each region. For structures with only two regions, the convention simplifies to R\(M\)-B\(N\). Using this notation, the two example structures shown in Fig. 2 are designated as R7-B10 and R5-B6-R5. Our approach combines tight-binding (TB) calculations with Hubbard interaction terms to model the magnetic properties of multilayer graphene systems. To parametrize and validate this TB+\(U\) model, we utilize reference data from recent hybrid density functional theory (DFT) calculations. Specifically, we employ the magnetic moments calculated using PBE0 hybrid functionals by Pamuk et al. [35] for ABC multilayers and Campatella et al. [36] for ABA multilayers. Both studies systematically investigated 3-8 layer systems and examined the effects of different basis sets and functionals on the electronic and magnetic structure. The TB+\(U\) model employs a tight-binding Hamiltonian extended with on-site Hubbard interaction terms to account for electron-electron repulsion. The Hubbard-\(U\) parameter was determined by fitting the calculated magnetic moments to the hybrid DFT results from the aforementioned studies. The detailed formulation of the TB+\(U\) model, the fitting procedure, and the computational implementation are described in 5.1, 5.2 and 5.3.

3 Results↩︎

We now examine the magnetic properties of different multilayer graphene systems using our TB+\(U\) model with the fitted Hubbard parameter \(U=5.84\text{ eV}\). We begin with pure stacking configurations before investigating the behaviour of mixed stacking systems. The magnetic moments of ABC multilayers calculated from our TB+\(U\) model show agreement with the hybrid DFT results, as demonstrated in Fig. 1. The magnetic structure exhibits a consistent pattern across all layer numbers: within each layer, the system displays antiferromagnetic ordering, with magnetic moment magnitudes that decay toward the center of the multilayer. Additionally, the overall magnetic moments increase systematically with the number of layers. For systems with more than 12 layers the change of the maximum value of the magnetic moments between consecutive layer number is smaller than the convergence criteria.

The ABA multilayers exhibit fundamentally different magnetic behaviour compared to their ABC counterparts. While they also show antiferromagnetic ordering within each layer (Fig. 1), the spatial distribution of magnetic moments is inverted: the magnitude increases as we move toward the center of the system. Notably, although the hybrid functional calculations predict no magnetic moments for three, four, and seven layers, our TB+\(U\) calculations reveal non-zero magnetic moments for these systems as well. Importantly, a single Hubbard parameter adequately describes the magnetic properties across both stacking configurations and multiple layer thicknesses (3-8 layers), suggesting reasonable transferability of our TB+\(U\) approach.

Having validated our model against pure systems, we now turn to ABC-ABA mixed stackings, which represent a computationally prohibitive challenge for hybrid DFT methods due to their large system and configurational sizes. To investigate how the distinct magnetic behaviours of ABC and ABA regions interact, we calculated the magnetic structure of various ABC-ABA mixed stackings shown in Fig. 2. The results reveal that the magnetic moment patterns are largely preserved within each stacking region, with the ABC portions maintaining their characteristic decay toward the center while the ABA regions exhibit their typical center-enhanced behaviour. In the ABA regions located far from the ABC-ABA junction (typically 2-3 layers away), the magnetic structure remains virtually indistinguishable from that of a free-standing ABA system. However, as the number of ABC layers in the mixed stackings increases, the perturbation to the ABA magnetic structure becomes more pronounced, indicating a non-local influence of the ABC stacking on the overall magnetic properties.

3.1 Effects of mechanical distortion↩︎

Since the extent of the surface flat band, that gives rise to the magnetic instability is proportional to the interlayer hopping terms [26], it is expected that the magnetic moments in the LAF phase should scale with the interlayer hopping term. Having established the magnetic properties of unstrained systems, we now investigate how mechanical deformation affects the magnetic structure of multilayer graphene.

3.1.1 Uniaxial strain in the armchair direction↩︎

We applied uniaxial strain to study its effects on the magnetic structure of pure ABC and ABA multilayers as well as selected mixed stackings. The strain was applied in the armchair (\(x\)) direction, parallel to one of the carbon-carbon bonds, with magnitudes of \(\pm 1\%\). The calculations were performed on multilayer systems with N=4, 5, 6, 9, 11 layers for both pure ABC and ABA configurations, as well as the two selected mixed stackings. The results, shown in Fig. 3, display the maximum magnetic moment values for each system to provide a clear comparison of the strain effects across different configurations. Under tensile strain (stretching), the magnitude of magnetic moments increases in ABC, ABA, and mixed stacking systems. Conversely, compressive strain exhibits asymmetric effects on the two stacking types: in ABA systems, magnetism can be suppressed by sufficiently large compression, while ABC systems maintain their magnetic character throughout the experimentally accessible strain range [52].

3.1.2 Deformation of the van der Waals spacing↩︎

In addition to in-plane strain effects, we investigated how interlayer separation influences the magnetic properties by applying compression along the \(z\) direction (perpendicular to the graphene planes). This allows us to probe the role of interlayer coupling in stabilizing the magnetic states. We applied \(\pm 5\%\) changes to the interlayer distance, where \(+ 5\%\) corresponds to an increased layer separation of \(1.05~d_0\) relative to the equilibrium distance \(d_0 = 3.347\) Å. We chose higher strain values because compression values of up to 5%are available experimentally in van der Waals systems when hydrostatic pressure is applied [53]. The same multilayer systems (N=4, 5, 6, 9, 11 layers) were examined for both ABC and ABA stackings. The results presented in Fig. 4 show the maximum magnetic moment for each system configuration, allowing for straightforward comparison of how interlayer distance affects magnetic strength across different layer numbers and stacking types.

The calculations reveal a universal trend across all layer numbers and stacking configurations: increasing the interlayer distance consistently reduces the magnitude of magnetic moments, while decreasing the separation enhances them. This behaviour can be understood through the weakening of interlayer orbital interactions as the layers are moved further apart. The reduced overlap between \(\pi\) orbitals in adjacent layers diminishes the electronic coupling that drives the magnetic instability, leading to smaller magnetic moments. Conversely, when layers are brought closer together, the enhanced interlayer interactions strengthen the conditions for magnetic ordering. This \(z\)-direction response demonstrates that the magnetic properties of multilayer graphene are sensitive not only to in-plane deformation but also to the three-dimensional structural parameters.

4 Conclusion and outlook↩︎

We have shown that a universal Hubbard-\(U\) value can be used to capture the magnetic properties predicted by hybrid DFT calculations in defect free ABA and ABC stacked graphene multilayers. Our calculations showed that the magnetic structure of ABC graphene is robust when it is combined with ABA graphene into structures with mixed stacking. At the ABC–ABA interface, the magnetic moments within two layers of the boundary are reduced compared to those in free-standing ABA graphene, but farther away they converge to the same values. Interestingly the edge momenta of ABC regions are less perturbed and remain equal to their free-standing values when adjacent to ABA stacking. This suggests that "buried" flat-band [54] magnetic states could be prepared, where neighbouring hexagonal regions can protect the magnetism from charge disorder.

To incorporate the effect of different mechanical distortions we applied uniaxial strain in the armchair direction and pressure perpendicular to the graphene planes. By decreasing the layer distance, the orbitals within adjacent layers get closer to each other, therefore the interaction gets larger between the electrons. This causes the magnetic moments on each site to increase. On the other hand calculations with uniaxial strain showed opposite behaviour in the magnetic structure. By stretching the system the magnetic moments become larger, while by compressing it they become smaller. Based on the results from the calculations incorporating both strain and pressure, the absence of magnetism in 3, 4 and 7 layer ABA graphene in Ref. [36], could be explained with geometrical reasons, for details see Appendix 7.

The emergence of magnetic moments in the hexagonal phase, which our calculations also support is a noteworthy and often overlooked feature. This should motivate future experimental work and many-body calculations, since charge transport measurements already suggest signatures of correlation effects [55].

5 The tight-binding model↩︎

5.1 The model↩︎

To determine the electronic properties of the systems a distance dependent tight-binding model was used. The tight-binding Hamiltonian for electrons in a single sheet of graphene considering that electrons can hop between first, second and third nearest neighbours has the following form:

The operators \(a^{\dagger}_{i,\sigma}\;(a_{i,\sigma})\) create (annihilate) an electron on sublattice A, site \(i\) with spin \(\sigma\). The notation is the same for the \(b\) operators. \(\langle\cdots\rangle\) denotes summing over the first-, \(\langle\langle\cdots\rangle\rangle\) denotes summing over the second- and \(\langle\langle\langle\cdots\rangle\rangle\rangle\) denotes summing over the third-nearest neighbours. The hoppings have the following form [56]: \[\gamma_i(r) = \gamma_{i,0}\cdot \mathrm{e}^{\beta_i\left(\frac{r}{a_{\mathrm{CC}}}-1\right)} \label{eq:inplane95hopping}\tag{3}\] where the subscript \(i=1,2,3\) denotes the first, second and third neighbours, \(\gamma_{i,0}\) is the value of the hopping integral in the equilibrium position, \(\beta_i\) is the strength of the decay of the hopping integral and \(r\) is the distance between the two atoms. Both \(\gamma_{i,0}\) and \(\beta_i\) were obtained from fitting the tight-binding band structure to DFT results using the least-squares method (see Sec. 5.2). To include the hoppings between the different layers, the Hamiltonian 2 was extended with the following terms: \[\hat{\mathrm{H}}_{\mathrm{TB},\mathrm{interlayer}}=\sum_{i,j,\sigma}\gamma_{ij}c^{\dagger}_{i,\sigma}c_{j,\sigma}+\mathrm{H.c.}\] The operators \(c^{\dagger}_{i,\sigma},\;(c_{i,\sigma})\) create (annihilate) an electron on site \(i\) with spin \(\sigma\). Sites \(i\) and \(j\) are in adjacent layers. The interlayer hopping between atoms at positions \(\mathbf{r}_i\) and \(\mathbf{r}_j\) has the following form following the work of Slater and Koster [57]:

Here the angle between the vector \(\mathbf{r}_{ij}=\mathbf{r}_i-\mathbf{r}_j\) and the \(z\) axis is \(\theta_i\) as it can be seen in Fig. 5. The magnitude of vector \(\mathbf{r}_{ij}\) is \(r_{ij}\). To eliminate the very small hoppings a cut-off function \(\Theta(\mathbf{r}_{ij})\) was introduced that is 1 for \(r_{ij}<7\) Å  and 0 otherwise. By choosing this cut-off distance we eliminate the very small matrix elements. For a few given \(\theta\) values the distance dependence of the \(\gamma_{ij}\) hoppings can be seen in Fig. 6.

\(\tilde{V}_{pp\pi}\), \(\tilde{V}_{pp\sigma}\) and \(\alpha\) are parameters that had to be fitted to DFT calculations.

5.2 Fitting procedure↩︎

Fig. 7 shows the band structure of 6 layer ABC and ABA graphene along the full \(\Gamma-K-M-\Gamma\) high symmetry path in the Brillouin zone. The only part of the high symmetry that has energy bands in the relevant energy range (\(\pm\) 1 eV) is the small vicinity of the \(K\) point. Therefore the TB parameters were obtained by fitting the band structure of the TB Hamiltonian to the dispersion relation calculated using the SIESTA code [58]–[60] using the least squares method in a small vicinity of the \(K\) point.

First the nearest neighbour in-plane hopping were fitted to the monolayer graphene band structure. To include the distance dependence, biaxial strain was applied to the system, and the value of \(\gamma_{1,0}\) was obtained by fitting the TB band structure to the DFT dispersion relation. After that the a function of the form of Eq. 3 was fitted to the different hopping values to get the value of \(\beta_1\). To further improve the precision of the model, the second the third nearest neighbour hoppings were introduced. The \(\gamma_{2,0}\), \(\gamma_{3,0}\) hopping values and their respective \(\beta\) values were fitted to different DFT band structures when uniaxial strain was applied. Regarding the parameters describing the interlayer hoppings the DFT band structure of AB and AA bilayer graphene was used. The three parameters (\(\tilde{V}_{pp\pi}\), \(\tilde{V}_{pp\sigma}\), \(\alpha\)) were simultaneously fitted to the AB and AA band structure by using four low energy bands from the AB and four low energy bands from the AA band structure. Fig. 8. shows the fitting procedure the interlayer hopping parameters.

5.3 TB+U↩︎

To include the electron-electron interaction we extended the Hamiltonian with the Hubbard term [40]. Therefore the Hamiltonian has the following form: \[\hat{\mathrm{H}} = \hat{\mathrm{H}}_{\mathrm{TB}}+\sum_{i,\uparrow,\downarrow}U_in_{i\uparrow}n_{i\downarrow}\] Where \(\hat{\mathrm{H}}_{\mathrm{TB}}=\hat{\mathrm{H}}_{\mathrm{TB},\mathrm{interlayer}}+\hat{\mathrm{H}}_{\mathrm{TB},\mathrm{in-plane}}\). By applying the meanfield approximation, the Hamiltonian gets the following form: \[\hat{\mathrm{H}}_{\mathrm{MF}} = \hat{\mathrm{H}}_{\mathrm{TB}}+U\sum_{i,\sigma}\langle n_{i\sigma}\rangle n_{i\bar{\sigma}}\] Where instead of having different \(U_i\) values on each site there is one global \(U\) for all of them and \(\langle n_{\sigma}\rangle\) is the \(\sigma=\uparrow,\downarrow\) spin density on site \(i\). In the above Hamiltonian the only free parameter is the Hubbard \(U\). Its value was found by simultaneously fitting the magnetic moments of the five layer ABC and ABA systems obtained from the TB+U model to the results of the hybrid DFT calculations [35], [36]. The above Hamiltonian contains the expectation value of the spin-resolved density operator \(n_{\sigma}\). This expectation value is computed from the eigenvectors of \(\hat{\mathrm{H}}_{\mathrm{MF}}\), therefore a self-consistent field (SCF) method should be applied. During the SCF calculations a \(300\times 300 \times 1\) Monkhorst-Pack [61] Brillouin zone sampling was used. The temperature was et to be 10\(^{-5}\;k_BT\). The TB and TB+U calculations were performed using the sisl[62] and the hubbard[63]–[65] codes.

6 Details of the DFT calculations↩︎

The DFT calculations were performed using the SIESTA code. for the exchange-correlation functional the Perdew-Burke-Ernzerhof (PBE) [66] parametrization of the generalized gradient approximation (GGA) was used with double-\(\zeta\) basis set. The norm-conserving pseudopotential were collected from the PseudoDojo [67] project.
To ensure well converged results, a \(k\)-grid of \(300\times 300 \times 1\) was used using the Monkhorst-Pack sampling. A real space grid cut-off energy of 500 Ry was used. The structures were relaxed until the maximum force was less than 0.008 eV/Å. During the calculations vacuum separating the samples was chosen to be at least 30 Å.

7 Simultaneous uniaxial strain and compression↩︎

In [35] and [36] the effect of different geometries was not studied. In accordance with this we fixed the lattice parameters. With these settings the TB+U model could not reproduce the absence of magnetism in 3, 4 and 7 layer ABA graphene. Calculation incorporating both uniaxial strain and interlayer compression (Fig. 9) show that for a large set of possible lattice parameters the magnetic moments are negligibly small or zero for the ABA 3 layer graphene. The same effect occurs in 5 layer ABA but for a much smaller set of lattice parameters. Meanwhile the magnetism is robust for both ABC multilayers. Therefore an explanation for the absence of magnetism in the 3,4 and 7 layer ABA systems in [36] is the selected values of the lattice parameters.

References↩︎