Unconventional orbital-charge density wave mechanism in transition metal dichalcogenide 1T-TaS2

The transition metal dichalcogenide 1T-TaS2 attract growing attention because of the formation of rich density-wave (DW) and superconducting transitions. However, the origin of the incommensurate DW state at the highest temperature (~ 550 K), which is"the parent state"of the rich physical phenomena, is still uncovered. Here, we present a natural explanation for the triple-q incommensurate DW in 1T-TaS2 based on the first-principles Hubbard model with on-site U. We apply the paramagnon interference mechanism that gives the nematic order in Fe-based superconductors. The derived order parameter has very unique characters: (i) the orbital-selective nature, and (ii) the unconventional sign-reversal in both momentum and energy spaces. The present study will be useful for understanding rich physics in 1T-TaS2, 1T-VSe2, and other transition metal dichalcogenides.

The transition metal dichalcogenides (TMDs) provide a promising platform of the exotic low-dimensional electronic states with strong electron correlation. Among the TMDs, 1T-Ta(S,Se) 2 exhibits very interesting electronic properties, such as the metal-insulator transition with the David-star formation as well as exotic superconductivity. The electronic states are easily controlled by the gate-electric-field carrier doping [1][2][3], by changing the dimensionality [4][5][6], and by applying the pressure [7,8] and picosecond laser pulses [9].
In 1T-TaS 2 , at ambient pressure, the incommensurate charge-density-wave (IC-CDW) appears as the highest transition temperature at T = T IC ≈ 550K [10]. With decreasing T , the IC-CDW changes to the nearly commensurate (NC) CDW at T NC ≈ 350K, and finally David-star commensurate (C) CDW appear at T C ≈ 200K successively. This rich multistage CDW transition is suppressed under pressure, and the superconductivity emerges at T SC 10K. These exotic ordered states emerge under the IC-CDW state. That is, the IC-CDW is the parent electronic state of rich physics in 1T-Ta(S,Se) 2 [11][12][13][14][15][16][17][18][19]. Nonetheless, the understanding of the origin and nature of the IC-CDW is very limited at present.
Although phonon-driven CDW inevitably causes sizable lattice distortion (LD) in proportion to the transition temperature, the LD below T IC is much smaller than that below T C in 1T-TaS 2 [20,21]. Considering the importance of electron correlation in 1T-TaS 2 , it is important to investigate the electron-correlation-driven IC-CDW mechanism, although the derivation of "nonmagnetic IC-CDW order" is a very difficult theoretical problem. In fact, magnetic order is always obtained based on the Hubbard models with on-site U within mean-field theories. Thus, one may consider the existence of large off-site bare interactions (such as the off-site Coulomb interaction and RKKY interaction [22]) comparable to U . Therefore, the study of IC-CDW is very important to uncover the real Hamiltonian for 1T-TaS 2 , based on which the multistage CDW transition and the superconductivity should be studied.
In this paper, we present a natural explanation for high-T IC IC-CDW, which is "the parent electronic states" of the exotic multistage CDW and superconductivity in 1T-TaS 2 . The predicted IC-CDW is the correlationdriven "unconventional CDW", in which the CDW order parameter possesses strange orbital-momentum-energy dependences, in analogy to the unconventional superconductivity. The wavevectors of the IC-CDW state coincide with the Fermi surface (FS) nesting vectors q = q 1 , q 2 , q 3 in Fig. 1 (a) [40,41]. In addition, with the aid of the Ginzburg-Landau (GL) theory, we reveal that the tripleq CDW state is stabilized. This study provides necessary knowledge in resolving the mysterious C-CDW state [42][43][44][45][46][47][48][49]. This theory will be useful for understanding rich CDW states in 1T-TaS 2 , 1T-VSe 2 , and other TMDs. binding model of 1T-TaS 2 , H 0 , composed of five 5d 1 orbitals of Ta-ions and six 3p orbitals of S ions, using the Wien2k and Wannier90 softwares. The d-electron eigenfunctions in the S 6 octahedron under the trigonal distortion are composed of one a 1g -orbital, two e ′ g -orbitals, and two e g -orbitals. In this paper, we assign a 1g , e ′ g , e g -orbitals as orbitals 1, (2, 3), (4,5) in order. The FSs are mainly composed of orbitals 1-3. The wavefunctions of orbitals 1-5 and the model Hamiltonian based on experimental crystal structure [50] are given in the Supplemental Materials (SM) A [51].
We also introduce the on-site Coulomb interaction Hamiltonian, H U . It is composed of the intra (inter) orbital interaction U (U ′ ), and the exchange interaction J = (U − U ′ )/2. Below, we fix the ratio J/U = 0.1. (The obtained results are similar for J/U < 0.2.) These interaction are included into the spin (charge) channel interaction matrixΓ s(c) ; see SM A [51] for detail.
However, the IC-CDW without magnetization cannot be explained by the RPA because charge Stoner factor is always smaller than α S in the RPA [26]. To explain the IC-CDW state, we study the charge-channel susceptibility χ DW (q) due to the higher-order vertex corrections (VCs), based on the density-wave (DW) equation method [27,33]. (The VCs are dropped in the RPA.) Figure 1 (b) is an Aslamazov-Larkin (AL) type VC for χ DW (q), which is proportional to the convolution of paramagnons C q ∼ p χ s (q 1 +p)χ s (p). We will show that the AL type VC induces the IC-CDW order at q = q i (i = 1, 2, 3), since C q is large at q = q i due to the momentum conservation in Fig. 1 (c).
Here, we introduce the linearized charge-channel DW equation [27,33]: where k ± ≡ k ± q/2, k ≡ (k, ǫ n ) and p ≡ (p, ǫ m ) (ǫ n , ǫ m are fermion Matsubara frequencies). L ≡ (l, l ′ ) and M i represents the pair of d-orbital indices. λ q is the eigenvalue and f L q (k) is the Hermite form factor. The former represents the instability of the DW fluctuations at wavevector q, which reaches unity when the long-range order is established. The latter is the general chargechannel order parameter: f l,l ′ q (k) = σ c † k+,l,σ c k−,l ′ ,σ − c † k+,l,σ c k−,l ′ ,σ 0 . The DW equation (1) is interpreted as the "charge-channel electron-hole pairing equation" with the pairing interaction I L,M q (k, p). I L,M q at q = 0 is given by the Ward identity −δΣ L (k)/δG M (k ′ ) that is composed of one singlemagnon exchange Maki-Thompson (MT) term and two double-magnon interference AL terms; see Fig. 2 (a). Here, we set T = 0.04eV and U = 3.87eV (α S = 0.85). (Since d-electron weight in the DOS at the Fermi level is about 70%, U is reduced to ∼ 2.9eV in the d-orbital Hubbard model.) The analytic expression of the MT and AL terms are explained in the SM C [51]. The AL terms are proportional to the convolution of paramagnons, C q , so they become important when α S approaches unity [26,33]. Their essential role has been revealed by the functional-renormalization-group (fRG) study in which higher-order VCs are produced in an unbiased way [30,32]. In contrast, the MT term is important for the superconducting gap equation and for the transport phenomena [52].  shows the obtained eigenvalue of the DW equation λ q for T = 0.04eV and α S = 0.85; We see that λ q has six peaks at the nesting vector (q = ±q n , n = 1, 2, 3). As shown in Fig. 2 (c), the eigenvalue λ q1 increases with U , and it exceeds α S and reaches unity for α S 0.90. The susceptibility of the DW is given as χ DW (q) ∝ 1/(1 − λ q ). Thus, the obtained results are consistent with the IC-CDW order at q = ±q n without magnetization in 1T-TaS 2 .
We stress that T IC is enlarged by phonons in case that the form factor of the electron-phonon interaction is equal to the DW form factor in symmetry [53]. As discussed in Ref. [53], the "total DW susceptibility" is given as χ DW tot (q) = χ DW (q)/(1 − 2gχ DW (q)), where g(> 0) is the phonon-induced attraction. In fact, In Fe-based superconductors, the nematic transition temperature is raized by ∼ 50K due to the B 1g phonon mode. Similar phonon-assisted increment of T IC is expected in 1T-TaS 2 .
FIG. 3: (a) Real-part of obtained form factor f 1,2 q 1 (k), which is the most significant component for the IC-CDW at q = q1. (b) Fourier transformation of the form factor f 1,2 q 1 (r). The signal at r = 0 represents the local orbital+charge order with respect to orbitals 1 and 2, and the signal at r = 0 exhibits the non-local bond order. The strongest bond order is at r = (0, ±2) marked by pink broken circles.
Here, we discuss the nature of the form factor at q = q i . In Fig. 3 (a), we show the obtained real-parts of f 1,2 q1 (k) at the lowest Matsubara frequency. It is the largest intra-orbital component of the form factorf q1 (k). The largest intra-orbital component is f 1,1 q1 (k), which is exhibited in Fig. S3 It also induces the finite charge order (δn) since E 1 = E 2 in the present model, as shown Fig. S1 (b) in SM A [51]. Interestingly, Fig. 3 (a) has sign reversal, which is very different from the conventional CDW with nearly constant form factor. Similar "sign reversing form factor" is observed in the nematic phase in FeSe by ARPES measurement [54], and it is satisfactorily explained as the paramagnon interference mechanism [27].
Next, we examine the real space DW structure with q = q n . For this purpose, we perform the Fourier transformation of the form factor: where r i is the real space position of site i, and θ is arbitrary phase factor. Here,f l,m (r i , r i ) represents the local charge and/or orbital order at r i , andf l,m (r i , r j ) with i = j is the bond order (i.e., the modulation of the hopping integral) between r i and r j . Figure 3 (b) is the obtained f 1,2 q1 (r): Its large magnitude at r = 0 represents the orbital order with respect to |1 ± |2 . In addition, f 1,2 q1 (r) shows large values for r = 0. Thus, both orbital order and bond order coexist in the obtained IC-CDW state. The trace of the form factor, f tr q1 (r) ≡ l f l,l q1 (r), is shown in Fig. S3 (f) in the SM B. The large value of f tr q1 (r) at r = 0 means the emergence of the local charge order. To summarize, the IC-CDW in 1T-TaS 2 is identified as the combination of the charge/orbital/bond order in the present study.
In the obtained IC-CDW state, essentially all elements of the form factor f l,m qn in the a 1g +e ′ g orbital space (l, m = 1 ∼ 3) are large. This fact indicates that all a 1g + e ′ gorbital states cooperatively magnify the eigenvalue. In order to identify the major order parameter, we solve the DW equation (1) for q = q 1 by considering only f 1,2 q1 and f 2,1 q1 , by setting other elements zero. In this case, the obtained form factor f 1,2 q1 (k) is very similar to Fig.  S3 (a), and λ q1 is reduced just to 86% from Fig. 2 (b). In contrast, the eigenvalue becomes very small if only diagonal elements, f l,l q1 , are considered. Therefore, offdiagonal form factor f 1,2 q1 is the main order parameter of the IC-CDW state.
In the present mechanism, the orbital+charge DW order due to f 1,2 q1 = 0 originates from the interference among spin fluctuations at q = q 2 and q = q 3 , as depicted in Fig. 1 (c). Mathematically, f 1,2 q1 is derived from the kernel I L,M q1 . However, local net charge order (δn) is energetically unfavorable due to the mean-field term (∼ U δn) in the first term of Fig. 2 (a). In fact, in Fe-based and cuprate superconductors, the bond and orbital orders with δn = 0 appear since they are not prohibited by on-site U [27,32]. To understand why net charge order is obtained in 1T-TaS 2 model, we examine the energy-dependence of the form factor. Figure 4 (a) shows the frequency (ǫ n ) dependence of Ref 1,2 q1 (k, ǫ n ) near the van-Hove singular point k = (0, 0.5π). Similar sign reversal appears in other elements of the form factor. This sign reversing form factor is very similar to the sign reversing gap function in the s-wave superconductors with U = 0, known as the "retardation effect" that drastically reduces the depairing by U . Thus, net charge order in the IC-CDW state in 1T-TaS 2 , which is very unusual in metals with large U , is stabilized by the retardation effect. To summarize, the predicted "unconventional CDW state" in 1T-TaS 2 is characterized by the orbital selective form factor with strange sign reversals in the momentum and energy spaces.
Finally, we explain the "triple-q CDW state" in 1T-TaS 2 , which is the uniform coexisting state of the three order parameters with q = q 1 , q 2 , q 3 . For this purpose, we construct a simple Ginzburg-Landau Free energy for the CDW order Φ(k) = (η 1 f q1 (k), η 2 f q2 (k), η 2 f q2 (k)), where η i is the real order parameter and f qi (k) is the normalized form factor. Then, the free energy is given by The fourth-order coefficients b 0 and c 0 can be calculated by using the Green functions and form factors; see Fig. 4 (b). The derivation of a 0 , b 0 and c 0 in addition to the third order term F (3) = d 0 η 1 η 2 η 3 is given in the SM D [51]. Here, the single-q state and the triple-q state correspond to (η 1 , η 2 , η 3 ) = (η, 0, 0) and (η, η, η), respectively. It is easy to show that the triple-q condition is c 0 /b 0 < 2 if d 0 is negligible. As we show in the SM D [51], the ratio c 0 /b 0 = 1.1 is obtained by using the form factors in the present study. (f lm q1 is given in Fig. 3 (a) and Figs. S3.) Thus, the present IC-CDW state satisfies the triple-q condition. In contrast, in case of conventional CDW form factor with f l,m qn = δ l,m , the obtained ratio c 0 /b 0 = 3.2 does not satisfy the triple-q condition; see SM D [51]. Thus, the obtained unconventional form factor due to the AL processes is indispensable to explain the triple-q IC-CDW state in 1T-TaS 2 , which is schematically shown in Fig. 4 (c).
Here, we calculate the electronic states below T IC based on the 4 × 4 cluster tight-binding model with finite DW order given in Eq. (2). We make the wavevector of the DW order q 1 = (0.5π, 0) by introducing 20% holedoping. Although the folded FS under the CDW state is very complex, it can be "unfolded" to the original BZ by restoring the translational symmetry of the spectral function [55]. Figures 5 (a) and (b) show the obtained unfolded FS under the single-q and triple-q CDW states, respectively, by setting f max ≡ max l,m,k {f l,m q1 (k)} = 88meV. In the single-q 1 case, sizable Fermi arc appears in FS1,3 due to the band-folding byf q1 . In the same way, Fermi arc appears in FS1,2 (FS2,3) byf q2 (f q3 ). The expected charge density modulation by |f tr q1 (r = 0)| ∼ 0.5 in Fig. S3 (f) is δn ∼ (f max |f tr q1 (r = 0)|)N (0) ∼ 0.02. Counter-intuitively, the size of the Fermi arc in the triple-q case in Fig. 5 (b) is much reduced, where we set f max = (88/ √ 2)meV because b 0 ≈ c 0 in the present study; see SM D [51]. To understand the reason, we consider the hybridization between Fermi momenta k, k+q 1 and k−q 2 in Fig. 5 (c). In the FS reconstruction by two form factorsf q1 andf q2 , the state |k hybridizes with |k + q 1 and |k − q 2 at the same time. Sincef q1 ∼f q2 , one eigenstate |k + q 1 − |k − q 2 is unhybridized and therefore gapless. (For general hybridization potentials, one of the three bands always remains unhybridized.) For this reason, after the unfolding, the Fermi arc structure around k + q 1 and k − q 2 in Fig. 5 (a) is recovered, as shown in Fig. 5 (b). Also, the spectral recovery in the unfolded bandstructure is explain in the SM B [51]. This hallmark in the triple-q CDW state could be observed by high-resolution ARPES study. Figure 5 (d) shows the obtained density-of-states (DOS). The pseudogap at E F in the triple-q CDW is small by reflecting the short Fermi arc in Fig. 5 (b). This result is consistent with experimental good metallic state below T IC .
In summary, we succeeded in explaining the triple-q IC-CDW in 1T-TaS 2 in terms of the "unconventional CDW", in which the form factor has strange orbitalmomentum-energy dependences. Thanks to the present paramagnon interference mechanism, the triple-q IC-CDW state is naturally understood based on a simple Hubbard model with on-site U , without introducing any nonlocal interactions. The same mechanism would be applicable for 1T-VSe 2 and other TMDs. Based on the knowledge on the IC-CDW state obtained by this study, it would be useful to develop Ginzburg-Landau theory to understand the NC-and C-CDW states.
We are grateful to R. Tazai for useful discussions. This study has been supported by Grants-in-Aid for Scientific Research from MEXT of Japan. This work was supported by the "Quantum Liquid Crystals" No. JP19H05825 KAKENHI on Innovative Areas from JSPS of Japan, and JSPS KAKENHI (JP18H01175, JP17K05543, JP20K03858).

[Supplementary Material]
Unconventional orbital-charge density wave mechanism in transition metal dichalcogenide 1T- TaS  Here, we construct the multiorbital Hubbard model for 1T-TaS 2 . In this metal compound, the TaS 2 -plane gives the almost perfect two-dimensional bandstructure. Figure S1 (a) shows the TaS 6 octahedron that is the building block of the TaS 2 -plane. The 5d orbital levels in TaS 6 octahedron are shown in Fig. S1 (b). When the octahedron is regular, the 5d orbitals split into t 2g orbitals (1-3) and e g orbitals (4,5), respectively. However, t 2g orbitals split into orbital 1 and orbitals 2,3 due to the finite trigonal distortion in real 1T-TaS 2 . Then, the wavefunction of orbitals 1-5 are given as orbital 1 : |d 3z 3 −r 2 (S1) orbital 2 : a|d x 2 −y 2 + b|d xz (S2) orbital 3 : a|d xy − b|d yz (S3) orbital 4 : a|d yz + b|d xy (S4) where a 2 + b 2 = 1 and a > b > 0. Now, we calculate the bandstructure of 1T-TaS 2 using the Wien2k software, by referring the experimental crystal structure [1]. Next, we construct the 11 orbital tightbinding model using the Wannier90 software, by taking account of the 4d-orbital wavefunction in Eqs. (S1)-(S5). The orthogonality of the 5d orbital is maintained by setting a = 0.847 (b = 0.532). The wavefunctions of 1-3 orbitals are shown in Fig. S1 (c). Figure S1 (d) exhibits the bandstructure of the obtained model, which reproduces the original bandstructure by WIEN2k satisfactorily. The conduction bands are mainly composed of 1-3 orbitals. Figure S1 (e) shows the weight of each d-orbital on the FS. We see that three d-orbitals are heavily entangled in many parts of the FS. Due to this fact, all the form factors f l,m q (k) with l, m = 1 ∼ 3 become large in the DW equation analysis. This fact is favorable for high temperature transition temperature T IC .
We note that the small hole-pockets around Γ point in Figs. S1 (d) and (e) are not observed by recent ARPES study [2], and they do not exist in the recent first principles study [3]. Therefore, in order to eliminate their artificial contributions, we multiply the Green function G(k) with the Heaviside step function Θ(|k| − k 0 ) with energy (eV) k 0 = 0.15π in the present numerical study in the main text.
Next, we explain the multiorbital Coulomb interaction. The matrix expression of the spin-channel Coulomb in-teraction is The matrix expression of the charge-channel Coulomb interaction is Here, U (U ′ ) is the intra-orbital (inter-orbital) Coulomb interaction, J is the Hund's coupling, and J ′ is the pair hopping term. In the main text, we assume the relations U = U ′ − 2J and J = J ′ , and set the constraint J/U = 0.10. The obtained results are not sensitive to the ratio J/U .

B: NUMERICAL RESULTS, UNFOLDED BAND SPECTRUM
The spin susceptibility in the RPA, χ s l,l ′ ;m,m ′ (q), is given byχ where the element of the irreducible susceptibility is χ 0 l,l ′ ;m,m ′ (q) = −T k G l,m (k + q)G m ′ ,l ′ (k). G l,m (k) is the electron Green function.
In the present model, χ s l,l ′ ;m,m ′ (q) becomes large for l, l ′ , m, m ′ = 1 ∼ 3. Figure S2 (a) shows the total spin susceptibility χ s tot (q) ≡ l,m χ s l,l;m,m (q). We also display some important elements of χ s l,l ′ ;m,m ′ (q) in Figs. S2 (b)-(e) Next, we explain the form factor obtained by the DW equation. The major form factor in the present study is f 1,2 q1 , which is shown in Fig. 3 (b) in the main text. It mainly originates from the AL process due to large value of χ s 1,1;1,1 (q) at q ∼ q 1 and that of χ s 1,2;1,2 (q) at q ∼ 0. In Figs. S3 (a)-(e), we show other form factors f l,m q1 for l, m = 1 ∼ 3. Thus, both diagonal and offdiagonal form factors with respect to orbitals 1 ∼ 3 take large values in the present numerical study. In Fig. S3 (f), the signal at r = 0 represents the presence of the local charge modulation.
Thus, all elements of f l,m q1 with l, m = 1 ∼ 3 are large and important. The main form factor f 1,2 q1 induces not only the orbital polarization, but also net charge modulation, because orbitals 1 and 2 are not degenerated. For this reason, net charge modulation in Fig. S3 (f) exists in the present unconventional CDW state. In the triple-q CDW state, both the FS and the bandstructure are folded intricately into the folded BZ. They can be unfolded into the original size BZ, which correspond to the experimental results by ARPES measurements. In the main text, we show the unfolded FSs in both the single-q and triple-q CDW states. In Figs. S4 (a) and (b), we show the unfolded bandstructure in the (a) single-q CDW state and (b) triple-q CDW state, respectively, along the cut A-B in Figs. 5 (a) and (b). (The green lines gives the original spectra.) In the single-q CDW state, the spectrum weight around the Fermi level disappears due to the band hybridization. This missing spectrum is partially recovered in the triple-q CDW state, as we explain in the main text. The recovery spectrum is stressed as blue dots in Fig. S4 (b). This hallmark in the triple-q CDW state may be observed by high resolution ARPES study.

C: KERNEL FUNCTION IN THE DW EQUATION
Here, we explain the kernel function in the DW equation, I l,l ′ ,m,m ′ q , studied in the main text. It is given by the Ward identity from the one-loop fluctuation-exchange self-energy. In multiorbital models, it is given as is the matrix expression of the bare multiorbital Coulomb interaction for channel b.
The first term of Eq. (S9) corresponds to the Maki-Thompson term, and the second and third terms give Aslamazov-Larkin terms.

D: GINZBURG-LANDAU EQUATION
Here, we construct a simple Ginzburg-Landau Free energy for the CDW order Φ(k) = (η 1 f q1 (k), η 2 f q2 (k), η 3 f q3 (k)), where η i is the real order parameter and f qi (k) is the normalized form factor.
From now on, we derive the fourth-order coefficients b 0 and c 0 by using the Green functions and form factors. Here, we introduce six nesting vectors q n (n = 1 ∼ 6) with the relation q n+3 = −q n . We also introduce the scalar order parameter η n for the form factor f qn with the relation η n+3 = η * n . Note that the relation f qn = f * qn+3