Engineering Three Dimensional Moiré Flat Bands
Lede Xian, Ammon Fischer, Martin Claassen, Jin Zhang, Angel Rubio, Dante M. Kennes
EEngineering Three Dimensional Moiré Flat Bands
Lede Xian,
1, 2, ∗ Ammon Fischer, ∗ Martin Claassen, JinZhang, Angel Rubio,
2, 5, 6, † and Dante M. Kennes
3, 2, ‡ Songshan Lake Materials Laboratory,523808 Dongguan, Guangdong, China Max Planck Institute for the Structure and Dynamics of Matter,Center for Free Electron Laser Science, 22761 Hamburg, Germany Institut für Theorie der Statistischen Physik,RWTH Aachen University and JARA-Fundamentals ofFuture Information Technology, 52056 Aachen, Germany Department of Physics and Astronomy,University of Pennsylvania, Philadelphia, PA 19104, USA Center for Computational Quantum Physics,Simons Foundation Flatiron Institute, New York, NY 10010 USA Nano-Bio Spectroscopy Group, Departamento de Fisica de Materiales,Universidad del País Vasco, UPV/EHU- 20018 San Sebastián, Spain a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec bstract We demonstrate that the concept of moiré flat bands can be generalized to achieveelectronic band engineering in all three spatial dimensions. For many two dimensionalvan der Waals materials, twisting two adjacent layers with respect to each other leadsto flat electronic bands in the two corresponding spatial directions – a notion sometimesreferred to as twistronics as it enables a wealth of physical phenomena. Within thistwo dimensional plane, large moiré patterns of nanometer size form. The basic conceptwe propose here is to stack multiple twisted layers on top of each other in a predefinedpattern. If the pattern is chosen such that with respect to the stacking direction oflayers, the large spatial moiré features are spatially shifted from one twisted layer tothe next, the system exhibits twist angle controlled flat bands in all of the three spatialdirections. With this, our proposal extends the use of twistronic to three dimensions.We exemplify the general concept by considering graphitic systems, boron nitride andWSe as candidate materials, but the approach is applicable to any two-dimensionalvan der Waals material. For hexagonal boron nitride we develope an ab initio fittedtight binding model that captures the corresponding three dimensional low-energyelectronic structure. We outline that interesting three dimensional correlated phases ofmatter can be induced and controlled following this route, including quantum magnetsand unconventional superconducting states. ∗ These authors contributed equally † Corresponding author: [email protected] ‡ Corresponding author: [email protected]
2n the past few years twisting adjacent layers of van der Waals materials has emerged as aversatile route to control the ratio between kinetic, potential and vibrational energy of two-dimensional systems. Central to the idea of twistronics, the selective suppression of kineticenergy scales permits tuning materials into a regime dominated by electronic interactions,as well as precise control over electronic filling via gating [1, 2]. Early experimental and the-oretical studies concentrated on graphetic systems of different kinds, such as twisted bilayergraphene [3–7], twisted double bilayer graphene [8–12], trilayer rhombohedral graphene onhexagonal boron nitride [13, 14] and twisted mono-bilayer graphene [15, 16]. More recently,twisted transition metal dichalcogenides (TMD) moved into the center of attention as an-other important class of van der Waals materials, such as WSe [17–19], MoS [20–22], andTMD heterostructures [23–26], allowing access to new regimes, beyond graphitic systems.With more of these phenomena within experimental reach twisted van der Waals materialsare increasingly viewed as potential avenues towards solid-state based platforms of quantumcontrol and quantum materials with properties on demand [2, 27]. Furthermore, the twistangle allows to control those systems to such a high degree that moiré aided metrologyoutperforms the current gold standard regarding structural questions in van der Waals ma-terials [28]. One important question guiding theoretical and experimental research effortsconcerns the exploration of the tremendous combinatoric space of chemical compositions ofvan der Waals materials to shed light onto the basic question of which additional phenom-ena might be accessible using twistronics. In the quickly expanding cosmos of twisted vander Waals materials, one guiding principle remains the control of the low-energy degrees offreedom, which might realize prototypical models of condensed matter research in a morerigid, clean and penetrable context [2]. In this spirit, and in addition to the directions al-ready outlined above, twisted hexagonal boron nitride was shown to harbor entire familiesof flat bands [29, 30] and twisted two-dimensional magnets, such as CrI , were identified torealize moiré skyrmion lattices and noncollinear twisted magnetic phases [31, 32]. Beyondrealizing different phenomena in two-dimensions using twist as a control paradigm one mightask whether entirely different dimensionalities can (effectively) be addressed. Twisted twodimensional monolayers of monochalcogenides (e.g. GeSe) were shown to allow access to theone-dimensional limit with twist [33, 34] providing the same unprecedented level of controlas in the two-dimensional counterparts.However, a practical generalization of twistronics to three dimensions with suppressed3inetic energy scales remains an outstanding challenge (aside from the possibility of addingfurther synthetic dimensions [35, 36]). Here, we address this by showing that stackingvan der Waals materials in a predefined fashion allows to engineer three dimensional moiréflat bands. With this advance the realization of three dimensional systems controlled bytwistronics is no longer elusive, completing the list of systems controlled by twistronics inone, two, and now three dimensions. Our idea works generically and relies mainly on basicgeometric arguments. Importantly, it is achievable within recent experimental advancesto fabricate bulk-like artificial twisted materials [37–39] and can be applied to any of themany van der Waals materials, which we will exemplify here for three important materials:graphitic systems, WSe and hexagonal boron nitride. The last of these will be examinedin more details and we provide a full ab initio characterization for its three dimensionaltwist-dependent band structures. We then consider the effects of correlations on the low-energy bands and find three dimensional magnetic and superconducting states which can berealized as a function of twist angle.The general idea is shown in Fig. 1 (a) relying on successively stacking up lattice sitesdefined by the moiré potential. In panel (b) we define the Brillouin zone, the in-plane andout-of-plane directions k ⊥ and k k as well as a cartoon of the corresponding band widths usedin the following discussion. We are looking for ideal 3D flat bands that satisfy the followingconditions: (1) the band width is controllable by twist angles at all three dimensions; (2)the system is periodic along all three dimensions such that it remains a well-defined crystal.Therefore, our work is distinct from previous works that study intrinsic 3D flat bands in somesolids [40–42] with a fixed band width and little tunability. Our approach also is distinctlydifferent to proposals such as the 3D chiral twisted structure [43, 44] that render the system aquasicrystal. To this end we start by considering three different stacking patterns illustratedin Fig. 1 (c-e).Firstly, we consider stacking monolayers of van der Waals materials in an alternatingfashion as depicted in panel (c) of Fig. 1, meaning that every second layer is aligned per-fectly while adjacent layers have a relative twist angle between them. This can be regardedas stacking twisted bilayers repeatedly. Following this route moiré patterns form by the in-terference between adjacent layers. Viewed top-down, the moiré lattice sites where in-planecharge density localizes, align on top of each other. As a consequence, the electronic bandsbecome flat within the plane; just as is the case for two twisted sheets of van der Waals ma-4 igure Stacking ideas to construct 3D flat bands.
Panel (a) exemplifies the general ideaof stacking up lattice sites defined by the moiré potential in a successive way to approach the threedimensional limit. Panel (b) shows the Brillouin zone and defines the out-of plane k ⊥ and in-plane k k directions (left) as well as a cartoon of the corresponding band widths w ⊥ and w k along thesedirections (right). Panels (c)-(e) show three different configurations to create three dimensionalmaterials out of stacks of two dimensional twisted van der Waals materials. The orange areaindicate the moiré lattice sites where the charge density localized. When these moiré lattice sitesalign atop of each other (c) the bands become flat in two of the three directions only. When usingspacer layers (d) the band width can be reduced in all three dimensions, but the band width alongthe stacking direction cannot be controlled by the twist angle. If moiré lattice sites by geometricreasons (see main text) are shifted from layer to layer (c), three dimensional flat bands with theflatness in all three dimensions being controlled by the twist angle emerge. terials. Conversely, the alignment of moiré regions along the out-of-plane direction retainssubstantial band dispersion in this direction due to significant amount of hybridization be-tween moiré sites at adjacent layers (see Fig. S1 in the SI for an example of 3D twisted boronnitride with such stacking). While this allows to effectively engineer quasi-one-dimensionalsystems with very low residual coupling along the remaining two spatial directions – aninteresting opportunity of materials engineering in its own rights (for a similar quasi-one-dimensional system it was shown that in-plane confinement albeit imposed by an magneticfield gives rise to a 3D quantum Hall effect [45])– it does not allow to realize three dimensionalmoiré flat bands. 5econdly, one might consider the case in which twisted van der Waals materials with flatbands in their two spatial directions are stacked on top of each other with an insulatingbuffer layer in between as shown in panel (d) of Fig. 1. The properties and thickness ofthe buffer-layer could then be adjusted such that the hopping from one twisted sheet ofvan der Waals materials to the next sheet is suppressed substantially. This would leadto flat electronic bands in all three spatial dimensions. However, such an approach hasmultiple problems. First, the flatness of the bands in the out-of-plane direction is mainlydetermined by the residual coupling between neighboring moiré charge localization sitesacross the insulating layers. This limits available band structures that can be engineeredquite substantially compared to the flexible control that twist angle offers with respect tothe remaining two directions. Second, with a buffer layer, there is no guarantee to keepthe moiré sites across the buffer layer well aligned, i.e., the centers of the moiré sites ofneighboring twisted pairs can relocate to different in-plane positions when stacking up. Thiscould introduce significant amount of disorder along the out-of-plane direction such that thesystem is no longer a well-defined crystal.To remedy the shortcomings of the previous two stacking approaches, thirdly, we presentan idea using a stacking sequence in which the moiré charge localization sites simply dueto geometric considerations do not form atop of each other, but are shifted with respectto the out-of-plane directions of the van der Waals materials used. Such a configurationcan be constructed by expanding the basic stacking unit, e.g., from a twisted bilayer toa twisted double bilayer. This is visualized in panel (e) of Fig. 1. Compared to the firstapproach, the repeating unit along the out-of-plane direction is changed from layers 1,2 inFig. 1(a) to layers 1-4 in Fig. 1(c). In this approach, layers 2,3 and layers 4,1 remain attheir intrinsic Bernal AB stacking or AA’ stacking sequence as in the pristine bulk materialand the twisting takes place only between layers 1 and 2 as well as 3 and 4 in the notationof Fig. 1(c). Although the in-plane crystal axis of layer 2 is aligned with those of layer3, the atomic positions of the two layers are translated with respect to each other (as inintrinsic AB stacking), or flipped (as in intrinsic AA’ stacking). The same happens for layer4 and 1. This naturally displaces the moiré charge localization sites in these layers withrespect to each other, which are now separated by the moiré length scale. As the twist angleis decreased and the in-plane distance between moiré sites increases, so does the distancebetween sites on adjacent bilayers. Therefore, the idea of using natural (or intrinsic) bilayer6s a stacking unit to construct alternating patterns, allows to engineer robust flat bandsin all three dimensions, with the flatness being continuously controlled by the twist. Thus,it satisfy the condition (1) we set for a ideal flat band system. Moreover, as we will showbelow, such an approach will also generate local stacking regions that resemble the stackingsequence in the pristine bulk crystal, which can act as a low-energy stabilization centerto prevent disorder along the out-of-plane direction. Therefore, this approach also meetscondition (2) of a nearly ideal robust flat band crystallographic system.We put this very general idea to the test by first performing ab initio and tight bindingbased characterizations of such stacked materials using bilayers of graphene, WSe as wellas hexagonal boron nitride. All of these materials were successfully studied in the past forthe twisted single bilayer case rendering them ideal starting points to explore the idea weput forward here. The results are summarized for a twist angle of . ◦ for graphene and . ◦ for both WSe and boron nitride in Fig. 2. We show side and top views of the realspace stacking in panels (a)-(c) for graphene, WSe and boron nitride, respectively. Theunit cell for these bulk twisted systems is formed by a twisted double bilayer as highlightedby dashed lines in the top row and the solid brackets in the third row. The bottom panelsof (a)-(c) show the local stacking sequence in the three representative regions shown in themiddle row panels. The stacking sequence in region I is exactly the same as that in theintrinsic untwisted bulk crystal (AB stacking as in graphite, AA’ stacking as in 2H WSe and bulk boron nitride). This region is expected to pin the in-plane alignment of the layersand naturally prevent accidental layer displacements similar to what is discussed for the caseof twisted trilayer graphene [46]. In panels (d)-(f) we show the respective band structures.Since the angle is not very small for twisted WSe and twisted boron nitride, the bands stillshow substantial dispersion and the set of bands that is flattening, marked in red in Fig. 2,have not fully separated from other bands yet (a more detailed study of decreasing the twistangle further is given below). The important feature to note is the width of the red set ofbands, with respect to variations of both the out-of-plane k ⊥ ( A → Γ path in the Brillouinzone) and in-plane k k ( Γ → K → M → Γ path in the Brillouin zone) momentum. Both ofthese band widths in- and out-of-plane, called W k and W ⊥ respectively, decrease as the twistangle approaches smaller values (see below and SI). Eventually these bands separate fromthe rest giving rise to perfectly isolated flat bands, with comparable kinetic energy scales inall three spatial dimensions, whose magnitude can be tuned by the twist, as demonstrated7 igure
3D Flat bands for different materials, Graphene, WSe and boron nitride .Atomic structures of 3D twisted graphene (a), WSe (b), boron nitride (c). The top and the middlepanels show the perspective and the top views of the structures, respectively. The unit cells areindicated with dashed lines. The bottom panels show the local stacking sequence in the region I,II and III indicated in the middle panels. The repeating units along the out-of-plane direction areindicated with solid brackets. (d-f) The corresponding band structures for graphene at 1.3 degrees(d), for WSe at 5.08 degrees (e), and for boron nitride at 5.08 degrees (f). For smaller angles thebands become increasingly flat and detach from other bands. next.We stress that the idea we present here is general and allows three dimensional flat bandengineering also in other materials even beyond the ones discussed explicitly above. Howeverwe are going to illustrate the relevance of this new concept to electronic band engineeringtaking the specific case of hBN (see Fig. S2 in the SI for the case of WSe2). The choice of8 igure Low-energy model of three dimensional twisted boron nitride (thBN) . (a)Stacking pattern and moiré unit cell of thBN for θ = 5 . ◦ . Emerging charge accumulation points(olive green) form an effective lattice that resembles AA-stacked hexagonal multilayers, where onesite is shifted by D/ in z -direction (blue and gold spheres). (b) Results from our ab initio fittedTight-binding (TB) approach (see methods section for details) for different twist angles θ , taking upto next-nearest neighbor inter- and intralayer hopping terms t ,.., t into account. The band widthdecreases continuously with the twist angle and, eventually, the low-energy bands detach from therest of the spectrum. boron nitride is made for convenience (and for being widely used as protective 2D material),as the absence of sharp magic angles, makes it particularly feasible to large scale numericsproviding a full-fledged ab-initio characterization of the material’s band structure in threedimensions. As the relaxation of twisted boron nitride does not significantly alter the bandstructure according to the previous work [29], we fix the atomic structure in the large-scaleab-initio calculations. Our results are summarized in Fig. 3. In panel (a) we show top andside views of the charge localization within the moiré unit cell and the position of the Band N atoms. We choose intrinsic AA’ bilayers (as in pristine bulk hBN) as the buildingblocks of our three dimensional structure and therefore there are no ferroelectric domainsas recently reported for twisted bilayer systems [47, 48]. Generalizing this constitutes anintriguing avenue of future research. In panel (b) we summarize the ab-initio band structure9btained for different twist angles. Importantly, as we approach smaller values of the twistangle, both the in-plane and out-of-plane band width W k and W ⊥ , which are the same withsuch stacking, decrease and the flat bands detach from the other bands in the spectrum.Strong charge localization marked by blue and golden spheres in panel (a) highlight theemergence of a corresponding three dimensional effective low-energy tight binding model(see methods for details). The model we obtain for the case of boron nitride is an in-planetriangular lattice stacked out of plane, such that the lattice sites of one plane reside in thecenter of the triangular lattice of the next plane.The success of fitting the flat bands within such a low-energy tight binding model in-cluding only short ranged hoppings t , t and t on the moiré scale as denoted in panel (a)is demonstrated in panel (b). By fitting the three hopping parameters to the full ab-initioband structure for different angles, almost perfect agreement is achieved (consistent with theearlier study of a single twisted hBN bilayer [29]). The smaller panels left to panel (b) showthe extracted values of the hopping as well as the overall bandwidth (in this case W ⊥ = W k ),demonstrating the success of three dimensional flat band engineering by the twist proposedhere. Such an effective low-energy model is immensely useful as it can be treated much moreefficiently.We employ the effective low-energy tight binding model to outline a putative unexploredphase diagram that could be accessible via the three dimensional twistronics approach. Tothis end, we consider a local Hubbard interaction U added to the effective flat band modelfor θ = 3 . ◦ as discussed above. We note that a more realistic model should include longerranged interactions as well, which should be characterized from first principles. Such a studyis unfortunately beyond the scope of the present work and most likely requires a fundamen-tal methodological advance to treat the huge three dimensional unit cell (containing manytens of thousands of atoms at small twist angles). Here, however, we provide the first stepalong a characterization of elusive and exciting correlation effects and aim to identify inter-esting states of matter already present at the level of a Hubbard interaction. To achievethis we first perform a random phase approximation (RPA) study of the system (see sup-plemental material) and identify a plethora of magnetic instabilities. A putative magneticphase diagram is summarized in panel (a) of Fig. 4. As expected we find ferromagnetic(FM) ordering tendencies as the flat bands are either filled or empty, albeit with a ratherlarge critical U crit driving the transition. Due to the bipartiteness of the lattice (sublattice10 igure Correlated phases in three dimensional twisted hexagonal boron nitride fortwist angle θ = 3 . ◦ . (a-c) The RPA analysis reveals a variety of magnetic states including anti-ferromagnetic order (AFM) at charge neutrality, ferromagnetic order (FM) for strong electron/holedoping as well as general spin-density waves q -SDW with periodic patterns in all three spatialdimensions. (d) For µ = µ + 5 meV the Fermi surface of thBN is almost perfectly nested along thevector Q α resulting in a strong enhancement of particle-particle scattering between these sheets.In particular, the preferred superconducting gap symmetry (e) is two-fold degenerate and of type ( d xz , d yz ) . Below T c , the system will minimize its Ginzburg-Landau free energy by assuming thechiral linear combination ( d xz ± id yz ) and thus the gap parameter breaks time-reversal symmetry. A and B being charge localization sites marked by blue and golden spheres in Fig. 3(a)) wefind antiferromagnetic (AFM) ordering at half filling. In between these two phases a moregeneral spin density wave with filling dependent ordering vector q emerges. q is illustratedin panel (b) of Fig. 4, while panel (c) illustrates the magnetic ordering in real space for fourdifferent examples depending on the filling: FM, AFM or a spin density wave with a wavevector lying either in the in-plane ( q k ) or out-of-plane ( q ⊥ ) direction.For interaction values U < U crit there is no magnetic ordering and the system is param-agnetic. In this regime, spin and charge fluctuations may provide an effective pairing gluebetween the electrons leading to the formation of Cooper pairs. To pin down the pairinginstability mediated by spin and charge fluctuations, we take the RPA corrected interactionvertex in the fluctuation exchange approximation (see supplemental material) and linearize11he superconducting gap equation around the critical temperature T c for slight electron dop-ing µ = µ + 5 meV around half filling of the three dimensional flat bands µ = 2 . meV.In this scenario, only scattering events between Cooper pairs in the vicinity of the Fermisurface sheets C α,β , Fig. 4 (d), contribute notably to the formation of a superconductingstate with order parameter ∆ µ k . The fact that (i) the Fermi surface sheets C α are (almost)perfectly connected by the nesting vector Q α at which particle-particle scattering is strongestand (ii) the effective pairing glue contained in the spin-singlet channel is purely repulsivefor all scattering events, electron-electron pairing is conditioned on a relative sign changebetween the pairing form factors µ connected by the vector Q α i.e. ∆ µ k = − ∆ µ k + Q α . Thelinearized gap equation may be written as an eigenvalue problem where the eigenfunctions ∆ µ k corresponding to the largest eigenvalue λ yield the symmetry of the most prominentsuperconducting state, Fig. 4 (e). Our calculations reveal that the leading gap parameter isof spin-singlet type and is two-fold degenerate with symmetry classification ( d xz , d yz ) . Thetwo d -wave solutions are characterized by a nodal line along the k x - and k y -direction of theBrillouin zone and thus the system will minimize its Ginzburg-Landau free energy F below T c by assuming the chiral linear combination ( d xz ± id yz ) (see supplemental material) whichbreaks time-reversal symmetry.In conclusion our work generalizes the idea of two dimensional twistronics to the threedimensional realm. The main notion relies on cleverly stacking adjacent layers in such a waythat the hopping between adjacent charge puddles in all three dimensions gets successivelysuppressed as the twist angle is lowered. We argue that the proposed stacking methodis robust towards small twist angle imperfections and inhomogeneities which might varywithin one plane or between adjacent layers. Even more so, since our construction reliespurely on geometric arguments even in the presence of imperfections these would simplyreflect in slightly inhomogeneous hoppings, like they are present in any (not perfectly clean)crystal and which do not change the overall physics significantly. On the contrary, controlledvariations of the twist angle might allow to control the effective disorder and therefore providea long sought after inroad into tunable strongly disordered systems, from a condensed matterperspective. With this we add the three dimensional realm to the list of low-energy modelsthat can effectively be realized in moiré structures. As a side product we also provide analternative of engineering effectively quasi-one-dimensional structures, by placing the moirésites on top of each other (first scenario in Fig. 1). This is not at the center of attention in12his work but allows to access similar physics as discussed in the context of the quantumHall effect in Ref. [45]. We already reported on the rich behavior of correlation drivenphases in engineered three dimensional flat bands above, but another intriguing avenue offuture research should also address the question of three dimensional flat band engineeringfor purposes of controlling topological properties. This provides a very rich playground thatdirectly opens up by our approach and the future will tell which topological phenomena,such as Weyl physics, and correlated phases beyond the ones discussed here might be tunableby three dimensional twistronics. METHODS
The low-energy physics of twisted hBN (thBN) is captured by an effective three-dimensional tight-binding (TB) model that includes hopping terms between emerging chargelocalization points at Q = ( , , ) and Q = ( , , in the moiré unit cell. The coordi-nates are given with respect to the in-plane ( k ) and out-of plane ( ⊥ ) Bravais lattice vectors L k = ( L, , , L k = R ( π/ L and L ⊥ = (0 , , D ) . The lattice constants D and L aretwist-angle dependent and describe the spatial extent of the moiré pattern, see table 1.The effective structure defined by the charge accumulation points resembles AA-stackedgraphene multilayers, where one of the two inequivalent sites, i.e. Q , is shifted by D/ inz-direction. Hence, in each of the two "effective" planes with z -coordinate and D/ , thecharge puddles form a triangular lattice with lattice constant L .The simplest SU (2) symmetric TB model that can be constructed for this configuration istwist angle θ Hopping parameters(meV) Lattice constants(Å) t t t L D5.08 ◦ ◦ ◦ ◦ Table 1. Hopping parameters of the effective SU (2) -symmetric tight-binding model for differenttwist angles θ according to 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.
13 single-orbital two-band model that takes up to next-nearest neighbor intra- and interlayerhopping terms between the charge puddles into account H = t X h i,j i c † i c j + t X h i,j i ⊥ c † i c j + t X h i,j i k c † i c j . (1)Here, t denotes the hopping amplitude between neighboring Q - and Q -sites , whereas t and t denote hopping processes between two Q ( Q ) sites in either the same or differentlayers. The hopping parameters are determined by fitting the energy eigenvalues of H tothe flat bands of the ab-initio band structure of thBN. The single-particle spectrum for theperiodic system is then modeled by the following Bloch Hamiltonian H = X k h k = X k h ( k ) h ( k ) h ∗ ( k ) h ( k ) , (2)which is labeled in the order of the two charge localization points Q , Q . The matrixelements are obtained by a Fourier transform of the real-space hopping matrix Eq. (1) to(Bloch) momentum space h ( k ) = 2 t [ cos ( k · ( L − L )) + cos ( k · L ) + cos ( k · L )] + 2 t cos ( k · L ) ,h ( k ) = t (cid:2) e − i k · L + e − i k · L (cid:3) (cid:2) e − i k · L (cid:3) . (3)The matrix h k can then be diagonalized in orbital space for each momentum k to obtainthe bandstructure (cid:15) b ( k ) and orbital-to-band transformation u b r ( k ) , b = 1 ..N : H = X k ,b (cid:15) b ( k ) γ † k ,b γ k ,b with γ k ,b = u b r ( k ) c k , r . (4) [1] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nat. Phys. (2020), 10.1038/s41567-020-0906-9.[2] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N.Basov, A. Pasupathy, and A. Rubio, Nature Physics (2020).
3] R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U.S.A. , 12233 (2011).[4] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero,Nature , 43 (2018).[5] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi,K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature ,80 (2018).[6] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F.Young, and C. R. Dean, Science , 1059 (2019).[7] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe,T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Nature , 95 (2019).[8] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. H. Najafabadi, K. Watanabe,T. Taniguchi, A. Vishwanath, et al., Nature , 221 (2020).[9] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe,et al., Nat. Phys. , 1 (2020).[10] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, andP. Jarillo-Herrero, Nature , 1 (2020).[11] G. W. Burg, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc, Phys. Rev.Lett. , 197702 (2019).[12] C. Rubio-Verdú, S. Turkel, L. Song, L. Klebl, R. Samajdar, M. S. Scheurer, J. W. F. Venderbos,K. Watanabe, T. Taniguchi, H. Ochoa, L. Xian, D. M. Kennes, R. M. Fernandes, Ángel Rubio,and A. N. Pasupathy, arXiv preprint arXiv:2009.11645 (2020).[13] G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watan-abe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature , 215 (2019).[14] G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi,J. Jung, Y. Zhang, and F. Wang, Nat. Phys. , 237 (2019).[15] S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden,X. Xu, C. R. Dean, et al., Nature Physics (2020).[16] Y. Shi, S. Xu, M. M. A. Ezzi, N. Balakrishnan, A. Garcia-Ruiz, B. Tsim, C. Mullan, J. Barrier,N. Xin, B. A. Piot, et al., arXiv preprint arXiv:2004.12414 (2020).[17] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, . Bai, B. Kim, et al., Nat. Mater. , 1 (2020).[18] L. An, X. Cai, D. Pei, M. Huang, Z. Wu, Z. Zhou, J. Lin, Z. Ying, Z. Ye, X. Feng, et al.,Nanoscale Horiz. , 1309 (2020).[19] H. Pan, F. Wu, and S. D. Sarma, arXiv preprint arXiv:2004.04168 (2020).[20] M. Liao, Z. Wei, L. Du, Q. Wang, J. Tang, H. Yu, F. Wu, J. Zhao, X. Xu, B. Han, et al.,Nature Communications , 1 (2020).[21] M. H. Naik and M. Jain, Phys. Rev. Lett. , 266401 (2018).[22] L. Xian, M. Claassen, D. Kiese, M. M. Scherer, S. Trebst, D. M. Kennes, and A. Rubio,arXiv:2004.02964 (2020).[23] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. , 026402 (2018).[24] Y. Zhang, N. F. Yuan, and L. Fu, arXiv:1910.14061 (2019).[25] E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang,K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie,A. Zettl, and F. Wang, Nature , 359 (2020).[26] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald,J. Shan, and K. F. Mak, Nature , 353 (2020).[27] D. N. Basov, R. D. Averitt, and D. Hsieh, Nat. Mater. , 1077 (2017).[28] D. Halbertal, N. R. Finney, S. S. Sunku, A. Kerelsky, C. Rubio-Verdú, S. Shabani, L. Xian,S. Carr, S. Chen, C. Zhang, L. Wang, D. Gonzalez-Acevedo, A. S. McLeod, D. Rhodes,K. Watanabe, T. Taniguchi, E. Kaxiras, C. R. Dean, J. C. Hone, A. N. Pasupathy, D. M.Kennes, A. Rubio, and D. N. Basov, arXiv preprint arXiv:2008.04835 (2020).[29] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Nano Lett. , 4934(2019).[30] X.-J. Zhao, Y. Yang, D.-B. Zhang, and S.-H. Wei, Physical Review Letters , 086401 (2020).[31] Q. Tong, F. Liu, J. Xiao, and W. Yao, Nano Lett. , 7194 (2018).[32] K. Hejazi, Z.-X. Luo, and L. Balents, Proc. Natl. Acad. Sci. U.S.A. , 10721 (2020).[33] D. M. Kennes, L. Xian, M. Claassen, and A. Rubio, Nat. Commun. , 1 (2020).[34] T. Kariyado and A. Vishwanath, Phys. Rev. Research , 033076 (2019).[35] M. Lohse, C. Schweizer, H. M. Price, O. Zilberberg, and I. Bloch, Nature , 55 (2018).[36] O. Zilberberg, S. Huang, J. Guglielmon, M. Wang, K. P. Chen, Y. E. Kraus, and M. C.Rechtsman, Nature , 59 (2018).
37] F. Liu, W. Wu, Y. Bai, S. H. Chae, Q. Li, J. Wang, J. Hone, and X.-Y. Zhu, Science ,903 (2020).[38] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe,T. Taniguchi, A. Vishwanath, and P. Kim, arXiv preprint arXiv:2012.02773 (2020).[39] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv preprintarXiv:2012.01434 (2020).[40] M. Kang, S. Fang, L. Ye, H. C. Po, J. Denlinger, C. Jozwiak, A. Bostwick, E. Rotenberg,E. Kaxiras, J. G. Checkelsky, et al., arXiv:2002.01452 (2020).[41] Z. Liu, M. Li, Q. Wang, G. Wang, C. Wen, K. Jiang, X. Lu, S. Yan, Y. Huang, D. Shen, et al.,Nature communications , 1 (2020).[42] Y. Zhou, K.-H. Jin, H. Huang, Z. Wang, and F. Liu, Physical Review B , 201105 (2019).[43] T. Cea, N. R. Walet, and F. Guinea, Nano Letters , 8683 (2019).[44] F. Wu, R.-X. Zhang, and S. D. Sarma, Phys. Rev. Res. , 022010 (2020).[45] F. Tang, Y. Ren, P. Wang, R. Zhong, J. Schneeloch, S. A. Yang, K. Yang, P. A. Lee, G. Gu,Z. Qiao, and L. Zhang, Nature , 537 (2019).[46] S. Carr, C. Li, Z. Zhu, E. Kaxiras, S. Sachdev, and A. Kruchkov, Nano Letters (2020).[47] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv preprintarXiv:2010.06600 (2020).[48] C. Woods, P. Ares, H. Nevison-Andrews, M. Holwill, R. Fabregas, F. Guinea, A. Geim,K. Novoselov, N. Walet, and L. Fumagalli, arXiv preprint arXiv:2010.06914 (2020). upplemental Material forEngineering Three Dimensional Moiré Flat Bands Lede Xian,
1, 2, ∗ Ammon Fischer, ∗ Martin Claassen, JinZhang, Angel Rubio,
2, 5, 6, † and Dante M. Kennes
3, 2, ‡ Songshan Lake Materials Laboratory,523808 Dongguan, Guangdong, China Max Planck Institute for the Structure and Dynamics of Matter,Center for Free Electron Laser Science, 22761 Hamburg, Germany Institut für Theorie der Statistischen Physik,RWTH Aachen University and JARA-Fundamentals ofFuture Information Technology, 52056 Aachen, Germany Department of Physics and Astronomy,University of Pennsylvania, Philadelphia, PA 19104, USA Center for Computational Quantum Physics,Simons Foundation Flatiron Institute, New York, NY 10010 USA Nano-Bio Spectroscopy Group, Departamento de Fisica de Materiales,Universidad del País Vasco, UPV/EHU- 20018 San Sebastián, Spain ∗ These authors contributed equally † Corresponding author: [email protected] ‡ Corresponding author: [email protected] a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec I. COMPUTATIONAL DETAILS FOR 3D TWISTED GRAPHENE, WSE2 ANDBORON NITRIDE
For the calculations of 3D twisted graphene, we construct the unit cell with a twisteddouble bilayer graphene at twist angles close to 0 degree, and impose periodic boundarycondition along all three dimensions. As it is not realistic to optimize such a large systemwith density functional theory (DFT) calculations, we fix the lattice constant along theout-of-plane direction to be 13.415 Å, and set the in-plane lattice constant according tothe twist angles such that it corresponds to 2.46 Å for a 1x1 cell. The atomic structureis relaxed using the LAMMPS code [1] with the same parameters as described in [2]. Theintralayer interactions within each graphene layer are modeled via the second-generationreactive empirical bondorder (REBO) potential [3]. The interlayer interactions are modeledvia the Kolmogorov-Crespi (KC) potential [4], using the recent parametrization of [5]. Therelaxation is performed using the fast inertial relaxation engine (FIRE) algorithm [6].We calculate the band structures for 3D twisted graphene with the following tight-binding Hamiltonian [7]: H= P i , j t ij | i i h j | , where t ij is the hopping parameter be-tween pz orbitals at the two lattice sites r i and r j and it has the following form: t ij = n γ exp h λ (cid:16) − | r i − r j | a (cid:17)i + (1 − n ) γ exp h λ (cid:16) − | r i − r j | c (cid:17)i , a=1.42 Å, c=3.364 Å, γ = − . eV , γ = 0 . eV , λ = 3 . and λ = 7 . .For the calculations of 3D twisted WSe and boron nitride, we perform first principlescalculations 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 setsare employed with an energy cutoff of 550 eV for WSe and 400 eV for boron nitride. Theprojector augmented wave method (PAW) [11] is used to construct the pseudopotentials feltby the valence electrons. For the calculations of 3D twisted WSe , the exchange-correlationfunctionals are treated within the generalized gradient approximation (GGA) [12]. All theatoms are relaxed until the force on each atom is less than 0.01 eV/Å. Van der Waals interac-tions 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 aretreated 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 muchupon relaxation. Therefore, as the calculations for 3D twisted boron with twist angles down2 igure
1. Band structures of 3D twisted boron nitride with type I stacking at 5.08 degrees (a) and3.89 degrees (b). The band width of the top valence band along the in-plane at k z =0 decreaseswith twist angles, while the band width along the out-of-plane direction remains highly dispersive. Figure
2. Band structures of 3D twisted WSe with type III stacking at 5.08 degrees (a) and 3.89degrees (b). The band width along both the in-plane and the out-of-plane directions decreases withtwist angles. to 2.28 degree are very heavy, we perform these large scale calculations for 3D twisted boronnitride without relaxation. 3 II. FLUCTUATION EXCHANGE APPROXIMATION IN MULTI-ORBITALSYSTEMSA. 3D multi-orbital susceptibility
We define the free Matsubara Green’s function in orbital-momentum (frequency) spaceas g r , r ( iω, k ) = ( iω − ( H ( k )) r , r ) − = X b u b r ( k ) g b ( iω, k ) u b ∗ r ( k ) = r r ( b, k ) (1)where u b r are the orbital-to-band transformations that render the unperturbed Hamiltonian H and the free Green’s function g b ( iω, k ) = ( iω − e b ( k )) − diagonal. The orbital indices r = { Q , Q } are restricted to the same unit cell and the momenta k lie in the first Brillouinzone. To this end, we define the free polarization function χ ( q ) as χ r , r ( q ) = χ r , r ( q , iω ) = 1 N k β X k ,ω g r , r ( iω , k ) g r , r ( i ( ω + ω ) , k + q ) . (2)The Matsubara summation occuring in Eq. (2) can be evaluated analytically giving thewell-known Lindhard function for multi-orbital systems χ r , r ( q , iω ) = 1 N X k ,b,b n F ( (cid:15) b ( k )) − n F ( (cid:15) b ( k + q )) iν + (cid:15) b ( k ) − (cid:15) b ( k + q ) u b r ( k ) u b ∗ r ( k ) u b ∗ r ( k + q ) u b r ( k + q ) , (3)where n F ( (cid:15) ) = (1 + e β(cid:15) ) − 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 electron-electron interaction, we employ a repulsive Hubbard term for electrons with opposite spin σ ∈ {− , } with σ = − σ residing on site r in moiré supercell R V = 12 X R , r i ,σ U n R , r i ,σ n R , r i ,σ = 12 N X k , k , q X r ,σ U c † r k ,σ c † r k ,σ c r k − q ,σ c r k + q ,σ (4)4ere, the occupation number operator is defined as n R , r i ,σ = c † R , r i ,σ c R , r i ,σ . We calculatethe renormalized interactions within the random-phase approximation (RPA) to analyzethe electronic instabilities mediated by spin-fluctuation exchange between electrons to highorder in the bare coupling U . Admittedly, this approach is biased as it does not capture theinterwind fluctuations in different two-particle scattering channels, which would require theuse of unbiased fRG techniques. r , k , ↓ r , k , ↓ r , k + q, ↑ r , k + q, ↑ = r , k , ↓ r , k , ↓ r , k + q, ↑ r , k + q, ↑ U + r k ↓ r k + q ↑ r r U (5)The renormalized interaction in RPA approximation Eq. (5) is then given by ˆ V RPA ( q ) = U/ [1 + U ˆ χ ( q )] . Magnetic instabilities can be classified according to a generalized Stonercriterion: The effective (RPA) interaction diverges, when the smallest eigenvalue λ of ˆ χ ( q , iω ) reaches − /U , marking the onset of magnetic order for all interaction strengths U ≥ U crit. = − /λ . The corresponding eigenvector v (0) ( q , iω ) is expected to dominate thespatial structure of orbital magnetisation. C. Pairing Symmetry
We may write the general particle-particle scattering vertex between electrons with op-posite momenta ( k , − k ) → ( k , − k ) as V = 12 N X s,s X r ,..., r X k , k Γ r r → r r k , − k → k , − k c † r k s c † r − k s c r − k s c r k s = r , − k , s r , − k , s r , k , s r , k , s Γ r r → r r k , − k → k , − k (6)For interaction values U < U crit the magnetic instabilities prescribed by the RPA analysismight not be strong enough to actually occur. In this paramagnetic regime, spin and chargefluctuations contained in the transverse and longitudinal spin channel can give rise to aneffective interaction between electrons that may lead to the formation of Cooper pairs. The5iagrams can be separated into spin-singlet and spin-triplet contributions, depending onwhether pairing same/opposite spins, i.e. s = s (singlet) or s = s (triplet). In general, wemay separate the dependence of the gap parameter on momentum, spatial and spin degreesof freedom ∆ r r k s s = f ( k , r , r ) χ ( s , s ) . (7)Since for spin singlet gaps the spin function χ ( s , s ) is antisymmetric under exchangeof indices, i.e. χ ( s , s ) = − χ ( s , s ) , the spatial and momentum dependence must besymmetric in order to fulfill the Pauli principle. For spin triplet gaps we hence require f ( k , r , r ) = − f ( − k , r , r ) . Since the system is assumed to be paramagnetic, pairingsame/opposite spins yields the same result after explicitly symmetrizing/anti-symmetrizingthe interaction vertex in orbital-momentum space.Restricting the pairing to Kramer’s degenerate pairs ( k , ↑ ) and ( − k , ↓ ) , the particle-particle scattering vertex in FLEX approximation is given by transverse ( t ) and longitudinal( l ) spin fluctuations. For simplicity, we will use the abbreviation Γ r r → r r k , − k → k , − k = Γ r , r k , k inthe following. The diagrams contributing to these spin channels are shown below.The effective spin-mediated interaction in the opposite spin channel thus becomes Γ r r → r r k , − k → k , − k = δ r , r δ r , r (cid:20) ˆ U + U ˆ χ ( q l )1 − U ˆ χ ( q l ) (cid:21) + δ r , r δ r , r (cid:20) − U ˆ χ ( q t )1 + U ˆ χ ( q t ) (cid:21) (8)The spin-dependence of the susceptibilities occuring in the diagrammatic expansion abovecan be neglected due to the emergent SU (2) symmetry in the paramagnetic phase. Toobtain the effective interaction in the singlet ( s ) and triplet ( t ) channel, we symmetrize/anti-symmetrize the interaction vertex, i.e. Γ s/t = 12 r , − k , s r , − k , s r , k , s r , k , s + σ r , − k , s r , k , s r , k , s r , − k , s (9) D. Linearized Gap Equation
Assuming that spin- and charge fluctuation provide the superconducting glue in thesystem, we confine our considerations to the vicinity of the Fermi surface and only treat6 , − k , ↓ r , − k , ↓ r , k + q, ↑ r , k + q, ↑ =Γ t r , r k , k r k r k + q t r r U U + k k k + q t k + q t U U U + ... (a) Diagrams contributing to the transverse spin-fluctuation mediated pairing interaction Γ t r , r k , k . Themomentum transfer occurring in the polarization function in RPA is given by q t = k + k due to momentumconservation. r , − k , ↓ r , − k , ↓ r , k + q, ↑ r , k + q, ↑ =Γ l r , r k , k r , − k , ↓ r , − k , ↓ r , k + q, ↑ r , k + q, ↑ U + r , − k , ↓ r , − k , ↓ r , k + q, ↑ r , k + q, ↑ k + q l ↓ k ↓ rr k + q l ↑ k ↑ r r + ... (b) Diagrams contributing to the longitudinal spin-fluctuation mediated pairing interaction Γ l r , r k , k . Themomentum transfer occurring in the polarization function in RPA is given by q l = k − k due to momentumconservation. Only an even number of particle-hole bubbles is allowed in the diagrammatic expansion in orderto preserve the spin in the upper and lower leg of the pairing interaction. The diagrams that are resummedin the longitudinal channel are connected to the particle-hole susceptibility describing screening effects of thebare Coulomb interaction. scattering processes of a Cooper pair from state ( k , − k ) on fermi surface C b to the state( k , − k ) on fermi surface C b . To this end, we project the pairing vertex Eq. (8) from orbitalto band space and only take intra-band scattering into account Γ bb s/t ( k , k ) = Re " X r , r , r , r Γ s/t u b ∗ r ( k ) u b ∗ r ( − k ) u b r ( k ) u b r ( − k ) . (10)The momenta k and k are restricted to the various fermi surface sheets { C } , such that k ∈ C b and k ∈ C b with b and b being the band indices of the fermi sheets. Neglecting thefrequency dependence of Γ , we can proceed further by considering only the real part of thepairing interaction. We then solve the linearized gap equation in order to obtain strengthand pairing symmetry of the superconducting order parameter, which takes the form of a7eneralized eigenvalue problem − V BZ X b ˛ FS b Γ bb s/t ( k , k ) v b F ( k ) ∆ b ( k ) = λ ∆ b ( k ) . (11)Here, v bF ( k ) = |∇ (cid:15) b ( k ) | is the Fermi velocity at k in band b . The largest eigenvalue λ > for a given interaction kernel Γ bb s/t ( k , k ) , will lead to the highest transition temperature T c and the corresponding eigenfunction ∆ b ( k ) determines the symmetry of the gap. Theeffective lattice model obtained from the charge accumulation points has point group D h .The symmetry of the gap can thus by classified according to the irreducible representationsof D h that are listed in Table 1.The linearized gap equation (11) only accounts for the leading pairing symmetry atthe transition temperature T c of the superconducting phase. In the case of degenerateeigenvalues (e.g. d -wave instabilities { d xz , d yz } ) belonging to a two-dimensional irreduciblerepresentation, an arbitrary linear combination might be favored for T < T c . In order to findthe linear combination that is preferred by the system below the transition temperature, wecompute the free energy of the system F = E − T S = 1 N X k ,ν (cid:20) E ν ( k ) n F ( E ν ( k )) − | ∆ ν ( k ) | E ν ( k ) tanh (cid:18) E ν ( k )2 T (cid:19)(cid:21) + TN X k ,ν [ n F ( E ν ( k )) ln( n F ( E ν ( k ))) + n F ( − E ν ( k )) ln( n F ( − E ν ( k )))] . (12)Here, E ν ( k ) is the energy of the Bogoliubov quasi-particles resulting from diagonalizationof the BdG Hamiltonian H BdG = X k ,ν ψ † ν k (cid:15) ν ( k ) − µ ∆ ν ( k )∆ † ν ( k ) − (cid:15) ν ( − k ) + µ ψ † ν k = X k ,ν ψ † ν k [ δ ν ( k ) · τ ] ψ † ν k , (13)where δ ν ( k ) = ( < [∆ ν ( k )] , = [∆ ν ( k )] , (cid:15) ν ( k ) − µ ) T and τ are the Pauli matrices. In the ex-pression of the free energy Eq. (12), we only account for states at the Fermi surface ascontributions from k points far away from the Fermi surface are negligible (cid:15) ν ( k ) (cid:29) | ∆ ν ( k ) | .At the filling µ ≈ µ + 5 meV studied in the manuscript, the leading pairing symmetry isthe d -wave which belongs to a two-dimensional irreducible representation. To minimize the8 θ ( π )00.51.01.52 φ ( π ) sin( θ ) d xz + cos( θ ) e iφ d yz − . − . − . − . − . − . − . − . − . F / e V Figure
3. Free energy of the linear combination ∆ ν ( k ) = sin( θ ) d xz ( k ) + cos( θ ) e iφ d yz ( k ) corre-sponding to the leading pairing symmetry at µ ≈ µ + 5 meV. The system minimizes its free energyby choosing the linear combination d xz ( k ) ± id yz ( k ) . free energy of the system we make the ansatz ∆ ν ( k ) = sin( θ ) d xz ( k ) + cos( θ ) e iφ d yz ( k ) , (14)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 systemEq. (12). In Fig. 3 we show that the linear combination ∆ ν k ∝ [ d xz ( k ) ± id yz ( k )] =[ sin ( k x ) sin ( k z ) ± i sin ( k y ) sin ( k z )] is generally preferred for the given filling.singlet triplet s p z ( d x − y , d xy ) · p z ( d x − y , d xy )( d xz , d yz ) ( p x , p y ) f x ( x − y ) · p z f x ( x − y ) f y ( y − x ) · p z f y ( y − x ) Table 1. Pairing symmetries for the effective lattice model of thBN separated into contributions tospin singlet and triplet channel.
1] S. Plimpton, Journal of computational physics , 1 (1995).[2] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio,Physical Review B , 235137 (2018).[3] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott,Journal of Physics: Condensed Matter , 783 (2002).[4] A. N. Kolmogorov and V. H. Crespi, Physical Review B , 235415 (2005).[5] W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod, Nano letters , 6009 (2018).[6] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Physical review letters ,170201 (2006).[7] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano Letters , 804 (2010).[8] G. Kresse and J. Hafner, Phys. Rev. B , 558 (1993).[9] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes,Y. Bai, B. Kim, et al., Nat. Mater. , 1 (2020).[10] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Nano Lett. , 4934(2019).[11] P. E. Blöchl, Phys. Rev. B , 17953 (1994).[12] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters , 3865 (1996).[13] A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. , 073005 (2009)., 073005 (2009).