Engineering Three-Dimensional Moiré Flat Bands

Twisting two adjacent layers of van der Waals materials with respect to each other can lead to flat two-dimensional electronic bands which enables a wealth of physical phenomena. Here, we generalize this concept of so-called moiré flat bands to engineer flat bands in all three spatial dimensions controlled by the twist angle. The basic concept is to stack the material such that the large spatial moiré interference patterns are spatially shifted from one twisted layer to the next. We exemplify the general concept by considering graphitic systems, boron nitride, and WSe2, but the approach is applicable to any two-dimensional van der Waals material. For hexagonal boron nitride, we develop an ab initio fitted tight binding model that captures the corresponding three-dimensional low-energy electronic structure. We outline that interesting three-dimensional correlated phases of matter can be induced and controlled following this route, including quantum magnets and unconventional superconducting states.


S I. COMPUTATIONAL DETAILS FOR 3D TWISTED GRAPHENE, WSE2 AND BORON NITRIDE
For the calculations of 3D twisted graphene, we construct the unit cell with a twisted double bilayer graphene at twist angles close to 0 degree, and impose periodic boundary condition along all three dimensions. As it is not realistic to optimize such a large system with density functional theory (DFT) calculations, we fix the lattice constant along the out-of-plane direction to be 13.415 Å, and set the in-plane lattice constant according to the twist angles such that it corresponds to 2.46 Å for a 1x1 cell. The atomic structure is relaxed using the LAMMPS code [1] with the same parameters as described in [2]. The intralayer interactions within each graphene layer are modeled via the second-generation reactive empirical bondorder (REBO) potential [3]. The interlayer interactions are modeled via the Kolmogorov-Crespi (KC) potential [4], using the recent parametrization of [5]. The relaxation is performed using the fast inertial relaxation engine (FIRE) algorithm [6].
We calculate the band structures for 3D twisted graphene using the tight-binding parametrization proposed in Ref. [7] Here, the operator c ( †) i annihilates (creates) an electron in the p z orbital of the carbon atom at site r i . The p z electrons are coupled via Slater-Koster hopping parameters t ij = t(r i r j ) t(d) = t k (d) + t ? (d) Due to the internal twist between adjacent graphene sheets, a sufficient description of the interlayer hopping must include contributions from pp⇡ bonds pp⇡ = 2.8 eV as well as from pp bonds pp⇡ = 0.48 eV [7]. To this end, the factor n = d·êz |d| captures the out-of plane component of the electron transfer integral. Furthermore, e z is a unit vector which points perpendicular to the graphene sheets, c = 3.364 Å is the interlayer spacing of graphite, a = 1.42 Å is the distance between neighboring carbon atoms and 1 = 3.15 and 2 = 7.462 describe the exponential cutoff of the electron hopping. For the calculations of 3D twisted WSe 2 and boron nitride, we perform first principles calculations based on DFT as implemented in the Vienna Ab initio Simulation Package (VASP) [8] following similar methods used in previous works [9,10]. Plane-wave basis sets are employed with an energy cutoff of 550 eV for WSe 2 and 400 eV for boron nitride. The projector augmented wave method (PAW) [11] is used to construct the pseudopotentials felt by the valence electrons. For the calculations of 3D twisted WSe 2 , the exchange-correlation functionals are treated within the generalized gradient approximation (GGA) [12]. All the atoms are relaxed until the force on each atom is less than 0.01 eV/Å. Van der Waals interactions are included using the method of Tkatchenko and Scheffler [13] during the relaxation.
For the calculations of 3D twisted boron nitride, the exchange-correlation functionals are treated within the local density approximation (LDA). As shown in the previous work [10], the flat bands near the top of the valence band of twisted boron nitride do not change much upon relaxation. Therefore, as the calculations for 3D twisted boron with twist angles down to 2.28 degree are very heavy, we perform these large scale calculations for 3D twisted boron nitride without relaxation.   The effective structure defined by the charge accumulation points resembles AA-stacked graphene multilayers, where one of the two inequivalent sites, i.e. Q 1 , is shifted by D/2 in z-direction. Hence, in each of the two "effective" planes with z-coordinate 0 and D/2 , the charge puddles form a triangular lattice with lattice constant L.
The simplest SU(2) symmetric TB model that can be constructed for this configuration is a single-orbital two-band model that takes up to next-nearest neighbor intra-and interlayer hopping terms between the charge puddles into account Here, t 1 denotes the hopping amplitude between neighboring Q 1 -and Q 2 -sites , whereas t 2 and t 3 denote hopping processes between two Q 1 (Q 2 ) sites in either the same or different layers. The hopping parameters are determined by fitting the energy eigenvalues of H 0 to the flat bands of the ab-initio band structure of thBN. The single-particle spectrum for the periodic system is then modeled by the following Bloch Hamiltonian which is labeled in the order of the two charge localization points Q 1 ,Q 2 . The matrix elements are obtained by a Fourier transform of the real-space hopping matrix Eq. (3) to The matrix h k can then be diagonalized in orbital space for each momentum k to obtain the bandstructure ✏ b (k) and orbital-to-band transformation u b r (k), b = 1..N :  Fig. 3 in the main text. The structure constants D and L (see Fig. 3(a)) describe the spatial extent of the moiré cell in in-plane and out-of-plane direction, respectively.

S III. FLUCTUATION EXCHANGE APPROXIMATION IN MULTI-ORBITAL SYSTEMS
A. 3D multi-orbital susceptibility We define the free Matsubara Green's function in orbital-momentum (frequency) space where u b r are the orbital-to-band transformations that render the unperturbed Hamiltonian H 0 and the free Green's function g b (i!, k) = (i! e b (k)) 1 diagonal. The orbital indices r = {Q 1 , Q 2 } are restricted to the same unit cell and the momenta k lie in the first Brillouin zone. To this end, we define the free polarization functionˆ 0 (q) = 0 r,r 0 (q) as The Matsubara summation occuring in Eq. (8) can be evaluated analytically giving the well-known Lindhard function for multi-orbital systems where n F (✏) = (1 + e ✏ ) 1 is the Fermi function.

B. Random-phase approximation for multi-orbital systems
To study correlated states of matter in thBN that arise due to the presence of electronelectron interaction, we employ a repulsive Hubbard term for electrons with opposite spin Here, the occupation number operator is defined as n R,r i , = c † R,r i , c R,r i , . We calculate the renormalized interactions within the random-phase approximation (RPA) to analyze the electronic instabilities mediated by spin-fluctuation exchange between electrons to high order in the bare coupling U . Admittedly, this approach is biased as it does not capture the interwind fluctuations in different two-particle scattering channels, which would require the use of unbiased fRG techniques.

C. Pairing Symmetry
We may write the general particle-particle scattering vertex between electrons with op- For interaction values U < U crit the magnetic instabilities prescribed by the RPA analysis might not be strong enough to actually occur. In this paramagnetic regime, spin and charge fluctuations contained in the transverse and longitudinal spin channel can give rise to an effective interaction between electrons that may lead to the formation of Cooper pairs. The diagrams can be separated into spin-singlet and spin-triplet contributions, depending on whether pairing same/opposite spins, i.e. s 6 = s 0 (singlet) or s = s 0 (triplet). In general, we may separate the dependence of the gap parameter on momentum, spatial and spin degrees of freedom r 1 r 2 Since for spin singlet gaps the spin function (s 1 , s 2 ) is antisymmetric under exchange of indices, i.e. (s 1 , s 2 ) = (s 2 , s 1 ), the spatial and momentum dependence must be symmetric in order to fulfill the Pauli principle. For spin triplet gaps we hence require f (k, r 1 , r 2 ) = f ( k, r 2 , r 1 ). Since the system is assumed to be paramagnetic, pairing same/opposite spins yields the same result after explicitly symmetrizing/anti-symmetrizing the interaction vertex in orbital-momentum space.
Restricting the pairing to Kramer's degenerate pairs (k 1 , ") and ( k 1 , #), the particleparticle scattering vertex in FLEX approximation is given by transverse (t) and longitudinal (l) spin fluctuations. For simplicity, we will use the abbreviation r 1 r 2 !r 3 r 4 k 1 , k 1 !k 2 , k 2 = r 1 ,r 2 k 1 ,k 2 in the following. The diagrams contributing to these spin channels are shown below.
The effective spin-mediated interaction in the opposite spin channel thus becomes The spin-dependence of the susceptibilities occuring in the diagrammatic expansion above can be neglected due to the emergent SU(2) symmetry in the paramagnetic phase. To obtain the effective interaction in the singlet (s) and triplet (t) channel, we symmetrize/antisymmetrize the interaction vertex, i.e. s/t = 1 2 r 2 , k 1 , s 0 r 4 , k 2 , s 0 r 1 , k 1 , s r 3 , k 2 , s + r 2 , k 1 , s 0 r 3 , k 2 , s

D. Linearized Gap Equation
Assuming that spin-and charge fluctuation provide the superconducting glue in the system, we confine our considerations to the vicinity of the Fermi surface and only treat . The momentum transfer occurring in the polarization function in RPA is given by q t = k 1 +k 2 due to momentum conservation.
... . The momentum transfer occurring in the polarization function in RPA is given by q l = k 1 k 2 due to momentum conservation. Only an even number of particle-hole bubbles is allowed in the diagrammatic expansion in order to preserve the spin in the upper and lower leg of the pairing interaction. The diagrams that are resummed in the longitudinal channel are connected to the particle-hole susceptibility describing screening effects of the bare Coulomb interaction.
scattering processes of a Cooper pair from state (k, k) on fermi surface C b to the state (k 0 , k 0 ) on fermi surface C b 0 . To this end, we project the pairing vertex Eq. (14) from orbital to band space and only take intra-band scattering into account bb 0 The momenta k and k 0 are restricted to the various fermi surface sheets {C}, such that Here, v b F (k) = |r✏ b (k)| is the Fermi velocity at k 0 in band b. The largest eigenvalue > 0 for a given interaction kernel bb 0 s/t (k, k 0 ), will lead to the highest transition temperature T c and the corresponding eigenfunction b (k) determines the symmetry of the gap. The effective lattice model obtained from the charge accumulation points has point group D 3h .
The symmetry of the gap can thus by classified according to the irreducible representations of D 3h that are listed in Table 2.
The linearized gap equation (17) only accounts for the leading pairing symmetry at the transition temperature T c of the superconducting phase. In the case of degenerate eigenvalues (e.g. d-wave instabilities {d xz , d yz }) belonging to a two-dimensional irreducible representation, an arbitrary linear combination might be favored for T < T c . In order to find the linear combination that is preferred by the system below the transition temperature, we compute the free energy of the system Here, E b (k) is the energy of the Bogoliubov quasi-particles resulting from diagonalization of the BdG Hamiltonian where the form factors are are given by d xz (k) = sin(k x )sin(k z ) and d yz (k) = sin(k y )sin(k z ).
The free parameters ✓ and are extracted by minizing the free energy of the system Eq. (18). In Fig. 3 we show that the linear combination b k / [d xz (k) ± id yz (k)] = [sin(k x )sin(k z ) ± isin(k y )sin(k z )] is generally preferred for the given filling.
singlet triplet f y(y 2 3x 2 ) · p z f y(y 2 3x 2 ) Table 2. Pairing symmetries for the effective lattice model of thBN separated into contributions to spin singlet and triplet channel.