Deconfinement of Mott Localized Electrons into Topological and Spin-Orbit Coupled Dirac Fermions

The interplay of electronic correlations, spin-orbit coupling and topology holds promise for the realization of exotic states of quantum matter. Models of strongly interacting electrons on honeycomb lattices have revealed rich phase diagrams featuring unconventional quantum states including chiral superconductivity and correlated quantum spin Hall insulators intertwining with complex magnetic order. Material realizations of these electronic states are however scarce or inexistent. In this work, we propose and show that stacking 1T-TaSe$_{2}$ into bilayers can deconfine electrons from a deep Mott insulating state in the monolayer to a system of correlated Dirac fermions subject to sizable spin-orbit coupling in the bilayer. 1T-TaSe$_{2}$ develops a Star-of-David (SoD) charge density wave pattern in each layer. When the SoD centers belonging to two adyacent layers are stacked in a honeycomb pattern, the system realizes a generalized Kane-Mele-Hubbard model in a regime where Dirac semimetallic states are subject to significant Mott-Hubbard interactions and spin-orbit coupling. At charge neutrality, the system is close to a quantum phase transition between a quantum spin Hall and an antiferromagnetic insulator. We identify a perpendicular electric field and the twisting angle as two knobs to control topology and spin-orbit coupling in the system. Their combination can drive it across hitherto unexplored grounds of correlated electron physics including a quantum tricritical point and an exotic first-order topological phase transition.

The interplay of electronic correlations, spin-orbit coupling and topology holds promise for the realization of exotic states of quantum matter. Models of strongly interacting electrons on honeycomb lattices have revealed rich phase diagrams featuring unconventional quantum states including chiral superconductivity and correlated quantum spin Hall insulators intertwining with complex magnetic order. Material realizations of these electronic states are however scarce or inexistent. In this work, we propose and show that stacking 1T-TaSe2 into bilayers can deconfine electrons from a deep Mott insulating state in the monolayer to a system of correlated Dirac fermions subject to sizable spin-orbit coupling in the bilayer. 1T-TaSe2 develops a Star-of-David (SoD) charge density wave pattern in each layer. When the SoD centers belonging to two adyacent layers are stacked in a honeycomb pattern, the system realizes a generalized Kane-Mele-Hubbard model in a regime where Dirac semimetallic states are subject to significant Mott-Hubbard interactions and spin-orbit coupling. At charge neutrality, the system is close to a quantum phase transition between a quantum spin Hall and an antiferromagnetic insulator. We identify a perpendicular electric field and the twisting angle as two knobs to control topology and spin-orbit coupling in the system. Their combination can drive it across hitherto unexplored grounds of correlated electron physics including a quantum tricritical point and an exotic first-order topological phase transition.

I. INTRODUCTION
Prospects of quantum information technologies have motivated an intense search for systems which intertwine topology and electronic correlations [1][2][3]. Strongly interacting and spin-orbit coupled electrons on the honeycomb lattice, as theoretically described by the Kane-Mele-Hubbard model [4], feature quantum spin Hall (QSH), Mott-Hubbard, collective magnetic and chiral superconducting states. While several weakly interacting QSH systems are now well-studied experimentally [2,5], material realizations of honeycomb Kane-Mele-Hubbard fermions with strong correlations and spin-orbit coupling are rare [6,7]. Here, we introduce a new "van der Waals engineering" platform to serve this purpose.
We show that stacking of 1T-TaSe 2 into bilayers can deconfine electrons from a deep Mott insulating state realized in the monolayer to a system of correlated Dirac fermions subject to sizable spin-orbit coupling. Central to this transition is the possibility of van der Waals materials to stack in different configurations. For a spe-cific honeycomb arrangement (Fig. 1), the kinetic energy associated with the electronic hopping t turns out to be of the same order of magnitude as the effective local Coulomb repulsion U . The system features therefore electronic correlations, which turn out to put the system right on the verge between QSH and correlated antiferromagnetic insulating states at charge neutrality and support chiral superconductivity under doping. We finally demonstrate that tuning the system via electric fields and twisting of the layers relative to each other sensitively affects the low-energy electronic structure in terms of emerging Dirac mass and spin-orbit coupling terms and leads to completely unexplored regimes of correlated electrons.

II. FROM A CORRELATED INSULATOR TO EMERGENT DIRAC FERMIONS
Layered group-V transition metal dichalcogenides (TMDCs) such as 1T-TaSe 2 or 1T-TaS 2 feature a low-temperature commensurate charge density wave (CCDW) where Ta-atoms are displaced into "Star-of-David" (SoD) patterns (Fig. 1a) [8,9]. In this phase, the SoDs form a triangular metal-to-insulator transition when entering the CCDW phase [9]. In 1T-TaSe 2 , the bulk remains conductive till lowest temperatures, while the surface exhibits a Mott transition around 250 K [10,11]. Recently, 1T-TaSe 2 has been fabricated down to monolayer thickness [12][13][14] and a pronounced thickness dependence of the electronic structure has been reported [14]. As bonds between the layers are mainly of van der Waals type, different stacking configurations are observed in experiments [15,16] and have a strong impact on the electronic structure [17][18][19].
We compare the CCDW state in the monolayer to a bilayer with "honeycomb stacking", where the SoD centers form a buckled honeycomb lattice (Fig. 1b). The bottom layer SoD centers form sublattice A and the top layer sublattice B. In this kind of stacking configuration, the Ta atoms of the bottom layer are approximately beneath the Se atoms of the top layer (Fig. 1c). While this bilayer stacking is not the most commonly observed one of the T-phase TMDCs [14,16], a corresponding stacking sequence has been found in transmission electron microscopy studies of 1T-TaS 2 [16] and turns out to be metastable in density functional theory (DFT) simulations of 1T-TaSe 2 (Supplementary Section S1).
To study the electronic structure of such engineered stackings, we combine ab-initio calculations in the frameworks of DFT (Supplementary Section S1) and the random phase approximation (RPA, Supplementary Section S3) with effective low-energy models, which we investigate with dynamical mean-field theory (DMFT) and twoparticle self-consistent (TPSC) many-body approaches (Supplementary Section S6).
In the CCDW phase, the band structure of monolayer 1T-TaSe 2 obtained with non-spin-polarized DFT is characterized by a single (Ta) flat band at the Fermi level [9,[17][18][19], which has a bandwidth of less than 20 meV (Fig. 2a left). Hence, the CCDW formation largely quenches in-plane hopping of the electrons. In the honeycomb stacked bilayer (with no twist), two dispersive bands with a bandwidth of the order of 200 meV emerge from the low-energy flat band of the monolayer (Fig.  2a right). Comparison of the mono-and bilayer bandwidths shows that interlayer hopping effects must dominate over intralayer hoppings by approximately an order of magnitude. In this sense, CCDW TaSe 2 bilayers are the exact opposite of graphene bilayer systems, since in the latter out-of-plane coupling is an order of magnitude weaker than in-plane hopping [21,22]. Notably, for the 1T-TaSe 2 honeycomb bilayer, the upper and lower lowenergy bands touch as Dirac points at the Brillouin zone corners K and K'. In the undoped system, these Dirac points are exactly at the Fermi level.
We next construct a Wannier Hamiltonian to describe the Dirac bands with one Wannier function for each SoD center, i.e. two Wannier orbitals per bilayer CCDW superlattice unit cell. The resulting nearest-neighbor hopping from a sublattice A site in the bottom to a neighboring sublattice B site in the top layer amounts to t = −34 meV and is the leading term of the Wannier Hamiltonian. There are further terms in the Wannier Hamiltonian, which are, however, at least an order of magnitude smaller than t.
The effective Hubbard interaction U for the SoD Wannier orbitals of CCDW TaSe 2 calculated in RPA (Supplemental Section S3) is U ≈ 130 meV, which is also in line with the experimental estimates in Ref. [14] and calculations for TaS 2 [18]. The ratio of hopping to Coulomb interaction is decisive in determining the strength and kind of electronic correlation phenomena taking place. Our calculations yield U t ≈ 3.8.
To study the resulting electronic correlations, we performed simulations of the Hubbard model for the nontwisted CCDW 1T-TaSe 2 bilayer in the framework of Dynamical Mean-Field Theory (DMFT) and the Two-Particle Self-Consistent (TPSC) approach (Supplemental Section S6) [23,24]. The quasi-particle weight Z shown in Fig. (2c) is a measure of the electronic correlation strength. In the temperature range T = 60 − 230 K, both, DMFT and TPSC consistently yield essentially constant Z ≈ 0.75. Our system is thus at intermediate local correlation strength and far away from the paramagnetic Mott-Hubbard transition taking place above U Mott t ≈ 8.2. The DMFT and TPSC studies of the full Hubbard model for the non-twisted CCDW 1T-TaSe 2 bilayer are, thus, in line with studies of the idealized Hubbard model on the honeycomb lattice [25].
The TPSC calculations give insight to spinfluctuations taking place in the system. The inverse intra-and inter-sublattice terms of the static magnetic susceptibility at wave number q = 0 as well as the antiferromagnetic correlation length ξ AFM are shown in Fig. 2d. We observe that antiferromagnetic fluctuations with alternating spin orientation between the two sublattices (A, B) are dominant and strongly enhanced at temperatures T ≲ 100 K. These fluctuations indicate that the non-twisted CCDW 1T-TaSe 2 bilayer is close to a quantum phase transition from a Dirac semimetal to an antiferromagnetic insulator, which occurs for ideal Hubbard honeycomb systems exactly in the range of U c t ≈ 3.6 − 3.8 [5,26,27].

III. SPIN-ORBIT COUPLING
The aforementioned magnetic correlation phenomena are sensitive to details of the low-energy electronic structure and spin-orbit coupling (SOC), which we discuss in the following based on a symmetry analysis. The space group of non-twisted CCDW 1T-TaSe 2 bilayer in the honeycomb structure is P3 (#147), comprising inversion symmetry and three-fold rotations C 3 around an axis perpendicular to the bilayer. Imposing also timereversal symmetry, every band must be two-fold degenerate. Therefore, SOC-induced qualitative changes of the band structures can occur near the Dirac points at K and K'. A corresponding ⃗ k ⋅ ⃗ p expansion (Supplementary Sec- where the pseudospin ⃗ S describes the sublattice degree of freedom, ⃗ σ acts on the electron spin, and τ = ±1 labels the valley (K, K') degree of freedom. This Hamiltonian comprises three contributions; the first contribution is a two-dimensional massless Dirac term, with the sublattice-pseudospin playing the role of the "spin" inherent to the Dirac equation. This term is analogous to the massless Dirac term in graphene [5,21,22]. SOC is responsible for the second and third contributions to H 0 : a valley-spin-sublattice coupling λ SOC = 0.74 meV, which is often called Kane-Mele spin-orbit term [28], and a sublattice-staggered Rashba term α R2 , which belongs to the R2 class according to the classification form Ref. [29]. A finite Kane-Mele term λ SOC opens a gap and turns the system described by H 0 into a quantum spin-Hall (QSH) insulator [28]. Importantly, the Kane-Mele term here is enhanced in comparison to its counterpart in graphene by two orders of magnitude [21,22,30], and corresponds to a temperature T ≈ 10 K, which is well accessible in experiments. Given that U t ≈ 3.8, the nontwisted CCDW 1T-TaSe 2 bilayer implements a material realization of the Kane-Mele-Hubbard model in a regime very close to the topological quantum phase transition from a QSH insulator to an antiferromagnetic insulator (Fig. 2e).
Dirac fermions and their topology are affected by different kinds of mass fields. An energy difference between electrons localized in sublattice A and B leads to a socalled Semenoff mass term M , which would enter the Hamiltonian H 0 (Eq. 1) in the form M S z . This term breaks sublattice invariance and thereby inversion sym-metry and leads to a transition from a QSH to a band insulator at M = λ SOC [28,31]. In the non-twisted CCDW 1T-TaSe 2 bilayer, the intrinsic Semenoff mass M 0 = 0 is required to vanish by symmetry [32]. However, vertical electric fields E z , as realizable in field effect transistor geometries (Fig. 1c), break inversion symmetry, translate into staggered sublattice potentials, and therefore corresponding extrinsic Semenoff mass contributions ∆M (Fig. 2b). DFT calculations (Supplementary Fig.  S3) yield the approximate relation ∆M ≈ eE z d ⊥ , where d = 6.3Å is the interlayer distance, e is the elementary charge, and ⊥ = 3.66 plays the role of an effective dielectric constant. The QSH to band insulator transition is reached for E z ≈ 0.5 mV/Å= 50 kV/cm (Fig. 2b), which is well within reach of experiments [33].
Taken together, our calculations show that the honeycomb-stacked non-twisted CCDW 1T-TaSe 2 bilayer is located in a region of the phase diagram ( Fig. 2e) with three different phases (QSH insulator, band insulator and antiferromagnetic insulator) coming together. Electron correlation is known to change the order of the QSH-band insulator transition from second to first order [20]. Contrary to the standard non-interacting QSH to band insulator transition, where the gap closes and reopens continuously with vanishing gap at the transition point, the QSH gap remains finite and the system changes discontinuously to a band insulating state at the transition point. Application of vertical electric fields in the system at hand represents hence a possibility to realize this exotic interaction-induced first order transition.

IV. TWISTED CCDW 1T-TASE2 BILAYERS
Layered van der Waals systems uniquely allow to realize different stacking configurations via twisting, i.e. relative rotations between the layers. General twist angles θ lead to incommensurate moiré patterns superimposed to the CCDW lattice. Arguably the simplest case of twisting is a rotation angle of θ = 180 ○ , which leads to a system with identical Bravais lattice but different symmetry of the supercell basis (Fig. 1c): in case of honeycomb stacking of the CCDW, the space group of the 180 ○ twisted structure of CCDW TaSe 2 bilayer is reduced to P3, meaning that inversion symmetry is lost with respect to the non-twisted case. The resulting band structure (Fig. 3a) is qualitatively similar to the non-twisted honeycomb case regarding the overall shape and width of the low energy bands. Thus, also the 180 ○ twisted case will be far away from the paramagnetic Mott transition. However, the low-energy band structure is markedly different in the 180 ○ twisted case. First, the conduction band is almost flat between K and M . Second, inversion symmetry breaking lifts band degeneracies: our DFT calculations reveal a staggered potential and an associated intrinsic Semenoff mass term of M 0 = 8.55 meV, which opens a gap at the K and K' points already without SOC and without external electric fields. Furthermore, addi-tional SOC terms are now allowed by symmetry (Supplementary Section S4) and completely lift the remaining band degeneracies except for the time-reversal symmetric points Γ and M.
Based on a symmetry analysis, we obtain the following low-energy model in the vicinity of K and K': The additional terms with respect to the non-twisted honeycomb structure Eq. 1 are the intrinsic Semenoff mass (M 0 ), the spin-valley coupling term (B), the Kane-Mele-Rashba interaction (λ R ), and a Rashba interaction belonging to the R1 class [29] giving rise to pure Rashba spin-polarization patterns (α R1 ). The last term (with coupling constant λ D ) can also be seen as an effective k−dependent magnetic field parallel to the z axis. Our DFT calculations yield M 0 = 8.55 meV, B = −1.85 meV, λ SOC ≲ 0.05 meV and λ R = 3.21 meV. These terms affect the dispersion and imprint an intricate sublattice and spin structure to the low-energy bands, which can be manipulated by vertical electric fields as shown in Fig. 3b.
Except for E z = −4.2 meV/Å (∆M = −6.68meV and M = M 0 + ∆M = 1.87meV), the system is always gapped. We calculated the Z 2 topological invariant for the noninteracting 180 ○ twisted bilayer in comparison to the nontwisted bilayer case as well as for two cases in between where the ratio B λ SOC is varied (Fig. 3c, Supplemental Section S5). In the non-twisted bilayer, the system is in a QSH state unless an extrinsic sufficiently large Semenoff mass term or an additional Rashba SOC term λ R are added. In the 180 ○ twisted case, the situation is very different regardless whether or not the intrinsic Semenoff mass term is compensated by an external electric field and regardless of λ SOC . Indeed, many changes in the SOC terms suppress the QSH state in the 180 ○ twisted bilayer: The comparably large Rashba λ R and the spin-valley coupling terms B and a strong reduction in λ SOC . Each of these alone is sufficient to suppress the QSH state. At vertical electric field E z = −4.2 mV/Å the gaps at K and K' close, and a parabolic band touching point emerges.
While it is clear that there will be tendencies towards interaction-induced (quasi)ordered phases as well, here, the kind of ordering is likely different from the nontwisted case but largely unexplored. The band touching at E z = −4.2 mV/Å implements a situation similar to saddle points in a two-dimensional dispersion, where already arbitrarily weak interactions would trigger different kinds of magnetic or excitonic instabilities [34]. How these instabilities translate into the intermediately correlated and strongly spin-orbit coupled case of 180 ○ twisted TaSe 2 is a completely open question. For B = 0 the system shows the symmetric "onion"-like shape similar to Ref. [28]. When B increases, the QSH region shrinks and disappears when B = λ SOC . For B > λ SOC , the system behaves as a trivial band insulator. The yellow arrows in the leftmost and rightmost panels indicate the path in phase space accessible by varying the vertical electric field Ez.

V. CONCLUSIONS AND OUTLOOK
The field of "'twistronics"' with materials like bilayer graphene is based on the idea that weak interlayer coupling can flatten highly dispersive bands and thereby boost electronic correlations [35][36][37]. The system introduced here takes the opposite route of deconfining formerly Mott localized electrons. This approach should be generally applicable to interfaces of Mott localized electrons under two conditions: the interlayer coupling should substantially exceed the in-plane one and at the same time define a connected graph linking all sites of the system. Possible example systems range from stacking faults in the bulk of CCDW layered Mott materials [16] to molecular systems [38].
Especially the twisting degree of freedom opens new directions to experiments. Since interlayer hopping is the dominant kinetic term in deconfined Mott systems like bilayers of CCDW 1T-TaSe 2 , we expect incommensurability effects to be much more pronounced than in twisted graphene systems [35]. θ = 30 ○ twisted CCDW 1T-TaSe 2 bilayer should realize a quasicrystal with twelvefold rotation symmetry and provide an experimental route to correlated electrons and emerging collective states in a quasicrystalline environment.
The specific case of CCDW 1T-TaSe 2 bilayer demonstrates how deconfinement of Mott electrons leads to exotic states of quantum matter: the non-twisted bilayer approaches the quantum tricritical region of competing quantum spin Hall, trivial band insulating, and antiferromagnetic insulating states. At 180 ○ twist angle novel kinds of electrically controllable band degeneracies with associated many-body instabilities, hypothetically of excitonic type, emerge. Clearly, the phase space for manipulating deconfined Mott electrons is high dimensional. We here identified the combination of twist angle and perpendicular electric field as decisive for TaSe 2 bilayers. Further means to control emerging electronic states include dielectric engineering [39] and charge doping. Our calculations showed that the non-twisted CCDW 1T-TaSe 2 bilayer in honeycomb stacking approximates the (Kane-Mele) Hubbard model on the honeycomb lattice with U t ≈ 3.8 very well. In this regime, doping the system away from the Dirac point towards the van Hove singularity is expected to lead to chiral superconductivity [40][41][42][43], most likely of d + id-type.

S1. DENSITY FUNCTIONAL THEORY
The Density Functional Theory (DFT) [44,45] calculations presented here are obtained using the Vienna ab-initio simulation package (VASP) [46,47] with the generalized gradient approximation of Perdew, Burke, and Ernzerhof (GGA-PBE) for the exchange-correlation functional [48,49]. We obtain the electronic structure of undistorted bilayer 1T-TaSe 2 , as well as for other group V transition metal dichalcogenides (TMDCs) 1T-NbSe 2 and 1T-TaS 2 . We calculate the potential landscape for various possible stackings in undistorted bilayer 1T-TaSe 2 using a Γ-centered k-mesh of 15×15×1, and taking into account van der Waals (vdW) corrections within DFT-D2 method of Grimme [50], see Fig. S1. We show the results for non-and 180 ○ twisted cases. For non-twisted, Ta on top of Ta (stacking I) is the most stable configuration, but Ta on top of one Se (stacking III, Se above Ta plane) results to be a local minima. For 180 ○ twisted, the most stable situations are given for Ta on top of Se (stackings III and V are degenerated), so that the ideal honeycomb stacking can be easily realized.
We calculate the commensurate √ 13 × √ 13 charge density wave (CCDW) structure by relaxing first the monolayer (without including vdW corrections), where a = 12.63Å according to Ref. [14]. The ionic relaxation is done using the conjugate gradient algorithm, and the procedure stops when the total free energy reaches the resolution 0.02 eV/Å. We block off-plane displacements for Ta atoms, while allowing for in-plane displacements. Se atoms are allowed to freely relax in all three directions. We include then a second layer and optimize the interlayer distance. We find d = 6.3 A for 1T-TaSe 2 . We further relax the CCDW bilayer by following the same procedure described for the monolayer. In the bilayer, off-plane displacements for the Ta atoms in both layers are also blocked. The Γ-centered k-mesh was set to 9×9×1 during all the ionic relaxation and interlayer distance optimization procedures.
For the non-collinear magnetic calculation, i.e. when including the spin-orbit coupling (SOC), we set the net magnetic moment to zero in all atoms of the unit cell, and use a Γ-centered k-mesh of 6×6×1.

S2. WANNIER TIGHT-BINDING MODEL
The relevant subspace B for the low-energy bands of (distorted) CCDW 1T-TaSe 2 monolayer and bilayer contains only d z 2 orbitals from Ta A-type atoms, i.e. the central star-of-David atoms (see Fig. 1). We then construct a minimal tight-binding model using Wannier90 code [51]. if SOC is not included in the calculation, and Γ-centered k-mesh of 6×6×1 when SOC is included. In the latter case, an inner energy window covering states at K (and K') and M is considered, while leaving out the ones at Γ. We show in Fig. S2 the comparison between the DFT bands and our Wannier tight-binding model. Our model captures very well the low-energy bands of non-twisted and twisted CCDW 1T-TaSe 2 bilayers.
The parameters resulting from the k ⋅ p expansion of the non-twisted case and the 180 ○ twisted case as entering Eq. (1) and (2) of the main text are shown in Fig. S3.

S3. ESTIMATION OF THE SCREENED HUBBARD INTERACTION U
We estimate the local Hubbard interaction U for the flat bands around the Fermi level in the CCDW 1T-TaSe 2 bilayer from the ab-initio calculation of the screened Coulomb interaction, using the random phase approximation (RPA) for the undistorted bilayer. We follow a similar procedure as the one described in Ref. [52], which we summarize below: • We initially calculate the Wannier90 tight-binding model for the three low-energy Ta bands C, whose orbital character is mostly {d z 2 , d x 2 −y 2 , d xy } (see Fig. 2(a) in Ref. [52]) in the undistorted monolayer 1T-TaSe 2 .
• In the CCDW 1T-TaSe 2 bilayer, the d z 2 orbitals from Ta A-type atoms have the largest contribution for the bands around the Fermi level. Thus, for each q, we consider only the tensor element W (q) ≡ W αααα (q) with α = d z 2 .
• Then, the local Hubbard interaction U in a single star of David is calculated by averaging over the d z 2 orbital weight from each Ta atom (labeled by w d z 2 (R)) in the star of David: where U (R) is the discrete Fourier transform of W (q).

S4. k ⋅ p MODEL OF LOW ENERGY BANDS
We derive the effective model describing the low-energy band structure around the K and K' points by using the k ⋅ p perturbation theory [53,54] for the perturbed Hamiltonian H = H I + H II + H III , including Pauli-type spin-orbit interaction, where: The non-zero matrix elements of such perturbed Hamiltonian can be determined from group theory and exploiting the transformation properties of the one-electron wave functions under the symmetry operations of the high-symmetry point K. The knowledge of the little group of K allows to define a set of linear equations [55]: where h is the order of the group, u µ ν represents the ν-th component of the basis function of a representation µ and O α β is an operator that transforms like the β-th component of the basis function of a representation α, while i D j ′ j (R) is the (j ′ j) element in the matrix representative of the group element R in the i-th representation. The non-twisted (180 ○ -twisted) CCDW bilayer belongs to the P3 (P 3) space group; both structures display a three-fold rotation C 3 around an axis perpendicular to the bilayer, while the non-twisted stacked bilayer is also centrosymmetric, the two layers being inversion partners. At point K only C 3 symmetry is preserved, and non-relativistic bands can be grouped in a one-dimensional single-valued irreducible representation (IR) K 1 and in two-fold degenerate K 2 K 3 IR for the space group P3; when the symmetry is lowered to P 3, bands belonging to K 2 and K 3 are not degenerate anymore. The single-valued IRs, the corresponding characters and the basis functions are listed in Table S1. By introducing a general operator π with components each trasforming as K 1 , K 2 and K 3 , respectively, one can use Eq. (S5) to identify its symmetry-allowed non-zero expectation values on the basis functions φ 2 , φ 3 spanning the K 2 K 3 IR: a purely imaginary nearest-neighbor spin-conserving SOC hopping, whose coupling constant can be parametrized by λ D . We obtain all effective parameters entering in the low-energy models around K and K from a first-order Taylor expansion in k-space of the Wannier tight-binding models described in Supplementary Section S2. Here ω = e i 2π 3 , and P (A) stands for polar (axial) vector. Notice that polar and axial vectors transform in the same way under the symmetry operations of C3 point group.

S5. TOPOLOGY
The topological properties of the twisted and non-twisted TaSe 2 bilayers without the effects of correlation are studied using the Wannier charge center evolution, as described in [57]. Using the Wannier90 tight-binding Hamiltonians as input, we determine the Z 2 invariant to assess whether they describe a trivial or a quantum spin Hall insulator.
For the non-twisted structure we find that the Z 2 invariant equals 1 for all values of the electric field between 0 and 0.43 mV/Å while Z 2 =0 for all other E-field strengths. For the 180 ○ twisted bilayer we find Z 2 =0 for all calculated E-field strengths instead.
In order to understand why the 180 ○ twisted bilayer does not show any topologically non trivial phases we analyze the influence of the terms entering the Hamiltonian H 180 ○ from Eq.2 of the main text. We proceed from the Wannier Hamiltonian of the 0 ○ twisted bilayer and artificially add and vary parameters occurring in H 180 ○ from Eq. (2). Specifically, we consider a parameter space, where we vary the Semenoff mass M , the Rashba-spin orbit term λ R and the spin-valley coupling term / valley Zeeman B. To determine the topological properties of non-interacting Hamiltonians in this parameter space, it turns out to be numerically efficient to track the minimal gap in the band structure and to identify connected regions of the parameter space with strictly non-zero gap, which are bounded by a submanifold of the parameter space with zero gap. If one point inside has Z 2 ≠ 0 (Z 2 = 0) the whole region is topologically non-trivial (trivial). The resulting topological phase diagrams are shown in Fig. (3c) of the main text.
It is useful to compare our system to the well-known case of the ideal Kane-Mele model [28]. The characteristic shape of the phase-diagram with the Rashba coupling (λ R ) on the y-axis and the staggered potential M (referred to as "Semenoff mass") on the x-axis in the ideal Kane-Mele model resembles that of an onion [28]. We observe this behavior for the non-twisted case, see leftmost panel of Fig. (3c). The yellow horizontal line represents the trajectory of our Hamiltonian upon changing the electric field E z . The main effect of E z is to change M according to M ≈ M 0 + eE z d ⊥ (see Fig. S3) and hence E z drives the 0 ○ twisted system from inside the topological region to the outside.
Upon twisting, we reduce the symmetry from D 6h to C 3v consequently switching on various terms in H 180 ○ , as described in Ref. [30], among which the most important ones are the previously-mentioned Semenoff mass M , Rashbaspin orbit λ R and the spin-valley coupling B. The main effect of B is to shift two of the bands having the same sublattice but different spin character at the K-point towards each other. This yields a deformed onion as topological region, where the topological region on the λ R axis is reduced. (See Fig. (3c) of the main text for the case of B λ SOC = 0.9.) Eventually when B = λ SOC the topological non-trivial region is completely suppressed and only a vertical band touching line remains, which, however, does not separate a trivial from a non-trivial region. The latter case is similar to the phase diagram of the 180 degree twisted TaSe 2 bilayer. The effects of B and λ R can, thus, qualitatively explain the differences observed for the topological phases of the non-twisted and 180 ○ twisted structure. In the 180 ○ twisted structure, λ SOC is additionally strongly reduced, which further contributes to suppressing the QSH phase. Following this line of argumentation the electric field in the 180 ○ twisted case induces a gap closing and a reopening but does not induce a topological phase transition.
In the non-twisted case, where we find QSH states in absence of interactions, we study the impact of correlations on the topological phase diagram within the TPSC approach. These calculations confirm the generic shape of the schematic phase diagram shown in (Fig. 2c).
It is known [20] that sufficiently strong Coulomb repulsion U can change the nature of the topological phase transition between QSH and band insulator (Fig. 2c)  the gap (solid line in Fig. 2c). Since the electric field E z directly controls M , one can possibly tune the 0 ○ system through this exotic first-order transition in experiments by varying E z .

S6. DMFT AND TPSC
For non-twisted CCDW 1T-TaSe 2 bilayer, we constructed tight-binding-Hubbard Hamiltonians of the type from a Wannierization of the ab-initio DFT band structure (Supplementary Section S2) and estimated the on-site repulsion U to be about 130 meV by means of RPA calculations (Supplementary Section S3). We studied these effective low-energy models from a many-body perspective to judge the type of correlations in the system. In dynamical mean field theory (DMFT) the lattice Hamiltonian is mapped onto a self-consistently determined single impurity problem, solved -in our case -within numerical exactly quantum Monte Carlo in the hybridization expansion flavour (CT-HYB) (for a review, see [58]). The resulting sublattice-resolved self-energy is local (k-independent) but it contains frequency-dependent non-perturbative corrections beyond Hartree-Fock to all orders and can account for Mott Hubbard metal insulator transitions. All DMFT calculations have been performed using w2dynamics [59]. The double counting is accounted for using the fully-localized limit.
For two-dimensional systems it is important to estimate non-local effects at the level of the self-energy, not included in DMFT. To this goal, we apply the Two-Particle Self-Consistent (TPSC) method [60], which produces accurate results in the weak-to-intermediate coupling regime, if compared to lattice quantum Monte Carlo calculations in the single band Hubbard model. For non-twisted CCDW 1T-TaSe 2 bilayer, which is modelled by a multi-band system we use the multi-site formulation of TPSC [23] while neglecting the Hartree term to avoid double counting of correlation effects already accounted for in DFT. Moreover, to be able to apply TPSC to this system we project out spin offdiagonal terms and take only the diagonal spin-up contributions from DFT while still assuming a paramagnetic state. The combination of TPSC accounting for the k-dependence of Σ and DMFT, in which we can include all off-diagonal terms between spin-orbitals and access antiferromagnetic ordering at strong coupling consitutes a powerful tool to determine the many-body nature of 1T-TaSe 2 .
As shown in the main text (Fig. 2), TPSC and DMFT consistently yield quasiparticle weights Z ≈ 0.75 for temperatures in the range 60-230K for non-twisted CCDW 1T-TaSe 2 bilayer. Also double occupancies ⟨n ↓ n ↑ ⟩ ≈ 0.15 agree well between TPSC and DMFT in this temperature range. These results consistently put non-twisted CCDW 1T-TaSe 2 bilayer at intermediate local correlation strength and clearly far away from a paramagnetic Mott Hubbard transition.
The onset of non-local correlations can be inferred from the temperature dependent enhancement of antiferromagnetic susceptibilities and correlation lengths, which we obtained with TPSC and show in the main text (Fig. 2). These correlations manifest in sublattice off-diagonal contributions to the self-energy shown in Fig. S4.