# Magnetic ordering tendencies in hexagonal boron nitride-bilayer graphene
moiré structures

Maria Spethmann  Institute for Theoretical Solid State Physics, RWTH Aachen
University, and JARA Fundamentals of Future Information Technology, 52056
Aachen, Germany  Carsten Honerkamp  Institute for Theoretical Solid State
Physics, RWTH Aachen University, and JARA Fundamentals of Future Information
Technology, 52056 Aachen, Germany  Dante M. Kennes  Institute for Theory of
Statistical Physics, RWTH Aachen University, and JARA Fundamentals of Future
Information Technology, 52056 Aachen, Germany  Max Planck Institute for the
Structure and Dynamics of Matter, Center for Free Electron Laser Science,
Luruper Chaussee 149, 22761 Hamburg, Germany

(May 29, 2021)

###### Abstract

When hexagonal boron nitride (hBN) and graphene are aligned at zero or small
twist angle, a moiré structure is formed due to the small lattice constant
mismatch between the two structures. In this work, we analyze magnetic
ordering tendencies, driven by onsite Coulomb interactions, of encapsulated
bilayer graphene (BG) forming a moiré structure with one (hBN-BG) or both hBN
layers (hBN-BG-hBN), using the random phase approximation. The calculations
are performed in a fully atomistic Hubbard model that takes into account all
$\pi$-electrons of the carbon atoms in one moiré unit cell. We analyze the
charge neutral case and find that the dominant magnetic ordering instability
is uniformly antiferromagnetic. Furthermore, at low temperatures, the critical
Hubbard interaction $U_{c}$ required to induce magnetic order is slightly
larger in those systems where the moiré structure has caused a band gap
opening in the non-interacting picture, although the difference is less than
6%. Mean-field calculations are employed to estimate how such an interaction-
induced magnetic order may change the observable single-particle gap sizes.

##  I Introduction

hBN is frequently used as a substrate for graphene devices, because it offers
a flat and clean surface and enables a high charge carrier mobility [1]. When
single layer graphene (SG) and hBN are stacked at no or small twist angle to
each other, they form a moiré structure due to a small lattice constant
mismatch. hBN has a lattice constant of around $a_{hBN}=0.2504\,\mathrm{nm}$
[2] whereas SG has a lattice constant of $a=0.246\,\mathrm{nm}$ [3]. The
formed moiré structure, with lattice constants up to $14\,\mathrm{nm}$ [4],
changes the electronic properties of the material and has been analyzed in
several theoretical [5, 6, 7, 8] and experimental [9, 10, 11] papers. For
example, a band gap opens at the central Dirac points of SG [11]. The size of
the central band gap was analyzed in calculations, which stressed the
importance of structural relaxation of the carbon atoms within the moiré
structure, but also included many-body perturbation theory and electron-
electron interactions [12, 6]. Further, secondary, gapped Dirac points have
been reported to form at the mini Brillouin zone edges [1, 8]. In addition,
the wavelength of the moiré structure and the underlying potential modulation
is well suited to observe the Hofstadter butterfly, a self-similar pattern of
the energy levels as a function of a magnetic field [10, 11].

In contrast to hBN-SG systems, literature on hBN-BG and hBN-BG-hBN moiré
structures is less extensive. In Ref. [8] the band structure of hBN-BG was
calculated in a tight binding model that explicitly took the boron and
nitrogen atoms into account. Then an effective continuum model was used to
calculate the spectrum featuring Hofstadter’s butterfly. In Ref. [13]
$n$-layer structures built of graphene and hBN were analyzed theoretically,
including interaction effects and a vertically applied electrical field, which
contributes to the formation of flat bands and changes the Chern number.
Further, correlated states such as quantum anomalous Hall insulating states
and fractional anomalous quantum hall effect states were described. In the
experimental work of Ref. [14] band gap widths were measured at the central
Dirac point and the secondary Dirac points both in hBN-SG and hBN-BG samples.
It was found that in hBN-BG the secondary Dirac points are not fully gapped,
and the particle hole asymmetry is less pronounced compared to hBN-SG. The
measurements were accompanied by theoretical calculations and it was concluded
that the observed behavior can be explained qualitatively in a non-interacting
picture. Another experimental study created doubly aligned hBN-BG-hBN with
close-to-zero rotation angle [15]. They measured a non-local resistance which
they attribute to a “valley current” near the band gap at charge neutrality.
In Ref. [16] doubly aligned hBN-BG-hBN was produced with a small twist in one
hBN-layer and resistivity measurements displayed secondary Dirac points from
both moiré structures.

In this work we want to understand hBN-BG and hBN-BG-hBN moiré structures more
comprehensively by analyzing interaction-driven magnetic. This magnetic order
is relevant as a proxy for interaction-induced gaps in BG and is a promising
candidate for the cases where such interaction-induced gaps may have been
observed. In fact, in a number of experiments on freely suspended fewlayers
[17, 18, 19, 20, 18] and multilayers [21, 22, 23] of graphene, all in Bernal
stacking, evidence for an interaction-induced gap opening was found. This was
interpreted theoretically as due to layered-antiferromagnetic order [24, 25,
26, 27, 28, 29, 30]. In theory, also SG hosts this ordering tendency [31], but
the density of states at low energies is not as high as in the layered
systems, and hence the ordering instability does not occur under normal
circumstances. Clearly, it is of interest how this ordering tendency of
layered graphene is altered in presence of an hBN environment. The influence
of the aligned hBN is twofold: first, it modifies the band structure, which is
addressed in this work. Second, the hBN environment leads to additional
screening of the electron-electron interactions. This effect is e.g. dealt
with by Roesner et al. [32] and can be incorporated in our reasoning by
choosing an appropriately renormalized interaction constant.

Our calculations for the interaction effects are restricted to charge
neutrality and zero external magnetic field. First we extend the hBN-SG model
of Ref. [33] to a model of hBN-BG and hBN-BG-hBN. In this process we also
distinguish between different stacking configurations. Next, we calculate the
parameters for which the interacting spin susceptibility diverges in the
random phase approximation (RPA). This provides information about the dominant
magnetic order of these systems. Finally, we compare our RPA results to mean
field calculations, where we are able to calculate the band gap widths of the
magnetic phases.

##  II Methods

We set up a Hubbard model that describes the $\pi$-electrons of the carbon
atoms, where each atom provides one orbital:

$\displaystyle H=$
$\displaystyle\sum\limits_{\begin{subarray}{c}\boldsymbol{R}_{1},\boldsymbol{R}_{2},\\\
\boldsymbol{r}_{1},\boldsymbol{r}_{2},\sigma\end{subarray}}-t(\boldsymbol{d})c_{\boldsymbol{R}_{1},\boldsymbol{r}_{1}\sigma}^{\dagger}c_{\boldsymbol{R}_{2},\boldsymbol{r}_{2}\sigma}+U\sum\limits_{\boldsymbol{R},\boldsymbol{r}}n_{\boldsymbol{R},\boldsymbol{r}\uparrow}n_{\boldsymbol{R},\boldsymbol{r}\downarrow}$
$\displaystyle+\sum_{\boldsymbol{r},\boldsymbol{R}}\epsilon^{\text{hBN}}_{\boldsymbol{r}}c_{\boldsymbol{R},\boldsymbol{r}\sigma}^{\dagger}c_{\boldsymbol{R},\boldsymbol{r}\sigma}$
(II.1)

The field operators $c_{\boldsymbol{R},\boldsymbol{r}\sigma}^{(\dagger)}$
annihilate (create) an electron with spin $\sigma$ at a site with lattice
vector $\boldsymbol{R}$ and basis vector $\boldsymbol{r}$. The atomic
positions and the tight binding part including the hopping parameter
$t(\boldsymbol{d})$, which is a function of the distance $\boldsymbol{d}$, is
constructed as explained in the Refs. [34, 35], a model originally used to
describe twisted bilayer graphene. We fix the interlayer distance to be
$d_{AB}=0.335\,\mathrm{nm}$ [3]. The interaction is a Hubbard on-site
interaction of strength $U$ and the operators
$n_{\boldsymbol{R},\boldsymbol{r}\sigma}$ are the spin-dependent electron
density operators. Although screening effects are less pronounced in graphene
materials, and therefore nearest neighbor interaction is comparable to onsite
interaction [36], a description using on-site interaction only is justified
when $U$ is properly renormalized [37, 38].

The boron and nitrogen electrons are not taken into account explicitly in our
model. Instead, the moiré structure is modeled by an effective potential term
$\epsilon^{\text{hBN}}_{\boldsymbol{r}}$ on the carbon atoms, as proposed for
hBN-SG in Ref. [33]. There, a commensurate moiré structure is assumed with
adjusted lattice constants to form a well-defined moiré unit cell. A close
approximation would be $\frac{a}{a_{\text{hBN}}}\approx\frac{55}{56}$ [2, 8].
But since this requires a very large unit cell (with a lattice constant a
factor $N_{\text{moir\'{e}}}=56$ larger than in graphene) and extensive
computational resources, we perform most calculations on a simplified model
with $\frac{a}{a_{\text{hBN}}}\approx\frac{19}{20}$
($N_{\text{moir\'{e}}}=20$). This means each unit cell in the moiré structure
contains $N_{\text{moir\'{e}}}\times N_{\text{moir\'{e}}}\times 4$ carbon
atoms. In Ref. [33] an effective potential on the carbon atoms from the hBN-
layer was proposed by performing a fit to DFT results. We repeated the last
part of the fit using the same fit data points but a different function that
keeps the threefold rotational symmetry:

$\displaystyle\Delta_{\boldsymbol{r}}=A\,\Big{(}$
$\displaystyle\sin(\boldsymbol{G}_{1}\cdot\boldsymbol{r}+\phi)+\sin(\boldsymbol{G}_{2}\cdot\boldsymbol{r}+\phi)$
$\displaystyle+\sin(-(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})\cdot\boldsymbol{r}+\phi)\Big{)}+C$
(II.2)

Here, $\boldsymbol{G}_{1}$, $\boldsymbol{G}_{2}$ are the reciprocal primitive
lattice vectors, and $A$, $\phi$ and $C$ are fit parameters, with best fit
results at $A=P_{A}=21.9\,\mathrm{meV}$, $\phi=P_{\phi}=1.655$, and
$C=P_{C}=-4.7\,\mathrm{meV}$. In the following, $\Delta_{\boldsymbol{r}}$ will
be called hBN-potential and is shown in Fig. 1. The effective potential
incorporates $\Delta_{\boldsymbol{r}}$ with alternating sites on the two
carbon atoms of the honeycomb lattice,

$\displaystyle\epsilon_{\boldsymbol{r}}^{\text{hBN}}=s_{\boldsymbol{r}}\Delta_{\text{r}},$
(II.3)

with $s_{\boldsymbol{r}}=+1/\text{-}1$ for $\boldsymbol{r}$ being at an
A-/B-site. To set up a model of hBN-BG from this model of hBN-SG, we assume
that the substrate only affects the first graphene layer [14], thus the hBN-
potential is only added on one graphene layer. hBN-BG has two different
stacking configurations, as pointed out in Refs. [8, 14]. These two structures
are shown in Fig. 2, which depicts one region of the moiré structure. For hBN-
BG-hBN moiré structures we need to differentiate between more stacking
configurations and we concentrate on those configurations where we can
determine a clear stacking axis, as shown in Fig. 3. The hBN-potential is
modified on the second graphene layer as to match a given configuration. For
example, in hBN-BG-hBN case 3 the potential on the upper layer is multiplied
by minus one and shifted by two third of the moiré cell diagonal.

Figure 1: The ‘hBN-potential’ $\Delta_{\boldsymbol{r}}$ acting on the carbon
sites at position $\boldsymbol{r}=(r_{x},r_{y})$, as adapted from Ref. [33].
The carbon sites are shown within one moiré unit cell with
$N_{\text{moir\'{e}}}=20$. Figure 2: Different stacking configurations of hBN-
BG. The second layer can be positioned such that at the corners of the moiré
unit cell (see Fig. 1) the stacking axis either goes through a nitrogen atom
(case i) or a boron atom (case ii). Figure 3: Different stacking
configurations of hBN-BG-hBN. In the cases 1–3 the upper hBN layer is parallel
to the lower hBN layer, in the cases 4–7 the upper hBN layer is rotated by
$60^{\circ}$.

After having properly defined the non-interacting model, we can now take into
account the Hubbard interaction and identify the possibility of magnetic order
in the material. We calculate the magnetic susceptibility $\chi^{m}$ in RPA by
summing up infinitely many diagrams via the exchange channel, analogous to
Ref. [39], resulting in

$\displaystyle\hat{\chi}^{m}=\hat{\chi}\left({\mathbbm{1}+U\hat{\chi}}\right)^{-1}$
(II.4)

Here, $\hat{\chi}$ is the polarization function

$\displaystyle\chi_{\boldsymbol{r},\boldsymbol{r^{\prime}}}(\omega,\boldsymbol{q})$
$\displaystyle=\frac{1}{\beta\Omega}\sum_{\omega^{\prime}}\int\text{d}\boldsymbol{k}\
g_{\boldsymbol{r},\boldsymbol{r^{\prime}}}(\omega^{\prime},\boldsymbol{k})g_{\boldsymbol{r^{\prime}},\boldsymbol{r}}(\omega+\omega^{\prime},\boldsymbol{k}+\boldsymbol{q}),$
(II.5)

a matrix of the size of the number of atoms in the moiré cell. The functions
$g$ are the bare green’s functions, $\beta=\frac{1}{T}$ is the inverse
temperature ($k_{\text{B}}=1$), $\omega$, $\omega^{\prime}$ are bosonic and
fermionic Matsubara frequencies, the integral is taken over wave vectors
$\boldsymbol{k}$ in the Brillouin zone with area $\Omega$, and
$\boldsymbol{q}$ is the wave vector of the corresponding magnetic spin density
wave. The magnetic susceptibility diverges when the ’Stoner denominator’, i.e.
the matrix $\left({\mathbbm{1}+U\hat{\chi}}\right)$ to be inverted, develops a
vanishing eigenvalue. This is the case for certain critical values of $T$ and
$U$. Since the interaction strength is a model parameter that is not known
exactly, we can calculate the critical interaction strength $U_{c}$ for fixed
temperature via the relation $U_{c}=-\frac{1}{x_{0}}$. Here, $x_{0}$ is the
most negative eigenvalue of the polarization function. $U_{c}$ marks the
minimum interaction strength needed for the system to be in a magnetic phase
at fixed temperature. The eigenvector $\boldsymbol{x}_{0}$ corresponding to
the eigenvalue $x_{0}$ is proportional to the magnetization on each atom,
whose magnitude cannot be obtained from this analysis. This is because all
other eigenvectors become irrelevant when the spin susceptibility diverges. As
in Ref. [39] we evaluate the Matsubara sum in $\chi$ numerically and set
$\omega=0$ for the static case and $\boldsymbol{q}=0$ for identical moiré
cells.

We also calculate the magnetic order in mean field theory to compare our
results with the RPA calculation. For this, we replace the spin density
operator
$S_{\boldsymbol{R},\boldsymbol{r}}^{z}=\frac{1}{2}(n_{\boldsymbol{R},\boldsymbol{r}\uparrow}-n_{\boldsymbol{R},\boldsymbol{r}\downarrow})$
by its mean field $S_{\boldsymbol{R},\boldsymbol{r}}^{z}\rightarrow\langle
S_{\boldsymbol{r}}^{z}\rangle+\delta S_{\boldsymbol{R},\boldsymbol{r}}^{z}$
and neglect higher orders of fluctuations $\mathcal{O}((\delta
S_{\boldsymbol{R},\boldsymbol{r}}^{z})^{2})$ [40]. Further, we neglect terms
in the Hamiltonian describing charge density waves, because we want to focus
on spin density fluctuations [31]. The resulting non-interacting Hamiltonian
(up to a constant)

$\displaystyle H$
$\displaystyle=\sum_{\begin{subarray}{c}\boldsymbol{r}_{1},\boldsymbol{r}_{2},\\\
\sigma,\boldsymbol{k}\end{subarray}}\Big{(}H_{\boldsymbol{r}_{1},\boldsymbol{r}_{2}}(\boldsymbol{k})-\sigma
U\left\langle
S_{\boldsymbol{r}}^{z}\right\rangle\delta_{\boldsymbol{r}_{1}\boldsymbol{r}_{2}}\Big{)}c^{\dagger}_{\boldsymbol{k},\boldsymbol{r}_{1}\sigma}c_{\boldsymbol{k},\boldsymbol{r}_{2}\sigma}$
(II.6)

is then solved self-consistently. This allows us to obtain mean field
estimates for interaction-induced energy gaps. These can then be compared with
spectral gaps that derive from the ’single-particle’-description, here
primarily from the substrate and stacking.

##  III Results

Before we analyze the magnetic order, it is insightful to first understand the
impact of the hBN-potential defined in Eq. (II) on the non-interacting system.
In Fig. 4 we compare the band structure of hBN-BG, with the bands of (free-
standing) BG downfolded into the mini Brillouin zone of the composite
structure.

Figure 4: Band structure of hBN-BG ($N_{\text{moir\'{e}}}=56$, red lines) and
free-standing BG (black lines) downfolded to the same mini Brillouin zone. At
the $K$ and $K^{\prime}$ points the hBN-potential (Eq. (II)) creates
secondary, gapped Dirac points at band energies
$\epsilon\sim+0.06\,\mathrm{eV}$ and $\epsilon\sim-0.05\,\mathrm{eV}$. In
addition, a band gap opens at charge neutrality (see Fig. 5 for a zoom-in).

In agreement with Ref. [14], secondary Dirac points form at the new mini
Brillouin zone edges and are split-up in hBN-BG. We observe only a small
particle hole asymmetry. In addition, a band gap opens around zero energy.

Figure 5: Low-enery band structure near the $K$-point of the mini Brillouin
zone. The band structure shows a band gap at charge neutrality for hBN-BG and
for hBN-BG-hBN cases 1–3. This band gap originates from an energy difference
at neighboring carbon sites. Although alignment of the hBN layer and
development of the moiré structure is necessary for the band gap opening, the
gap itself can be modelled in a small unit cell corresponding to $A=0$ in the
hBN-potential $\Delta_{\boldsymbol{r}}$.

The low-energy band structure mostly determines our results of the magnetic
order and is therefore shown as a zoom-in in Fig. 5. In this figure we see two
of the four Dirac points around $K$ in free-standing BG that are a result of
trigonal warping [36]. In both configurations of hBN-BG a band gap opens. For
the studied hBN-BG-hBN configurations the band structures split into two
subgroups. In the cases 1–3 the band gap opens further compared to the one-
sided cases whereas in the cases 4–7 the band gap closes again. The good
agreement within the subgroups is a result of the fact that the $C$ parameter
of the hBN-potential is mostly responsible for the band gap opening, while the
$A$ parameter does not have a significant influence. The $C$ parameter causes
an energy difference between A and B carbon sites in the cases 1–3 and an
energy difference between dimer and non-dimer sites in the cases 4–7 [41]. In
Fig. 3 we see that in the cases 1 –3, the upper hBN layer is parallel to the
lower layer (“parallel alignment”), whereas in the cases 4–7 the upper hBN
layer is rotated by $60^{\circ}$ (“$60^{\circ}$-alignment”).

Next, we turn to the possibility of magnetic order in our structures, within
the RPA approach described above. In all cases studied, the dominant magnetic
order is interlayer and intralayer antiferromagnetic, as shown in Fig. 6. No
noticeable variation of this antiferromagnetic order was found within the
moiré cell as the relevant eigenvectors turned out to be very homogeneous in
their magnitudes through the moiré cell. This is quite different from the
twisted bilayer cases studied, e.g. in Ref. [39, 42], where the eigenvectors
vary strongly in the moiré cell.

Figure 6: In all of our calculations the dominant magnetic order is found to
be antiferromagnetic with opposite signs of the magnetization at neighboring
carbon sites. Besides the relative spatial ordering pattern, the orientation
of the magnetization is not fixed as long as the spin-orbit interaction is not
taken into account.  Figure 7: Critical Hubbard interaction $U_{c}$ for
magnetic order as a function of the temperature. The different stacking
configurations differ in that for structures with a band gap (e.g. hBN-BG) at
temperatures smaller than the band gap, $U_{c}$ is larger as compared to
structures without a band gap (e.g. free-standing BG).

Now, we calculate $U_{c}$ for the different structures as a function of the
temperature in Fig. 7. We find that the low-energy band structures of the non-
interacting system fully determine the differences and the similarities
between free-standing BG, hBN-BG and hBN-BG-hBN. For free-standing BG and
$60^{\circ}$-aligned hBN-BG-hBN, structures with band crossings at the Fermi
energy, $U_{c}$ continuously increases with the temperature. This is expected
as thermal fluctuations at some point destroy the magnetic order. For
‘substrate-gapped’ structures such as hBN-BG or parallelly aligned hBN-BG-hBN,
$U_{c}$ differs from free-standing and $60^{\circ}$-aligned BG at temperatures
below the size of the band gap. In this regime, $U_{c}$ saturates at an
approximately constant value. This behavior is consistent with the
expectations that band structure gaps become visible at temperatures below
these energy scales. Similar as in the analysis of the band structure, we see
that the $A$ parameter of the hBN-potential does not play a role in our
results. This means that, at least at charge neutrality, the moiré length
scale does not influence the ordering tendency. In total, interaction effects
are weakened by the substrate-induced band gap, although the difference in
$U_{C}$ is small.

Next, we evaluate the magnetic order of free-standing BG, hBN-BG-hBN case 1,
and hBN-BG-hBN case 2 in mean field theory at a temperature of
$T=10^{-4}\,\mathrm{eV}$. In Fig. 8 we show the magnetization $\langle
S_{\boldsymbol{r}}^{z}\rangle$ in its magnitude and plot its absolute value
averaged over the moiré cell as a function of $U$. The critical interaction
strength $U_{c}$ needed to obtain a finite magnetization can be read off to be
$5.0\,\mathrm{eV}<U_{c}<5.2\,\mathrm{eV}$ for free-standing BG and hBN-BG-hBN
case 2, and $5.2\,\mathrm{eV}<U_{c}<5.4\,\mathrm{eV}$ for hBN-BG-hBN case 1.
These values of $U_{c}$ are in agreement with the RPA results described above.
Next, we calculate the band structure of the converged mean field Hamiltonian.
In Fig. 9 we extract the band gap widths at charge neutrality, again for
different values of the interaction strength. As commonly observed in Hubbard
models [40], for parameters with nonzero magnetization an interaction-induced
band gap is opened and we can see that the size of this band gap increases
with the magnetization. In the non-magnetic phase the structures are described
by the non-interacting picture, which explains the band gap for low values of
$U$ in hBN-BG-hBN case 1. At $U=5.2\,\mathrm{eV}$ this substrate-induced band
gap in case 1 is still larger than the interaction-induced band gap in free-
standing and hBN-BG-hBN case 2. At $U=5.4\,\mathrm{eV}$ the interaction-
induced band gaps of free-standing BG and hBN-BG-hBN case 2 exceed that of
hBN-BG-hBN case 1.

Figure 8: Staggered magnetization in self-consistent mean field theory as a
function of the Hubbard interaction $U$ at a temperature of
$T=10^{-4}\,\mathrm{eV}$.  Figure 9: Band gap sizes of the converged mean
field theory corresponding to Fig. 8 as a function of the Hubbard interaction
$U$. For low $U$, the model can be described by the non-interacting system and
only hBN-BG-hBN case 1 has a band gap. For large $U$ the magnetic state
results in larger ‘magnetic’ band gaps in all structures.

##  IV Discussion

Both, the RPA analysis and the mean field calculation show an upward shift of
$U_{c}$ for substrate-gapped structures at low temperatures. Quantitatively,
this shift is less than 6%, putting constraints on the possibility that
instabilities of BG may be observed in free-standing bilayers but not in the
heterostructures. Even in this small window, the systems that do not have an
interaction-driven ordered ground state show a substrate-induced gap. Thus, if
no gap is observed, one should deal with a case without substrate-induced gap
where the effective onsite repulsion is too small to trigger a ground state
change at the given temperature. Lowering the temperature could still open a
gap, but this gap might be too small to be observed. Note that the trigonal
warping terms split any quadratic band crossing points of simple bilayer band
structures into Dirac points, such that at zero temperature there should be a
nonzero threshold value in $U$ for ordering.

We also discuss the justification and magnitude of the onsite repulsion used
in the Hubbard modelling of this work. Clearly, a more realistic model should,
beyond the onsite term, also include non-local interaction terms. This non-
local interaction can approximately be incorporated in an effective onsite
repulsion. In Ref. [37] such a simplification was analyzed in ab-initio
calculations employing the constrained random phase approximation and it is
suggested to estimate the effective onsite interaction via $U\approx
U_{0}-U_{1}$, where $U_{0}$ is the ‘pure’ onsite interaction and $U_{1}$ is
the nearest-neighbor interaction. For free-standing BG a value of $U\approx
9.3\,\mathrm{eV}-4.7\,\mathrm{eV}=4.6\,\mathrm{eV}$ would be obtained [32].
However, the hBN-substrate may screen and thus reduce both, the ‘pure’ onsite
interaction and the nearest neighbor interaction as calculated in Ref. [32].
This screening from the substrate may approximately cancel out for the
effective $U$ and we get, again, $U\approx
8.3\,\mathrm{eV}-3.7\,\mathrm{eV}=4.6\,\mathrm{eV}$. In all these
approximations we stay well below our threshold values for ordering, i.e. on
the side where potential gaps would be induced by the hBN-layers. This is
however not consistent with the experimental reports of gaps in freely
suspended Bernal-stacked systems [17, 18, 19, 20, 21, 22, 23]. Thus, we
conclude that a quantitative description of these systems remains a challenge.

From the mean field calculation we find that, while the band gaps of single-
particle nature are less variable, the interaction-induced gap is a strongly
varying function of $U$ and largely independent of the putative non-
interacting band structure once the interaction threshold is exceeded. This
observation makes it difficult to compare our different materials, because the
hBN-substrated materials may have a different effective $U$ compared to free-
standing BG due to the additional screening effects. Our mean field results
imply that the precise value of $U$ will be the determining component for the
size of the magnetic band gap. Further, we see that hBN-BG and parellely
aligned hBN-BG-hBN are gapped for all values of $U$. A case where the
substrate destroys the ordering instability and the interaction-induced gap
without producing another single-particle gap does not appear in our study.

##  V Summary

In summary, we identified several stacking configurations and differentiated
between parallelly aligned and $60^{\circ}$-aligned hBN-BG-hBN. The influence
of hBN on the electronic structure of the BG was captured by an effective
single-particle potential in analogy to earlier DFT work [33]. hBN-BG and
parallelly aligned hBN-BG-hBN have a band gap of ‘single-particle’ origin
already in the non-interacting picture. For temperatures smaller than this
band gap, the critical Hubbard interaction $U_{c}$ required for a magnetic
ordering is larger than for the $60^{\circ}$-aligned hBN-BG-hBN cases and
free-standing BG. Thus interaction effects are weakened by the substrate-
induced band structure changes. However, the shift in $U_{c}$ is less than 6%.
Mean field theory allowed us to calculate the band gap sizes of the
interacting systems. Our results demonstrate that the interaction-induced band
gap is a strongly varying function of the onsite interaction $U$.

### Acknowledgements

The Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) is
acknowledged for support through RTG 1995, within the Priority Program SPP
2244 “2DMP” and under Germany’s Excellence Strategy-Cluster of Excellence
Matter and Light for Quantum Computing (ML4Q) EXC2004/1 - 390534769. We
acknowledge support from the Max Planck-New York City Center for Non-
Equilibrium Quantum Phenomena. Calculations were performed with computing
resources granted by RWTH Aachen University under project rwth0510.

## References

  * Yankowitz  _et al._ [2019] M. Yankowitz, Q. Ma, P. Jarillo-Herrero, and B. J. LeRoy, Nature Reviews Physics 1, 112 (2019).
  * Liu  _et al._ [2003] L. Liu, Y. P. Feng, and Z. X. Shen, Physical Review B 68, 104102 (2003).
  * McCann and Koshino [2013] E. McCann and M. Koshino, Reports on Progress in Physics 76, 056503 (2013).
  * Yankowitz  _et al._ [2012] M. Yankowitz, J. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod, and B. J. LeRoy, Nature Physics 8, 382 (2012).
  * Slotman  _et al._ [2015] G. J. Slotman, M. M. van Wijk, P.-L. Zhao, A. Fasolino, M. I. Katsnelson, and S. Yuan, Physical Review Letters 115, 186801 (2015).
  * Jung  _et al._ [2015] J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nature Communications 6, 1 (2015).
  * Song  _et al._ [2013] J. C. W. Song, A. V. Shytov, and L. S. Levitov, Physical Review Letters 111, 266801 (2013).
  * Moon and Koshino [2014] P. Moon and M. Koshino, Physical Review B 90, 155406 (2014).
  * Woods  _et al._ [2014] C. Woods, L. Britnell, A. Eckmann, R. Ma, J. Lu, H. Guo, X. Lin, G. Yu, Y. Cao, R. Gorbachev, _et al._ , Nature Physics 10, 451 (2014).
  * Ponomarenko  _et al._ [2013] L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods, J. Wallbank, _et al._ , Nature 497, 594 (2013).
  * Hunt  _et al._ [2013] B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Science 340, 1427 (2013).
  * Bokdam  _et al._ [2014] M. Bokdam, T. Amlaki, G. Brocks, and P. J. Kelly, Physical Review B 89, 201404 (2014).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Physical Review B 99, 075127 (2019).
  * Kim  _et al._ [2018] H. Kim, N. Leconte, B. L. Chittari, K. Watanabe, T. Taniguchi, A. H. Macdonald, J. Jung, and S. Jung, Nano Letters 18, 7732 (2018).
  * Endo  _et al._ [2019] K. Endo, K. Komatsu, T. Iwasaki, E. Watanabe, D. Tsuya, K. Watanabe, T. Taniguchi, Y. Noguchi, Y. Wakayama, Y. Morita, and S. Moriyama, Applied Physics Letters 114, 243105 (2019).
  * Kuiri  _et al._ [2021] M. Kuiri, S. K. Srivastav, S. Ray, K. Watanabe, T. Taniguchi, T. Das, and A. Das, Physical Review B 103, 115419 (2021).
  * Velasco Jr  _et al._ [2012] J. Velasco Jr, L. Jing, W. Bao, Y. Lee, P. Kratz, V. Aji, M. Bockrath, C. N. Lau, C. Varma, R. Stillwell, D. Smirnov, F. Zhang, J. Jung, and A. H. MacDonald, Nature Nanotechnology 7, 156 (2012).
  * Bao  _et al._ [2011] W. Bao, L. Jing, J. Velasco Jr, Y. Lee, G. Liu, D. Tran, B. Standley, M. Aykol, S. B. Cronin, D. Smirnov, M. Koshino, E. McCann, M. Bockrath, and C. N. Lau, Nature Physics 7, 948 (2011).
  * Freitag  _et al._ [2012] F. Freitag, J. Trbovic, M. Weiss, and C. Schönenberger, Physical Review Letters 108, 076602 (2012).
  * Veligura  _et al._ [2012] A. Veligura, H. J. van Elferen, N. Tombros, J. C. Maan, U. Zeitler, and B. J. van Wees, Physical Review B 85, 155412 (2012).
  * Grushina  _et al._ [2015] A. L. Grushina, D.-K. Ki, M. Koshino, A. A. L. Nicolet, C. Faugeras, E. McCann, M. Potemski, and A. F. Morpurgo, Nature Communications 6, 6419 (2015).
  * Nam  _et al._ [2018] Y. Nam, D.-K. Ki, D. Soler-Delgado, and A. F. Morpurgo, Science 362, 324 (2018).
  * Myhro  _et al._ [2018] K. Myhro, S. Che, Y. Shi, Y. Lee, K. Thilahar, K. Bleich, D. Smirnov, and C. N. Lau, 2D Materials 5, 045013 (2018).
  * Zhang  _et al._ [2011] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. MacDonald, Physical Review Letters 106, 156801 (2011).
  * Throckmorton and Vafek [2012] R. E. Throckmorton and O. Vafek, Physical Review B 86, 115447 (2012).
  * Kharitonov [2012] M. Kharitonov, Physical Review B 86, 195435 (2012).
  * Nilsson  _et al._ [2006] J. Nilsson, A. H. Castro Neto, N. M. R. Peres, and F. Guinea, Physical Review B 73, 214418 (2006).
  * Lang  _et al._ [2012] T. C. Lang, Z. Y. Meng, M. M. Scherer, S. Uebelacker, F. F. Assaad, A. Muramatsu, C. Honerkamp, and S. Wessel, Physical Review Letters 109, 126402 (2012).
  * Scherer  _et al._ [2012a] M. M. Scherer, S. Uebelacker, and C. Honerkamp, Physical Review B 85, 235408 (2012a).
  * Scherer  _et al._ [2012b] M. M. Scherer, S. Uebelacker, D. D. Scherer, and C. Honerkamp, Physical Review B 86, 155415 (2012b).
  * Honerkamp [2008] C. Honerkamp, Physical Review Letters 100, 146404 (2008).
  * Rösner  _et al._ [2015] M. Rösner, E. Şaşıoğlu, C. Friedrich, S. Blügel, and T. O. Wehling, Physical Review B 92, 085102 (2015).
  * Sachs  _et al._ [2011] B. Sachs, T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein, Physical Review B 84, 195414 (2011).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano Letters 10, 804 (2010).
  * Moon and Koshino [2012] P. Moon and M. Koshino, Physical Review B 85, 195458 (2012).
  * Rozhkov  _et al._ [2016] A. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, Physics Reports 648, 1 (2016).
  * Schüler  _et al._ [2013] M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Physical Review Letters 111, 036601 (2013).
  * Tang  _et al._ [2015] H.-K. Tang, E. Laksono, J. N. B. Rodrigues, P. Sengupta, F. F. Assaad, and S. Adam, Physical Review Letters 115, 186602 (2015).
  * Klebl and Honerkamp [2019] L. Klebl and C. Honerkamp, Physical Review B 100, 155145 (2019).
  * Claveau  _et al._ [2014] Y. Claveau, B. Arnaud, and S. D. Matteo, European Journal of Physics 35, 035023 (2014).
  * Mucha-Kruczyński  _et al._ [2010] M. Mucha-Kruczyński, E. McCann, and V. I. Fal'ko, Semiconductor Science and Technology 25, 033001 (2010).
  * Klebl  _et al._ [2020] L. Klebl, D. M. Kennes, and C. Honerkamp, Physical Review B 102, 085109 (2020).