Moir´e Superstructures in Marginally-Twisted NbSe 2 Bilayers

,

The niobium (or tantalum) atoms of metallic TMDs host one fewer d-electron in the valence shell.The chemical potential then lies in the valence band, in contrast to semiconducting TMDs, allowing superconducting and charge-density wave (CDW) phases in monolayer and bulk [14][15][16][17][18][19].Here, interlayer effects are also found to be important, moderating the critical temperatures of the superconducting [20,21] and CDW [22,23] transitions, enabling van der Waals Josephson junctions through interlayer twisting [24,25] and leading to spatial variation of the Gibbs free energy of hydrogen absorption [26].Despite these fascinating properties, there have been no detailed theoretical studies of the structural relaxation of twisted bilayers of metallic TMDs, and the consequent effect of stacking on electronic properties.Here we address these questions for twisted bilayers of NbSe 2 using ab-initio density functional theory (DFT) modelling combined with multiscale analysis of lattice relaxation.
Ab-initio analysis of adhesion energy -DFT calculations were performed as implemented in the Quantum ESPRESSO code [27,28].Core electrons were approximated using Vanderbilt ultrasoft pseudopotentials [29] under the generalized gradient approximation (GGA), as parameterised by Perdew, Burke and Ernzerhof [30].A plane-wave cut-off of E w = 50 Ry and charge density cut-off of E ρ = 600 Ry were applied.[31] and Fermi-Dirac smearing of width σ = 0.01 eV is applied to aid convergence.Interlayer dispersion is implemented through the optB88-vdW functional [32][33][34][35].Band structure calculations were performed using E w = 80 Ry and a 21 × 21 × 1 k -point grid.In both band structure & adhesion parameterisation calculations, we fix the in-plane lattice constant to the bulk value.
To assess the preferred stacking order, bilayers of NbSe 2 were structurally relaxed in a variety of stacking configurations illustrated by sketches in Fig. 1a.Bilayer structures are divided into two classes: anti-parallel (AP), where a 180 o rotation of the top layer restores centrosymmetry, which is absent in a monolayer, and parallel (P) where monolayers have the same orientation and centrosymmetry is absent.Comparison of different stackings shows that the lowest-energy stacking polymorph is the metal-overlapping, MM (also called 2H a [36,37]) configuration (see Fig. 1b), for which DFT bulk lattice constants, a = 3.45 Å and c = 12.74 Å [38], and elastic moduli (see Supplemental Information (SI) section S2), are all found to be in good agreement with available experimental values [38][39][40].2H-stacking, which is preferred for semiconducting TMDs, is reasonably close in energy.While it has been found with a smaller areal density in CVD-grown NbSe 2 bilayers [41], it is a metastable state, whereas the MM configuration is most common and thermodynamically stable.We note that the bulk structure of NbSe 2 has frequently been referred to (and sometimes even modelled) as 2H stacking [42][43][44][45], therefore to avoid confusion, below we employ nomenclature of relevant stacking configurations as shown in Fig. 1.
In twisted bilayers, moiré superlattice reconstruction is fully determined by the stacking-dependent variation of adhesion energy between the constituent layers.Following the approach of Ref. [46] we employ an interpolation formula to fit the adhesion energy as a function of local disregistry, r 0 .This is parameterised in a coordinate system where XX stacking corresponds to zero displacement, i.e. r 0 = (0, 0), and interlayer distance d: The computed value is in a good agreement with the experimental value ω BM ≈ 34 cm −1 [47] , which validates the DFT-fitted shape of the adhesion energy (1) around d * .
Lattice reconstruction in twisted NbSe 2 bilayers is modelled using a multiscale approach implemented earlier in Ref. [46] for the analysis of semiconducting TMDs.This incorporates the microscopic expression for adhesion energy (1) and elasticity theory applied to the mesoscale strains across a long-period moiré superlattice, which reduces to minimization of total (adhesion and elastic) energy over the moiré supercell.In the elastic energy, characterised by elastic moduli λ and µ, we take into account only in-plane strains in top (t) and bot- ji , neglecting minor bending energies of layers (see SM in [46]) due to adjustment of local optimal interlayer distance with corresponding stacking across the moiré superlattice.To find the latter, we expand f (d) around its extremum (f 2 ), exponential functions in (1) up to linear order in d − d 0 , and, then, minimise the resulting expression as a function of d, to obtain an expression for the optimal interlayer distance for every stacking, d(r 0 ).Here, r 0 (r) = θẑ × r + u t − u b , where the first term describes contribution of geometrical twist between the layers, and in-plane displacements, u t/b , responsible for local deformations in t/b-layer, are found from minimisation to the total energy.
In Fig. 2 we display the reconstructed moiré superlattice resulting from minimisation of total energy in twisted NbSe 2 bilayers.Similarly to semiconducting TMDs, at [48] θ ≤ θ * P ≈ 2.9 • this results in the formation of arrays of triangular domains with rhombohedral stacking (XM and MX), separated by a domain wall network.Each domain wall in the network is a partial screw dislocation with dominating shear strain and a Burgers vector length a/ √ 3, shown by red arrows for several DWs merging into a DW network node.The magnitude of shear strain reaches 1% in the middle of domain wall, see Fig.

2, inset.
For AP-NbSe 2 bilayers with θ ≤ θ * AP ≈ 1.2 • , Fig. 3, lattice reconstruction leads to the expansion of the lowest energy MM domains, and the formation of a hexagonal superlattice domain wall network.The other highsymmetry (2H and XX) stackings occupy corners of the domains, linked by perfect screw dislocations, characterized by a stacking shift of a single translation vector along zigzag axes.We note that the critical angles, θ * P,AP are slightly higher than those for semiconducting TMDs [46] because of the softer NbSe 2 lattice.Electronic structure -The Fermi surface of an NbSe 2 monolayer is shown in Fig. 4a.There are three distinct hole pockets across the Brillouin zone: one Γcentred pocket and a pair of triangular K ± pockets.The dispersion around each pocket (i = Γ/K ± ) is where p = k − K i is Bloch state momentum relative to the pocket centre, m * i is effective mass, C i parameters account for trigonal and hexagonal warping, σ z is a Pauli matrix operating on spin, ϕ = arctan(k y /k x ) [49], and E 0,i is the energy difference between the band edge and the Fermi level.The shape of the metallic band is overall extremely similar to the conduction band in semiconducting TMDs, with similar effective masses m * i and a high degree of trigonal warping around the K-pockets (see Table 1) [16].Significant Ising spin-orbit coupling (SOC) is evident across the BZ, with splitting ∝ β(k)σ z .A map of β(k) across the entire BZ is displayed in Fig. 4a, with a maximum value of 78 meV deep in the K ± pockets.In the Γ-pocket SOC vanishes along Γ-M lines, β Γ = λ Γ |p| 3 cos(3ϕ), while it is approximately constant, β K = ±λ K in the K ± pockets [50].
Stacking modulation of interlayer hybridization results in variation of the electronic Fermi surface around each pocket.To quantify this, we relate DFT-calculated en- ergy eigenvalues across the entire Brillouin zone (BZ) to a momentum-dependent hybridization parameter α(k), assuming a simple layer-hybridized wavefunction, , where E M L (k) are monolayer energy eigenvalues.[51] Plots of this parameter calculated across the BZ are overlaid with the DFT-calculated Fermi surface in Fig. 4b.We observe that in centrosymmetric AP configurations, interlayer hybridization leads to a pair of spin-layer locked, layer-hybridized bands which cross the Fermi level in all stacking configurations.In contrast, broken centrosymmetry of the P-bilayers leads to four Fermi lines, which are essentially a pair of layer-hybridized copies of the monolayer bands.DFT calculations in mirror-reflected supercells also demonstrate the presence of a ferroelectric charge transfer, resulting in ≈11 meV potential energy drop between the layers (see SI section S3), which is notably smaller than for semiconducting TMDs [52].We note that this charge transfer in principle allows for a layer-dependent potential contribution to interlayer band splitting of P-bilayers, however, explicit incorporation of this term does not have a significant ef- fect on interlayer hybridization due to the small degree of charge transfer between the layers (see comparison between Fermi surface fitting with and without an explicit layer-polarization term, SI section S4).
In both bilayer orientations, there is a substantial six/three-fold modulation of the hybridization parameter around the Γ/K ± pockets, as shown in Fig. 4b (see also SI section S5).This modulation is proportional to the out-of-plane d 2 z -orbital component, leading to maxima along the Γ-M and K-M lines around each pocket.The resulting interlayer hybridization, and the associated shape of the Fermi lines, is moderated by the interlayer coordination of Nb atoms.Consequently, there are also distinct modulations of the interlayer hybridization at domain walls, which is evident as a two-fold "squeezing", of both the Γ and K-pockets, along distinguished directions of the Brillouin zone for the intermediate stacking configurations occurring at domain walls (see Fig. 4b, P-DW).For example, at the domain walls of a P-oriented moiré we find that there is significantly stronger hybridization of Γ-pocket electrons with crystal momenta parallel to the dislocation line compared to the perpendicular direction, while the opposite is the case for the K-pocket.A comparison of the degree of asymmetric warping, and the corresponding positions of Nb atoms in both layers, is shown in Fig. 5b for various intermediate stacking configurations of P and AP bilayers, which demonstrates qualitatively similar distortion between the two bilayers as a function of Nb-interlayer offset.Furthermore we note that, in both pockets, the degree of interlayer coupling can be lifted below or above that of the spin-orbit term, depending on pocket index and crystal momentum.(see Fig. 5c,d).
Using DFT data, we quantify effective interlayer hopping around each pocket by introducing an extra index which acts in layer space.This is encoded in a bilayer Hamiltonian, where d † /d are electron creation/annihilation operators, i is pocket index, p is Bloch state momentum relative to the pocket centre, s and l are spin and layer indices.η x,0 are Pauli matrices operating on the layer index, where the parameter a = 1, 2 for a AP, P bilayers respectively, and α i is the strength of hybridization around the corresponding pocket, where modulation is accounted for using appropriate periodic terms in α i .The hybridization term then takes the form, where ti is the average hybridization around pocket.We find that hybridization is approximately constant for momenta p ≈ |p F |, with angular variation ∝ ϕ incorporated via ti,n periodic terms.Further details of these terms in high-symmetry stackings are given in SI sections S6 and S7.
To facilitate simple comparison between different domains in the moiré superstructure, we analyse the average interlayer hopping, ti (the first term in Eq. 4), at each pocket as a function of interlayer disregistry.We numerically extract this parameter from DFT calculations by averaging interlayer splitting around each pocket along the high-symmetry directions (Γ-M, Γ-K, K-M, K-Γ).This procedure was repeated for relaxed bilayers with different interlayer offsets, and the extracted value of ti was found to fit to an expansion using only the first star of reciprocal lattice vectors with reasonable accuracy (see SI section S8 and S9).Fig. 5a shows this expansion for both pockets in AP-and P-bilayers.
Interestingly, we find distinct variation of interlayer hybridization between the Γ and K pockets in AP bilayers.For the former, which is most sensitive to the interlayer distance between Se atoms, it is maximal at XX corners, while interlayer hybridization at the K-pocket is strongest inside the lowest-energy MM domains.Interactions are weakest at the 2H corners for both pockets.An opposite trend was found for the moiré of a P-bilayer.In this case, the interlayer separation of chalcogen atoms is rather larger within domains due to trigonal interlayer coordination and there are correspondingly weaker interlayer interactions, which becomes larger at corners and domain walls where the interlayer Se distance is reduced.The overall variation is qualitatively similar at both pockets.
CDW Modulation -Lastly, we examine the ability of the reconstructed moiré superlattice to impact the relative phase of the CDW distortion in adjacent layers, as a function of the interlayer disregistry.CDW order in NbSe 2 monolayers induces a 3 × 3 reconstruction [20,[53][54][55][56].In a monolayer, we find two low energy triangular reconstructions of the lattice [57], which are characterised by one of two structural distortions, classified as "Hollow" or "Filled", based on the distortion pattern of atoms in the Nb sublattice [58] (see SI section 10), and that the Hollow phase has lowest overall energy of the two.
Both monolayer reconstructions are then used to create bilayer cells with different interlayer offsets of the CDW phase between the layers.The difference in DFT total energy between the normal (no reconstruction) and Hollow CDW phase on both layers, ∆E = E CDW − E N ormal , is extracted to quantify the realtive strength of the CDW phase between different structures.In agreement with the monolayer case, we find that structures with the Hollow reconstruction on both layers are lowest in energy.
We then calculate the relative strength of the CDW distortion at each of the 18 MX/XM (P) and 9 MM (AP) possible stacking configurations of the CDW in adjacent 3 × 3 reconstructed layers.Fig. 6a shows these three unique structural configurations in the AP case, where only the Nb-sublattice is shown to improve clarity (see SI section 10 for additional figures of the monolayer CDW reconstructions, and the different interlayer configurartions of these structures considered in this work.) We observe a notable trend in this energy difference with interlayer disregistry, shown in the bottom left (right) of Fig. 6b for the AP (P) bilayer.While it is possible to reach another low energy mutual orientation of the two CDW distortions in each layer via a full dislocation in the AP bilayer, this is not possible for P bilayers, where shifting by a partial dislocation always leads to a higher energy configuration.In a moiré superlattice, this suggests that the CDW order in adjacent layers will lock-in inside a domain.Glide of the lattice by a partial or full screw dislocation across domain boundaries will rigidly shift the relative orientation of the CDW on each layer, and, in order to attain a low-energy interlayer configuration, CDW discommensurations should nucleate at moiré dislocation boundaries [59,60], leading to CDW triplet domains in an AP-moiré [61] and isolated CDW domains in a P-moiré.
In conclusion, multiscale relaxation and electronic structure calculations were performed on twisted bilayer NbSe 2 .Compared to semiconducting TMDs, monolayer NbSe 2 is relatively soft, and significant relaxation begins at slightly larger twist angles.The resulting domain structure for AP-bilayers has hexagonal MM stacking domains with domain walls featuring XX stacking and seeds of 2H stacking in the alternating corners.The Fermi surface undergoes significant modulation across the moiré, with variation of the interlayer coupling at the K points on the order of 10-20 meV (25-55 meV) between different domains for an AP(P)-bilayer.Notably, there is significant anisotropy in the interlayer coupling is observed along domain walls.For AP-bilayers the interlayer coupling of metallic bands is strongest inside domains, in contrast to our finding regarding P-bilayers.In P-bilayers, the triangular MX/XM domains have smaller interlayer coupling of electron bands than along the triangular domain wall network.
A further observation arising from this work is that discommensurations [60] of the CDW phase should occur due to the rigid displacement of the relative orientation of the CDW in each layer which occurs at superlattice domain walls.CDW discommensurations are known to enhance the superconducting order in TiSe 2 [62][63][64].This suggests a promising application for moiré superlattices in metallic TMDs systems as a method to control the CDW, and potentially, other order parameters.The potential for distinct modulation of interlayer hopping at each pocket, in addition to control of the CDW phase, along with anisotropic hopping at dislocation boundaries, suggests that marginally-twisted NbSe 2 would be an interesting test system to further probe interlayer effects on correlated superconducting and CDW phases of the metallic TMDs [20].
Finally, we remark on possible implications of our results for Ising superconductivity in marginally-twisted bilayers.At small angles, the interlayer coupling in both pockets across a significant area of the moiré supercell is of a similar magnitude to the MM bilayer, for both twisted AP and P bilayers.In multilayer NbSe 2 [20] both the critical temperature, T c , and in-plane magnetic field H c of the superconducting transition show a significant dependence on the number of layers, suggesting the importance of interlayer electronic effects to the magnitude and form of the superconducting gap.If the breaking of spin-layer locking through increased interlayer hopping is the active mechanism in determining the superconducting transition temperature as a function of layer number, this suggests that the significant modulation of interlayer hopping by the moiré superstructure could be utilised to engineer superconductivity in specific domains of the superlattice.In an applied field, superconductivity would preferentially persist in regions of the moiré with smaller interlayer coupling, which are MX/XM domains (P-bilayer), or domain walls and corners (AP-bilayer).C1 (eV nm 2 ) C2 (eV nm 6 ) C3 (eV nm 10 ) A1 (eV/nm 2 ) A2 (eV/nm 2 ) ρ1 (nm) ρ2 (nm) d0 (nm) ε (eV/nm We model spin and interlayer splitting at each pocket with a Hamiltonian The interlayer term α varies with interlayer displacement between the layers, which also has periodic dependence on crystal momentum around a pocket, and the 3/6 cosine terms capture periodic dependence inside domains, while the 2/1 terms are due to nematic distortion of pockets in the vicinity of dislocations.The parameters tΓ,K are the average interlayer hybridization around a pocket.Spin-orbit terms are also included, such that spin splitting vanishes along Γ-M directions, and with fixed splitting in the K ± -pockets, S7. Full hybridization parameters around Γ and K pockets in domains and at domain boundaries.We have evaluated the average interlayer hybridization term, t, at individual pockets using DFT.We do this by performing band structure calculations along each of six individual high-symmetry lines in momentum space around each pocket, using fully-relaxed cells (at fixed disregistry -i.e.fixed planar coordinates).We then average the value of α(p) at the Fermi level along these directions.The obtained DFT data (as a function of disregistry) is fitted to a Fourier expansion using the first star of reciprocal lattice vectors around the lattice site r ′ , ti (r 0 ) = V 0 + g V (g)e ig•(r0−r ′ ) , (10) where the vectors g are the six TMD reciprocal lattice vectors in the first star, g n = ±g 1,2,3 , and

1 2HFIG. 1
FIG. 1. a) Side view of anti-parallel 2H stacking, antiparallel MM stacking and parallel MX stacking orders of NbSe2 bilayers.Additional stacking configurations are shown in Supplemental Information (SI) section S1. b) Adhesion energy versus interlayer separation for MM, MX, XX and ⟨W⟩ (configuration-averaged) stacking configurations.

FIG. 2 .
FIG. 2. Reconstructed moiré superlattice in P-NbSe2 bilayers at θ = 0.5 • .Red arrows show directions of Burgers vectors characterizing shift of atomic register across domain walls.On the upper inset we demonstrate distribution of the only nonzero shear strain across a single domain wall.Lower inset shows map of interlayer distance around XX area given by d(r0(r)).

FIG. 4
FIG.4.a) Two-dimensional Fermi lines for an NbSe2 monolayer.The colour overlay maps the magnitude of the spin-orbit splitting.The maximum magnitude of band splitting from spin-orbit effects only is ≈ 156 meV (β = 78 meV).b) Bilayer Fermi lines and magnitude of interlayer-hybridization parameter (α(k)) in selected stacking configurations.Splitting around K± pockets is determined by SOC effects, while Γ pocket splitting is largely due to interlayer interactions.DW corresponds to the P-bilayer "MX-XM" structure in SI S1.

FIG. 5
FIG. 5. a) DFT-calculated variation of the average interlayer hopping parameter ( ti)with disregistry, calculated by averaging splitting at the Fermi-level along six distinct crystallographic directions.Energy values are in meV.b) Modulation of the shape of the Γ-pocket Fermi lines for different intermediate stacking configurations in AP & P bilayers, and corresponding interlayer orientation of Nb atoms.c) Variation in the magnitude of the spin-orbit and interlayer contributions vs angular oientation relative to pocket centre, for the Γ and d) K pockets at a domain wall in a P-bilayer, extracted from DFT calculations.

FIG. 6
FIG.6.a) DFT-relaxed bilayer structures with different interlayer offsets of the CDW order between layers.For clarity, only the Nb sublattice is shown in red (blue) for the top (bottom) layer, and the centre of the CDW distortion is marked as the centre of the "large" triangle in each layer.b) Energy of the CDW phase, relative to the most favourable configuration, at different orientations of an AP-bilayer (left) and a P-bilayer (right), in units of meV/unit cell.

FIG. 7 .
FIG. 7. Plan view of the main stacking orders occurring in twisted NbSe2 domains.

TABLE I .
Fitted parameters of NbSe2 bands calculated at the Fermi level (brackets: close to centre of pocket).E0 is pocket depth, m is electron effective mass, C 3/6 accounts for trigonal/hexagonal warping at the K/Γ pocket, respectively, and λSOC is the magnitude of SOC.

TABLE IV .
AP-bilayer: Fit of αΓ at selected high-symmetry stackings.

TABLE V .
AP-bilayer: Fit of αK at selected high-symmetry stackings.
TABLE VI.P-bilayer: Fit of αΓ at selected high-symmetry stackings.

TABLE VIII .
) Parameters of Fourier expansion of the interlayer coupling as a function of disregistry.