Multi-flat bands and strong correlations in Twisted Bilayer Boron Nitride

TL;DR

This paper demonstrates that twisted bilayer boron nitride naturally exhibits multiple flat bands without precise angle tuning, leading to correlated insulating and superconducting phases, making it a promising platform for studying two-dimensional correlation physics.

Contribution

It reveals that TBBN hosts multiple flat bands and correlated phases without fine-tuning, offering a more robust alternative to TBG for exploring strongly correlated phenomena.

Findings

Multiple flat bands emerge in TBBN without fine-tuning.

Doping TBBN leads to correlated insulating and superconducting states.

TBBN exhibits well-separated, degenerate bands within the gap.

Abstract

In a groundbreaking experimental advance it was recently shown that by stacking two sheets of graphene atop of each other at a twist angle close to one of the so called "magic angles", an effective two-dimensional correlated system emerges. In this system the kinetic energy of the low-energy electrons is much reduced and consequently interactions become very relevant, providing a new platform into the physics of two-dimensional correlated materials. Evidence of a proposed Mott insulating as well as superconducting state in these highly tunable systems has spurred much attention as they could pave the way to understanding long-standing questions of high-$T_c$ superconductivity or provide candidate systems for topological chiral superconductors; key to highly relevant quantum technologies. Here, we demonstrate that twisted bilayer boron nitride (TBBN) is an exciting and even richer alternative to twisted bilayer graphene (TBG). Crucially, we show that in TBBN multiple flat bands emerge without having to fine tuning close to a "magic angle" that upon doping lead to correlated phases of matter (insulating and superconducting). TBBN could thus be much less sensitive to small deviations in the twist angle and therefore provide a particularly suited experimental platform to study correlation physics in two dimensions. Furthermore, we find that in marked contrast to TBG at small twist angle families of 2,4 and 6-fold degenerate, well separated, bands emerge within the gap, considerably broadening the addressable physics.

Figures

Equations

$$U=\frac{1}{N}\sum _m⟨mm|V_{ee}|mm⟩,$$

$$J=\frac{1}{N(N-1)}\sum _{m\ne m^{\prime }}⟨m{m}^{\prime }|V_{ee}|m^{\prime }m⟩,$$

$$U^{\prime }=\frac{1}{N(N-1)}\sum _{m\ne m^{\prime }}⟨mm|V_{ee}|m^{\prime }{m}^{\prime }⟩,$$

$$\bar{U}=\frac{{\sum }_m{\sum }_{\sigma }{\bar{P}}_{mm}^{\sigma }{\bar{P}}_{mm}^{-\sigma }(mm|mm)}{{\sum }_m{\sum }_{\sigma }{n}_{mm}^{\sigma }{n}_{mm}^{-\sigma }},$$

$$\overline{U^{\prime }}=\frac{{\sum }_{\{m\}}^{\ast }{\sum }_{\sigma {\sigma }^{\prime }}{\bar{P}}_{m{m}^{\prime }}^{\sigma }{\bar{P}}_{m^{\prime }{m}^{\prime \prime }}^{{\sigma }^{\prime }}(m{m}^{\prime }|m^{\prime \prime }{m}^{\prime \prime \prime })}{{\sum }_{m\ne m^{\prime }}{\sum }_{\sigma }[n_{mm}^{\sigma }{n}_{m^{\prime }{m}^{\prime }}^{\sigma }+n_{mm}^{\sigma }{n}_{m^{\prime }{m}^{\prime }}^{-\sigma }]},$$

$$\bar{J}=\frac{{\sum }_{\{m\}}^{\ast }{\sum }_{\sigma }{\bar{P}}_{m{m}^{\prime }}^{\sigma }{\bar{P}}_{m^{\prime }{m}^{\prime \prime }}^{\sigma }(m{m}^{\prime \prime \prime }|m^{\prime \prime }{m}^{\prime })}{{\sum }_{m\ne m^{\prime }}{\sum }_{\sigma }{n}_{mm}^{\sigma }{n}_{m^{\prime }{m}^{\prime }}^{\sigma }},$$

$$H^0=-t\sum _{⟨i,j⟩}\sum _{b=1}^2\sum _{\sigma =\uparrow ,\downarrow }{c}_{i,b,\sigma }^{†}{c}_{j,b,\sigma }-2t\sum _i\sum _{b=1}^2\sum _{\sigma =\uparrow ,\downarrow }{n}_{i,b,\sigma },$$

$$H_i^U=$$

$$-J^{\prime }\left(c_{i,1,\downarrow }^{†}{c}_{i,2,\uparrow }^{†}{c}_{i,2,\downarrow }{c}_{i,1,\uparrow }+c_{i,2,\uparrow }^{†}{c}_{i,2,\downarrow }^{†}{c}_{i,1,\uparrow }{c}_{i,1,\downarrow }+H.c.\right)$$

$$\sum _{k,q}={\tilde{\Gamma }}^{{\Lambda }_C}(k,q)[{\hat{\Delta }}_k,{\hat{\Delta }}_q]$$

Full text

Multi-flat bands and strong correlations in Twisted Bilayer Boron Nitride

Lede Xian

These authors contributed equally

Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany

  

Dante M. Kennes

These authors contributed equally

Dahlem Center for Complex Quantum Systems and Fachbereich Physik, Freie Universität Berlin, 14195 Berlin, Germany

  

Nicolas Tancogne-Dejean

Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany

  

Massimo Altarelli

Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany

  

Angel Rubio

angel.rubio@mpsd.mpg.de

Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149, 22761 Hamburg, Germany

Center for Computational Quantum Physics (CCQ), The Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

(March 17, 2024)

**Effectively two-dimensional materials rank among the most fascinating systems of contemporary condensed matter research and beyond. They can exhibit interesting phases of matter, such as the quantum hall state Stern (2008), are vital in quantum computing, due to their unique braiding properties and can host anyons Ivanov (2000); Nayak et al. (2007), which fundamentally refines our classification of particles in fermions and bosons. They, furthermore, still pose a grand challenge to theoretical tools; which lack a flagship method describing two-dimensional correlated systems. This has, e.g., given rise to the field of quantum simulations of two-dimensional systems in cold gases Kondov et al. (2015); Gross and Bloch (2017). In a groundbreaking experimental advance Cao et al. (2018a, b); Yankowitz et al. (2018) it was recently shown that by stacking two sheets of graphene atop of each other at a twist angle close to one of the so called “magic angles” Bistritzer and MacDonald (2011); Tarnopolsky et al. (2018), an effective two-dimensional correlated system emerges. In this system the kinetic energy of the low-energy electrons is much reduced and consequently interactions become very relevant, providing a new platform into the physics of two-dimensional correlated materials. Evidence of a proposed Mott insulating as well as superconducting state in these highly tunable systems has spurred much attention Cao et al. (2018a, b) as they could pave the way to understanding long-standing questions of high-$T_{c}$ superconductivity Anderson (2013) or provide candidate systems for topological chiral superconductors; key to highly relevant quantum technologies Stern (2008). Here, we demonstrate that twisted bilayer boron nitride (TBBN) is an exciting and even richer alternative to twisted bilayer graphene (TBG). Crucially, we show that in TBBN multiple flat bands emerge without having to fine tuning close to a “magic angle” that upon doping lead to correlated phases of matter (insulating and superconducting). TBBN could thus be much less sensitive to small deviations in the twist angle and therefore provide a particularly suited experimental platform to study correlation physics in two dimensions. Furthermore, we find that in marked contrast to TBG at small twist angle families of 2,4 and 6-fold degenerate, well separated, bands emerge within the gap, considerably broadening the addressable physics. **

Graphene is considered a very special 2D material, as its symmetries give rise to gapless linear dispersions for the low-energy bands around the Fermi level at the corner of the Brillouin zone Suárez Morell et al. (2010); Kim et al. (2017). These low energy states can be well described by a massless Dirac Hamiltonian, which gives rise to many of graphene’s unusual properties, such as lack of back-scattering and Klein tunneling Neto et al. (2009). So far, it is unclear whether this special gapless and massless features are crucial in the observed strongly-correlated-like phenomena Cao et al. (2018a, b); Yankowitz et al. (2018) in twisted bilayer systems Wu et al. (2018a). Here, we elaborate on this question by studying a similar system, twisted bilayer boron nitride (TBBN), which however has a large band gap. Hexagonal boron nitride (hBN) is also called “white graphene”, because it has similar hexagonal structure as graphene Wirtz et al. (2006); Blase et al. (1995). While graphene has two identical carbon atoms in the unit cell, they are replaced by one nitrogen and one boron atoms in hBN. This difference gives rise to a large band gap and finite effective mass for the low-energy electrons in hBN. As a consequence, hBN is widely used as an ideal insulator in 2D electronic devices (and to enhance the properties of graphene by isolating it from its environment). Interestingly, our density functional theory (DFT) calculations reveal that flat bands develop at both the top of the valence bands and bottom of the conduction bands in twisted bilayer hBN with a small twist angle. In marked contrast to the TBG systems, the band width of these flat bands decrease monotonically with twist angle and there is no occurrence of “magic angles” Bistritzer and MacDonald (2011). We further evaluate the Coulomb interaction for those localized states in the flat bands using a self-consistent ab inito DFT+$U$ Tancogne-Dejean et al. (2017) and functional renormalization group (FRG) Metzner et al. (2012) method to show that the flat band at the top of the valence bands can host exotic strongly-correlated physics, such as the appearance of a Mott insulator phase and unconventional superconductivity for slight amount of doping or photo-doping. Those results provide strong evidence that TBBN could also be utilized to study unconventional superconductors, much like TBG. Our results suggest that the gapless-ness and massless-ness of graphene is not essential in the observed exotic strongly-correlated phenomena motivating studies of many more two-dimensional materials stacked with a small twist angle.

Starting from a regular bilayer hBN with stacking as in the bulk, TBBN can be formed by twisting the top or the bottom layer with angle $\theta$. As shown in Fig. 1(a) and (b), there are two possible, distinct configurations of TBBN, called, $\alpha$ and $\beta$, that lead to the same supercell. Their twist angles differ by $60^{\circ}$. In TBBN, there are regions where nitrogen atoms on the top layer are approximately on top of those of the bottom layer (the N/N region, as highlighted by dash red circles) and regions where the same holds for the boron atoms (the B/B region, as highlighted by solid blue circles). In configuration $\alpha$, the two regions are located at different sites of the supercell, while in configuration $\beta$, they are located at the same sites (see Fig. 1(a) and (b)).

To study the effect of twisted stacking, we first investigate the band structures of TBBN with fixed interlayer separation $d$, which is set to be the optimized interlayer distance of bilayer hBN without twist (3.23 Angstrom, that is very close to the experimental value). As shown in Fig. 1(c) and (f), similar to the case of TBG, flat bands appear at the top of the valence band in TBBN at small twist angles. These flat band structures are significantly different from those of monolayer and normal bilayer hBN shown in Fig. 1(d) and (e), respectively. In contrast to the case of graphene, the flat bands are isolated from the other bands with a relatively large gap of over 100meV, and intriguingly, more than one set of them develop as twist angle decreases. For the smallest twist angle we studied in (Fig. 1(f)), three set of flat bands can be clearly identified (see the labels on the right side of Fig. 1(f)) and the number of bands in each set is 2, 4 and 6. The band width of these flat bands decreases significantly with the twist angle. As shown in Fig. 1(g), the band width of the top of the valence bands (the first set of flat bands) decreases from 101.4 meV at $7.34^{\circ}$ to 1.6 meV at $2.64^{\circ}$. Modeling the top of the valence bands of hBN by a massive Dirac Hamiltonian suggests that indeed the band width of these flat bands will decrease monotonically with twist angles without the existence of any “magic angles” (see supplemental materials). Similarly, flat bands also appear at the bottom of the conduction bands as the twist angle decreases (see Fig. 1(f) at around $1.87\sim 1.90$ eV).

The calculated DFT band structures for the two configurations $\alpha$ and $\beta$ of TBBN with fixed interlayer separation are almost identical, showing only subtle differences, while they are noticeably different after relaxation. To further understand the effect of the structural relaxation in the two configurations and the nature of the flat bands, we plot the charge density distribution for the low-energy states in the two configurations of TBBN in Fig. 2. Consistent with their narrow band width, those low-energy states are well localized in real space, with different characters in different flat bands. As expected, the low-energy states at the top of the valence bands are dominated by nitrogen $p_{z}$ orbital and those at the bottom of the conduction bands are dominated by boron $p_{z}$ orbital (see Fig. 2(g) and (f)). Therefore, the flat-band states of the valence bands are localized around the N/N region, while those of the conduction bands are localized around the B/B region. The relaxation increases the interlayer separation in the N/N region for both configurations and the valence bands of the two configurations remain very similar. However, we found that the response to the relaxation in the B/B regions is different for the two configurations, because of their different relative position with respect to the N/N region. This difference leads to different energy ordering of these flat band states in the conduction bands for the two configurations. Along with the different response to the relaxation, the band gaps of the two configurations are also different. In the total density of states, shown in Fig. 2(d) and (f), the low energy flat bands manifest as a few sharp peaks near the band edges, which are significantly different from those of the monolayer and untwisted bilayer hBN.

Next, in order to characterize the electronic correlations in those flat bands (that tend to be beyond the capabilities of standard local and semi-local XC functionals) we investigate TBBN by means of self-consistent DFT+$U$ (see Methods). We aim at extracting the most ab initio information possible concerning the electron-electron correlations in TBBN, that afterwards will be used to motivate the functional renormalization group (FRG) treatment (see below). We focus on the bands labeled 3 in Fig. 2 hereafter. These are two-fold degenerated bands, well described by a two-band Hubbard model. Concentrating on local interactions in such a model, one needs to determine an on-site Hubbard $U$ in addition to the Hunds-coupling $J$. Assessing how important $J$ is, from ab initio calculations, is crucial to our model building of TBBN.

DFT+$U$ is routinely used to describe transition metal oxides. A study of TBBN is quite different, as the localization occurs at a much larger spatial scale, and the effective electronic parameters are expected to be much smaller. Extracting quantitative parameters at finite doping is currently numerically not feasible and we considern just the undoped scenario. Notwithstanding, in order to address the role of $J$ and $U$ at least in this simpler case, we use the recently proposed ACBN0 functional, which allows for a DFT-based consistent ab-initio calculation of the value of $U$ and $J$ (see Methods section). For each system, from $7.34^{\circ}$ to $2.64^{\circ}$, we constructed the “single shot” Wannier states corresponding to the flat bands, and then computed the effective $U_{\mathrm{eff}}=U-J$. As shown in Fig. 3(a), we obtained an $U_{\mathrm{eff}}$ that scales linearly with the twist angle, similar to what was proposed for TBG Cao et al. (2018b). The calculated scaling of $U_{\mathrm{eff}}$ is 34.9 meV/degree. Compared with the value of about 5 to 26 meV/degree proposed for TBG Cao et al. (2018b), this relatively large value suggests a quite strong electron-electron correlation effect in TBBN, probably due to the smaller dielectric constant in hBN.

In order to asses the role of $J$, we also computed the so-called Hubbard-Kanamori parameters $U$, $U^{\prime}$, and $J$, using the bare Coulomb interaction (see Methods section). Our results yield two important results: i) we found that $(U-U^{\prime})/2$ is always very close to $J$ (which is a relation strictly motivated only for the case of cubic symmetry for $d$ orbitals Werner and Millis (2007)). ii) the ratio of $J/U$ is found to be small (between 0.03 to 0.05 for the various angles). These results are consistent with previous assumptions made for TBG, but are here confirmed using ab initio simulations. At finite doping, screening will decrease the values of the interaction parameters, which unfortunately is beyond an ab initio analysis at present. We thus next, consider the effect of correlations, keeping $U$ and $J$ as free parameters, but assuming $J=(U-U^{\prime})/2$.

We now turn to the construction of an effective low-energy model to deal with the inclusion of correlations, which can drive the system into different phases of matter. We find that the dispersion relation of the bands labeled by 3 and 7 in Fig. 2 are fitted well by a two-band nearest neighbor hopping $t$ tight binding model on a triangular lattice (in agreement with the corresponding charge density plots in Fig 2). In particular, this fit works well irrespective of the twist angle. This implies that the twist angle only appears as an indirect parameter of our model, via $U$, $J$, and the hopping $t$ (which is directly taken as $1/9$ of the DFT bandwidth). Therefore, analyzing the results in terms of $U/t$, $J/t$ and band-filling covers all the twist angles at once.

The role of Coulomb interaction versus phonons in TBBN is unclear at the moment. If phonons dominate we would expect an instability to emerge in the paring channel driving the system into a (BCS-type) s-wave superconducting phase, similar to what has been proposed for TBG Lian et al. (2018); Choi and Choi (2018); Wu et al. (2018b). If on the other hand Coulomb repulsion is dominant, the fragile balance of competing orders gives rise to a zoo of different possible phases in dependence of interaction strength and filling Kennes et al. (2018); Kang and Vafek (2018); Yuan and Fu (2018); Koshino et al. (2018); Su and Lin (2018); Ochi et al. (2018); Chittari et al. (2019). To explore the later scenario, we employ a FRG approach (see Methods), which can be viewed as a renormalization group enhanced random phase approximation treatment of the interacting system. To this end we include an onsite Hubbard-type of repulsion $U$ as well as a weak Hund’s coupling Xu and Balents (2018) term $J$ (in agreement with the above DFT+$U$ analysis). The dispersion relation in the full Brillouin zone is illustrated in a false color plot in Fig. 3 (e). We define the chemical potential $\mu$ with respect to the van Hove filling (meaning that at $\mu=0$ the density of states exhibits a van Hove singularity cf. Fig. 2 (d) and (f)). At this filling the Fermi surface also exhibits nesting giving rise to a competition between a tendency for $d\pm id$ superconductivity and a spin density wave as the interaction is turned on. This competition can be studied in an unbiased fashion using the FRG for small to intermediate $U/t$ by considering the so-called effective two-particle interaction $\Gamma$. During the renormalization group flow this effective two-particle interaction flows to strong coupling (diverges), which can be seen as the onset of ordering in the system. The energy scale at which this divergence occurs can be identified roughly with the critical temperature scale $T_{c}$ of the corresponding order. The momentum structure of the effective two-particle interaction indicates the type of ordering. We here parameterize the three independent momenta of the effective two-particle interaction $\Gamma(k_{1},k_{2},k_{3})$ by their projection on the Fermi surface as well the angles $\phi_{i}$ of this projection. The calculated FRG phase diagram is summarized in Fig. 3 (c) for zero Hund’s coupling $J/t=0$. Turning on small Hund’s couplings of up to 18% changes $T_{c}$ as shown in Fig. 3 (d), but does not influence the phase-diagram qualitatively, reinforcing the idea that the Hund’s coupling $J$ does not play a major role here. Close to van Hove filling ($\mu=0$) we find topological $d\pm id$ superconductivity in this model for not too large $U/t$ whereas it turns into a spin-density wave (SDW) for larger values of $U/t$ Thomson et al. (2018). Further away from $\mu=0$ we also find metallic behavior which means that up to the point in energy space where the flow was stopped no divergence was encountered in the effective two-particle interaction. Fig. 3 (b) illustrate the momentum structure of the effective two-particle interaction in the phase of $d\pm id$ superconductivity. The dominant diagonal features with sign change (blue to red) indicate that $d\pm id$ superconductivity is the dominant ordering tendency in doped TBBN.

Finally, we note that our ab initio results show a linear dependence of $U$ in the twist angle, and a polynomial (order 2) variation of the bandwidth (and hence of $t$) in the twist angle for angle $>2^{\circ}$ (see Fig. S3 in the supplemental materials). These two results implies that the ratio $U/t$ decreases as the twist angle increases and the system can be mechanically tuned through a transition from the SDW to the d-wave superconducting phase by tuning the twist angle.

In sum, here we propose twisted bilayer boron nitride as a fascinating material in which a zoo of different electronic phenomena are expected. If the twist angle between the layers is small, groups of two-, three- and four-fold degenerate and energetically separated flat-bands show up. These do not rely crucially on how close the chosen angle is to a “magic” angle as in the recently studied twisted bilayer graphene. Furthermore, at small angles the flat bands are truly separated from the valence and conduction bands, rendering twisted bilayer boron nitride an ideal candidate to study the physics of strong correlations in these materials. A combined DFT+$U$ and FRG ansatz reveals that the physics of the group of bands that are two-fold degenerate should be very similar to that of twisted bilayer graphene, possibly supporting topological d-wave superconductivity upon hole (or electron) doping. The existence of multiple flat bands significantly enriches the physics that can be probed. E.g. the separation of the flat bands found is in the 100meV regime. Coupling the present findings with the excitement of hyperbolic phonons Dai et al. (2014) opens the possibility of having a strong phonon mediated coupling between the flat bands as the distance between them is in quasi-resonance with these phonon modes. Photodoping the system is another particularly intriguing avenue of future research. We predict that using photodoping the identified correlated phases (insulator and superconducting ones) can be induced in a transient, non-equilibrium manner, strongly motivating experimental advances along these lines. Our study also opens the door to similar theoretical work in other 2D materials.

Acknowledgements

Useful discussions with Martin Claassen and Mei-Yin Chou are acknowledged. This work was supported by the European Research Council (ERC-2015-AdG694097) and Grupos Consolidados (IT578-13). The Flatiron Institute is a division of the Simons Foundation. LX acknowledges the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 709382 (MODHET).

Competing financial interests: The authors declare no competing financial interests.

Methods

Details about DFT Treatment

The ground state DFT calculations are performed with the Vienna Ab initio simulation package (VASP) Kresse and Hafner (1993). Plane waves are employed as a basis with an energy cutoff of 400 eV. The pseudo potentials are generated with the projector augmented wave method (PAW) Blöchl (1994) and the exchange-correlation potentials are treated in the local density approximation (LDA) Perdew and Zunger (1981). As the system size is relatively large, a 1x1x1 kgrid is employed in the ground state and relaxation calculations. Lattice constants are chosen such that they correspond to a lattice constant of 2.512 Angstrom for a 1x1 unit cell of hBN. During the relaxation, all the atoms are relaxed until the force on each atom is less than 0.01 eV/Angstrom. The VESTA code Momma and Izumi (2011) is used for the visualization of the charge density distributions of the low-energy states of TBBN.

DFT+$U$ calculations

The ab initio evaluation of the effective $U_{\mathrm{eff}}$ and the DFT+$U$ calculations are performed using the Octopus code Marques et al. (2003); Castro et al. (2006); Andrade et al. (2015); Tancogne-Dejean et al. (2017). A real-space spacing of 0.45 Bohr is chosen, and we employ norm-conserving pseudopotentials. The LDA is used for describing the local DFT part, and we utilize the ACBN0 functional to evaluate self consistently the effective $U_{\mathrm{eff}}$ of DFT+$U$. The vacuum size, atomic coordinates and lattice constant are taken to be the same as for the DFT treatment, as described above. Again, only the Gamma point is considered here.

The localized subspace is constructed as a single-shot Wannier states, taking the flat bands from LDA calculations. As the flat bands are energetically separated from the other occupied bands, this implies that the Wannier states reduce to the LDA states at the Gamma point.

The Hubbard-Kanamori parameters $U$, $U^{\prime}$ and $J$ are commonly defined from the expression

[TABLE]

where $V_{ee}$ denotes the screened Coulomb interaction.

To evaluate these parameters, we extended the definition of the ACBN0 function which yields

[TABLE]

where $n_{mm^{\prime}}^{\sigma}$ denotes the density matrix of the localized subspace, $\bar{P}_{mm^{\prime}}^{\sigma}$ is the the renormalization density matrix (see Ref. Tancogne-Dejean et al., 2017 for details), and $(mm^{\prime}|m^{\prime\prime}m^{\prime\prime\prime})$ refers to the Coulomb integrals computed from the bare Coulomb interaction. In the sum over the orbitals $\sum_{\{m\}}^{*}$, the asterisk means that the sum goes over all the orbitals (two here) except for the case in which $m=m^{\prime}=m^{\prime\prime}=m^{\prime\prime\prime}$. Importantly, we do not use any symmetry for computing the Coulomb integrals, and evaluate them directly on the real space grid, in order to analyze if the relation $\bar{U}-\bar{U^{\prime}}=2\bar{J}$ is fulfilled or not.

We tested both the fully-localized-limit and the around-mean-field double counting terms, and found no sizable change in the value of $U_{\mathrm{eff}}$.

Treating Correlations

To treat correlations we introduce a tight binding model motivated from the DFT analysis. We want to concentrate on the bands labeled by 3 and 7 in Fig. 1. These can be approximated very well by a two band nearest neighbor tight binding model on a triangular lattice

[TABLE]

where $c_{i,b,\sigma}^{(\dagger)}$ annihilates (creates) a particle on site $i$ in the band $b$ and with spin $\sigma$ and $n_{i,b,\sigma}=c_{i,b,\sigma}^{\dagger}c_{i,b,\sigma}$. The second term (energy-shift) ensures that we measure the chemical potential with respect to van Hove filling for convenience.

We add local Coulomb repulsion and Hund’s coupling to this by including

[TABLE]

in the full Hamiltonian $H=H^{0}+\sum_{i}H^{U}_{i}$. To reduce the number of parameters we set $J=J^{\prime}$ and $U^{\prime}=U-2J$ (which strictly can be shown to hold for d-orbitals in free space Werner and Millis (2007)). From our DFT+$U$ estimates we find that this relation is also approximately fulfilled for twisted boron nitride although there seems to be no symmetry argument why this has to be the case. We choose $U$ and $J$ as the two independent parameters to vary in the following, which fixes $U^{\prime}$.

The equilibrium phases of the Hamiltonian $H$ can be analyzed in the regime of small to intermediate interaction strength by the well-established functional renormalization group approach Metzner et al. (2012). Similar in spirit to other renormalization groups schemes within this method high-energy degrees of freedom are successively integrated out, resulting in a renormalized, effective low-energy theory of the problem. In our approach we monitor the behavior of the effective two-particle interaction over the flow, which addressed the successive inclusion of high-energy processes into an appropriate effective low energy model. Precursors to ordering tendencies are then indicated by a divergence of the effective two-particle interaction, signaling an instability of the Fermi-surface and a flow to strong coupling.

It is a particular advantage of the functional renormalization group that the method solely takes the Hamiltonian as an input and no a priori identification of dominant interaction channel (e.g. magnetic, superconducting, …) needs to be made by hand. This allows to treat these different channels on equal footing and thus analyze their competition in an unbiased fashion. As energy scales are successively integrated out over the flow, approaching the chemical potential $\mu$ the effective two-particle interaction acquires a strong momentum dependence. It is this momentum dependence combined with the analysis of the diverging channel that allows to draw conclusions about the dominant ordering in the system (varying the parameters of the Hamiltonian). Within the implementation of the FRG we use the implementation of Refs. Metzner et al. (2012); Kennes et al. (2018) we start the renormalization group flow from the bare two-particle interaction of the Hamiltonian $H^{U}$, which shows no momentum dependence (owed to its locality). Therefore, initially ($\Lambda_{i}$ denoting the start of the flow) the effective two particle interaction $\Gamma^{\Lambda_{i}}(k_{1},k_{2},k_{3},\omega_{1},\omega_{2},\omega_{3})$, which in general depends on three independent momenta and frequencies, is momentum and frequency independent. The functional renormalization group approach then provides a recipe of how $\Gamma^{\Lambda}(k_{1},k_{2},k_{3},\omega_{1},\omega_{2},\omega_{3})$ at a different value of the flow parameter $\Lambda$ can be obtained along the line from $\Lambda_{i}$ all the way to $\Lambda_{e}$, which is the end of the flow where the effective low-energy theory of the initial model can be extracted. This is done in terms of rather cumbersome differential equation given in Ref. Metzner et al. (2012). When an ordering tendency is encountered during this successive flow the effective two-particle interaction $\Gamma^{\Lambda}(k_{1},k_{2},k_{3},\omega_{1},\omega_{2},\omega_{3})$ diverges at a finite $\Lambda_{c}$ and the flow needs to be stopped. This $\Lambda_{c}$ can roughly be identified with the critical temperature $T_{c}$ associated with the ordering tendencies.

During the flow (employing momentum and frequency conservation) the effective two-particle interaction depends on three independent momenta as well as three independent frequencies (it is associated with four Fermionic operators, like the bare two-particle interaction). To make the numerical calculation feasible we have to reduce these degrees of freedoms. Since we are interested in the low-energy (equilibrium) physics of the model we project all frequencies to the chemical potential and all momenta to the Fermi-surface, but keeping their angle dependence, which we discretize in 18 uniform steps from [math] to $2\pi$. For the equilibrium low-energy physics of a system, processes far away from the Fermi-surface are irrelevant, while the angle dependence of the on-shell processes are relevant; justifying this approximation. Therefore in our formulation of the functional renormalization group the low-energy effective two-particle interaction on the Fermi surface is parametrized by the three angles $\phi_{1},\phi_{2},\phi_{3}$ of momenta $k_{1},k_{2},k_{3}$ and thus it can be equivalently written as $\Gamma(\phi_{1},\phi_{2},\phi_{3})$.

When a divergence is encountered in a particular ordering channel of the effective two-particle interaction, with corresponding order parameter $\hat{\Delta}_{k}$, one can employ a mean-field decomposition

[TABLE]

of the dominant contribution to the vertex after the flow. The prefactor $\tilde{\Gamma}^{\Lambda_{C}}(k,q)$ can then be decomposed into the irreducible representation of the underlying lattice model to find the symmetry (e.g. s,p,d,…) of the dominant ordering.