Orbital selective Mott transition and magnetic moment in charge density wave heterostructures NbSe${}_{2}/$Ta$X_{2}$

Abstract

We investigate the electronic properties of charge density wave (CDW) heterostructures out of monolayers of 1T-NbSe₂ and 1T-Ta$X_{2}$ (where, $X=$ S and Se) using first-principles followed by dynamical calculations. The CDW-ordered crystal structures are simulated using $\sqrt{13}\times\sqrt{13}$ supercells of NbSe₂ and Ta$X_{2}$. These two-dimensional heterostructures are modeled by stacking monolayers of NbSe₂ and Ta$X_{2}$ along (001) direction. Our investigations reveal the presence of non-zero magnetic moments in NbSe${}_{2}/$TaS₂, albeit without a long-range magnetic order, raising the issue of a possible quantum spin liquid (QSL) as suggested for monolayer 1T-TaS₂ recently. In contrast, the NbSe${}_{2}/$TaSe₂ heterostructure exhibits no magnetic moment. In order to capture the dynamical effects of local correlation, we use DFT plus multi-orbital dynamical mean field theory (MO-DMFT). Our findings indicate that NbSe${}_{2}/$TaS₂ is considerably influenced by the dynamic corrections, whereas NbSe${}_{2}/$TaSe₂ shows minimal effects. Additionally, an orbital-selective Mott transition (OSMT) is observed in the NbSe${}_{2}/$TaS₂ bilayer heterostructure.

I Introduction

Flat bands in electronic structures indicate a significant localization of electrons, where the on-site Coulomb interaction ($U$) dominates over the bandwidth ($t$), i.e., $U/t>1$, leading to a Mott transition due to electron correlations. Materials such as NbSe₂, TaS₂, and TaSe₂ exhibit metallic, dispersive bands in their pristine bulk and monolayer forms. Notably, the restructured $\sqrt{13}\times\sqrt{13}$ star-of-David (SoD) clusters in 1T-TaS₂ and 1T-TaSe₂ host these flat bands [1, 2, 3, 4]. Monolayer SoD structures in 1T-NbSe₂ are also studied for their flat band and Mott insulation properties. ARPES measurements confirm the presence of a robust flat band persisting from 40K up to 450K [5]. For 1T-TaSe₂ and 1T-NbSe₂, ARPES studies reveal Mott half-gaps of 0.28 eV and 0.23 eV, respectively [5]. 1T-TaS₂, recognized as a Mott insulator, exhibits Mott gaps below the Fermi level ranging from 0.2 to 0.4 eV in photoelectron spectroscopy and ARPES experiments [6, 7, 8, 9, 10]. Theoretical investigations using density functional theory (DFT) and dynamical mean field theory (DMFT) calculations consistently support these experimental findings [7, 8, 10, 11, 12, 13, 14, 15].

In its bulk form, NbSe₂ crystallizes in a hexagonal structure and undergoes a charge density wave (CDW) transition at low temperatures [16, 17, 18, 19, 20]. This CDW phase transition results in a periodic modulation of the charge density, leading to an incommensurate CDW state at higher temperatures and a commensurate CDW state at lower temperatures. When reduced to a monolayer, NbSe₂ exhibits intriguing changes in its electronic characteristics due to quantum confinement. Superconductivity in NbSe₂ samples is widely investigated [21, 22, 23] and owing to a large spin-orbit coupling (SOC) Ising superconductivity is being looked into, recently [24, 25, 26]. The competition between CDW and superconductivity is well known in the literature [27, 28, 29, 30].

Likewise, TaSe₂ and TaS₂ have hexagonal lattice structures in their bulk from and are known for their CDW phase transitions at low temperatures and high pressures [31, 32, 33, 34]. The Fermi surface nesting has long been rationalized a consequence for the CDW instability in transition metal dichalcogenides (TMDs) [35, 36, 37]. However, the mechanism has been disregarded in experiments [38, 39, 40, 41, 42]. A correlated, preformed exciton liquid route is then proposed by Taraphder et al., using local density approximation plus MO-DMFT, in which the CDW state arises as Bose condensation of preformed excitons [43, 44].

In addition to the Mott insulator, 1T-TaS₂ in its CDW phase also exhibits quantum spin liquid (QSL) [45] behavior caused by the lack of long range magnetic ordering [46]. The spin liquid state was initially proposed by Anderson and Fazekas for triangular lattice due to magnetic frustration [47, 48]. The half-filled, localized band centred on the central Ta site of the star of David in the CDW state is likely to have a moment. Recently, the monolayer CDW phase of 1T-NbSe₂, has been predicted to be a QSL state, when placed on a $3\times 3$ modulated 1H-NbSe₂ monolayer [49]. The scanning tunneling spectra (STS) of the 1T/1H-NbSe₂ heterostructure (HS) shows a sharp Kondo peak at the Fermi level indicating a considerable interaction between spin and charge states in two layers. A similar Kondo peak in STS was also observed earlier for 1T/1H-TaSe₂ HS [50].

Figure 1: Crystal structures, (a-d) H1: 1T-NbSe₂/TaS₂ and (e-h) H2: 1T-NbSe₂/TaSe₂. (a,e) Top and (b,f) side views. (c,d,g,h) The formation of SoD are shown. The bond lengths are mentioned therein.

The discussion above suggests that the triangular lattices of TMDs are good candidates for the QSL state. By combining different 2D materials with monolayer NbSe₂, TaS₂ or TaSe₂ in HSs, one can create and modulate the quantum states. In the current investigation, we model 2D CDW HSs using 1T polymorph of NbSe₂ and Ta$X_{2}$, where $X=$ Se and S. Our first-principles DFT calculations suggest a possible QSL state in 1T-NbSe₂/TaS₂ due to the absence of long range magnetic ordering. The magnetic moments on all the Ta/Nb transition metals (TMs) are not uniquely distributed. On the contrary, the QSL state is not observed for 1T-NbSe₂/TaSe₂ experimentally. We also find no magnetic moment for this heterostructure. Moreover, an orbital selective Mott transition is noted for both the HSs, by incorporating dynamical correlations within DFT+DMFT approach.

II Computational Details and Methodology

II.1 DFT

The electronic structures have been calculated using a plane wave implementation of DFT within the Vienna ab-initio simulation package (VASP) [51, 52, 53, 54] which uses projected augmented wave (PAW) [55, 56] potentials. The generalized gradient approximation [57, 58] is used for the exchange-correlation functional, and dispersion corrections are included within the DFT-D2 approach [59]. Both the lattice parameters and the internal positions are optimized until an energy convergence of $10^{-4}$ eV and force convergence of $10^{-3}$ eV/ Å have been achieved. $\Gamma$-centered k-mesh of 7$\times$7$\times$3 and 8$\times$8$\times$1 has been used for bulk and 2D CDW phases, respectively. A cutoff energy of 520 eV has been used for the plane-wave basis.

II.2 DMFT

We perform the DMFT calculations in a basis set of projective Wannier functions, calculated using DFTTOOLS and soliDMFT packages [60, 61, 62, 63] interfaced with TRIQS libraries [64]. The Anderson impurity problems are solved using continuous-time quantum Monte Carlo algorithm in the hybridization (CTHYB) expansion [65] as implemented in TRIQS [66]. We use a fully rotationally invariant Kanamori Hamiltonian parametrized by Hubbard’s $U$ and Hund’s coupling $J_{H}$, where we set the intraorbital interaction to $U^{\prime}=U-2J_{H}$ [67]. For our calculations we used various $U$ values ranging from 1.0 to 3.0 eV and $J_{H}=0.8$ eV in order to see the effect of correlations. Real-frequency spectra are obtained using the maximum-entropy method of analytic continuation as implemented in the TRIQS/MAXENT application [68]. We carry out all the DMFT calculations at an inverse temperature $\beta=40\,eV^{-1}$.

Figure 2: Band structures without spin-polarization of (a-c) 1T-NbSe₂/TaS₂ and (d-f) 1T-NbSe₂/TaSe₂ with (a,d) GGA, (b,e) GGA$+U$, where an average value of $U=2.61$ eV is used, and (c,f) GGA$+U$, where $U_{\rm{Nb}}=2.95$ eV and $U_{\rm{Ta}}=2.27$ eV. The red bands are used for Wannier projections.

III Crystal Structures

The crystal structures are generated using 1T polymorphs of TMDs, $MX_{2}$, in which the transition metals (TMs), $M$ are sandwiched between the chalcogens, $X$, in octahedral symmetry. The $\sqrt{13}\times\sqrt{13}$ supercells are created using 13 TM atoms cluster in which 12 TM atoms surround a central TM atom. The surrounding TM atoms are deviated slightly from their original Wyckoff positions and we notice that after relaxation they form SoD like structure. These are well known CDW crystal phases. Next, the HSs are formed by deposition of single layer of SoD structure above another with an interlayer distance, $d_{M-M}\sim 6$Å. We choose 1T-NbSe₂/TaS₂ and 1T-NbSe₂/TaSe₂ HSs for our investigation which we refer H1 and H2, respectively, throughout the paper. For the 2D bilayer HSs a vacuum of 15Å  has been added along third direction in order to discard spurious interactions between the periodic unit cells. The crystal structures, interlayer separations between the nearest chlcogens of two layers and the bond lengths are shown in Fig. 1.

Table 1: Energy per formula unit cell (in eV) for H1 and H2.

Structural phase	for H1 (eV)	for H2 (eV)
1T Normal	-22.30	-21.53
1T CDW	-22.32	-21.48

In order to check the structural stability of the CDW phase in both the HSs, we calculate energy per formula unit cell ($E_{fu}$) in the normal phase and in the CDW phase. The energies values are provided in Table 1. $E_{fu}$ is evaluated by calculating the DFT ground state energies of the bulk (containing two $MX_{2}$ units) and the bilayer CDW structures (containing 26 $MX_{2}$ units), then dividing the totat energy by the number of $MX_{2}$ units in that particular system. We find that after the formation of SoD clusters, $E_{fu}$ for H1 becomes slightly more negative than $E_{fu}$ of the normal 1T bulk phase. Conversely, $E_{fu}$ becomes less negative in the case of H2. Therefore, it turns out that the H1 is more stable in its CDW phase while H2 destabilizes in the cluster form. The destabilization of H2 can be attributed to the large bond length differences from the outer Ta atoms to the central Ta atom of SoD as shown in Fig. 1(h), which deforms the SoD. The difference in the bond lengths is smaller for H1, unlike H2 (see Fig. 1(d)). Additionally, we calculated $E_{fu}$ for the monolayer CDW phases of 1T-NbSe₂, 1T-TaS₂ and 1T-TaSe₂, and the measured values are $-19.59$ eV, $-21.18$ eV, and $-21.25$ eV, respectively. Though the direct comparision of the energies has no physical meaning the average energy of the two layers is still higher than the bilayer system, indicating the bilayers HSs are more stable than the monolayers. The Van der Waals interactions could be playing a role in stabilizing the bilayer systems.

IV Results and Discussions

IV.1 DFT, DFT$\mathbf{+\textit{U}}$ and spin-polarization

We calculate the band structures for H1 (Fig. 2(a-c)) and H2 (Fig. 2(d-f)) under GGA and GGA+$U$ formalisms. We find that without on-site Hubbard $U$, the highest valence band (HVB) crosses the Fermi level (FL) along $\Gamma-L$, and the lowest conduction band (LCB) merely touch the FL in case of H1 (see Fig. 2(a)). In case of H2, the HVB and LCB do not cross but touch the FL as shown in Fig. 2(d). However, the HVB and LCB cross or touch the FL along differnt momentum directions in H1 and H2. For H1, the HVB cross the FL along $\Gamma-M$ and $K-\Gamma$, whereas the momentum direction is $M-K$ for H2, where the bands touch the FL. The HVB in H1 forms tiny hole pockets along $\Gamma-M$ and $K-\Gamma$ Brillouin zone (BZ) directions. In rest of the BZ, the bands are nearly flat, indicating strong correlations supported by the large on-site interactions and small kinetic energy between the electrons. The band widths of HVBs for H1 and H2 are $0.12$ and $0.09$ eV, respectively. Our atom and orbital-projected density of states (DOS) calculations suggests that the four bands near the FL (highlighted in red in Fig. 2) corresponds to the central Nb (Ta) $4d$ ($5d$) orbitals. Moreover, there is a strong hybridization of different $d$ orbitals.

In order to see the effect of static on-site Hubbard $U$ on electron localization, we perform GGA$+U$ calculations, with different $U$ values for those having $d$ orbitals, Ta and Nb TM atoms in our case. As shown in the band structure for H1, $U_{\rm{Nb/Ta}}=2.61$ eV has no significant effect on the overall band topology, except, the two valence bands close to the FL are shifted down by a negligible amount of $0.012$ eV. On the other hand if we use different values of $U$ for Ta ($U=2.95$ eV) and Nb ($U=2.27$ eV), the two conduction bands slightly move up in the energy with respect to the GGA bands, and the two valence bands shift down, only touching the FL, without crossing it. The hole pockets in this case disappear. The band gap values achieved for H1 with $U_{\rm{Nb/Ta}}=2.61$ and $U_{\rm{Nb\,(Ta)}}=2.95\,(2.27)$ eV are $9.84$ meV and $10.5$ meV, respectively. It is important to emphasize here that the monolayer 1T CDW phase of NbSe₂ [13, 69] and TaS/Se₂ [5] are considered to be Mott insulators, with $U=2.8$ eV. The Mott insulation phase is also achieved in the case of HSs but rather small band gap values. Though we have not checked, following the discussion above regarding varying $U$ and shifting bands, we believe a band gap can be induced if we increase $U$ for Ta. On the other hand, a band gap of $14.89$ meV is opened with $U_{\rm{Nb/Ta}}=2.61$ eV in the case of H2 (see Fig. 2(e)). Similar to H1, the LCB shifts up with $U_{\rm{Nb}}=2.95$ eV and $U_{\rm{Ta}}=2.27$ eV (see Fig. 2(f)).

Figure 3: The magnitudes of magnetic moments in respective layers of stars of H1 are shown by the lengths of arrows with $U_{\rm{Nb}}=2.95$ eV and $U_{\rm{Ta}}=2.27$ eV.

As discussed in the introduction, monolayer 1T-TaS₂ exhibit a QSL state, which regards a lack of long range magnetic ordering. In the following, we discuss the magnetic structure of the bilayer HSs. We perform spin-polarized DFT and DFT$+U$ calculations, and find that the 1T CDW phase comes out to be non-magnetic without $U$ for both the HSs. Incorporating $U$ in the calculations induce magnetic moment on TM atoms for the H1 HS, and most of the contribution of magnetic moment is concentrated on the central Ta and Nb atoms, as denoted by arrows in Fig. 3. The size of the arrow depicts the magnitude of the magnetic moment, which is maximum for the central atoms of SoD. The magnetic moment of the central Nb and Ta, with $U_{\rm{Nb}}=2.95$ eV and $U_{\rm{Ta}}=2.27$ eV, are $0.684$ and $0.311$ $\mu_{B}$, respectively, where, $\mu_{B}$ is Bohr magneton. It is interesting to see that in the SoD of NbSe₂ layer, the nearest Nb atoms from the central Nb show the secondary contribution in the magnetic moment while for the TaS₂ layer, the secondary contribution is dominated from the next nearest neighbors. The secondary contribution in the magnetic moments is of the order of $10^{-2}$ $\mu_{B}$, while the lowest contribution is either zero or of the order of $10^{-3}$ $\mu_{B}$. The magnetic moment on each species is not uniquely defined, causing the absence of long range magnetic ordering in a particular layer as well as in the bilayer, which might result into a QSL state. We perform spin-polarized calculations for $U_{\rm{Nb/Ta}}=2.61$ eV, and found that only central atoms dominantly contribute to the magnetic moment with the values, $0.117$ $\mu_{B}$ and $-0.064$ $\mu_{B}$ for Nb and Ta, respectively. The rest of the TMs exhibit the magnetic moment of the order of $10^{-3}$ $\mu_{B}$ or zero. Here, the negative value of moment for Ta represents spin down. Furthermore, the spin-polarized band structures clearly show spin-split bands with non-zero $U$ (not shown). Additionally, the spin-polarized bands of H1 with $U_{\rm{Nb/Ta}}=2.61$ eV show a band gap of $8.34$ meV, while the system is metallic with $U=0$ and $U_{\rm{Nb\,(Ta)}}=2.95\,(2.27)$ eV. H2 HS, on the contrary, is non-magnetic with and without $U$, and there is no spin-plitting in the spin-polarized bands. As reported earlier, and mentioned in the introduction, 1T-TaSe₂ shows the signature of QSL with 1T-TaS₂ [70, 50], 1T-TaS₂ itself being a QSL [46]. Evidently, our investigation finds that the HS containing 1T-TaS₂ (H1) could be a potential candidate for QSL, unlike the H2.

IV.2 DMFT calculations

In order to see the electron dynamics we perform DMFT calculations. It is evident from the DFT band structure calculations that H1 and H2 exhibit flat bands near the FL, postulating strong correlations effects in the electrons. Strong correlation can induce Mott insulating behavior in 2D limit. We extract four bands around the FL by adopting the maximally localized Wannier functions (MLWFs) [71] using Wannier90 package [72] as implemented in QUANTUM ESPRESSO (QE) [73, 74, 75]. As discussed above in Sec. IV.1, there are multiple 4$d$ and 5$d$ orbitals hybridizing near the FL. Our extracted Wannier bands are nearly good match with DFT bands, and also the spreads of Wannier functions are small. The slight mismatch can be attributed to the orbital-hybridization, and considering only four most contributing orbitals near the FL. We choose $4d_{xy}$, $4d_{z^{2}}$, $5d_{xy}$ and $5d_{z^{2}}$ for Wannier fitting. Different values of Hubbard parameter, $U$, in the fully rotationally invariant Kanamori Hamiltonian [67] are used to see the variations in the bands and to investigate the Mott insulating phase. We would like to mention that there is a little mismatch in the Fermi energy calculated using VASP and QE. Owing to that, the valence bands extracted with Wannier90, shifts up in the energy and cross the FL, which we are going to discuss shortly. Also, $U$ is not taken into account in the DFT bands from which the Wannier bands are extracted.

Figure 4: Upper panel: spectral functions with different values of $U$, and bottom panel: the corresponding Fermi surfaces for (a-d) H1 and (e-h) H2. The dashed line denotes the FL.

The momentum resolved spectral function (DMFT spectrum) is shown for H1 in Fig. 4(a-d) and for H2 in Fig. 4(e-h). The DFT bands are shown in black while the bright bands are spectral functions. As it can be seen that for H1 the lowest two conduction bands which are called upper Hubbard bands (UHBs) in the literature move significantly up with $U=1$ eV (see Fig. 4(a)) and separate from the highest two valence bands which are called lower Hubbard bands (LHBs). With increasing $U$ the two UHBs goes up in the energy very rapidly from $U=1$ to $U=3$ eV. Notably, the LHBs also start isolating from each other at $U=2$ eV (see Fig. 4(b)) and at $U=2.61$ eV (see Fig. 4(c)) they completely detach and open up a gap in the spectrum. This gap opening is a signature of orbital-selective Mott transition (OSMT). OSMT is defined as a phase in which some of the orbitals show strong correlations (localized electrons) while the others remain itinerant [76, 77]. The OSMT is robust with higher $U$ values (see Fig. 4(d)). In addition to OSMT, splitting in the DMFT bands can be seen at the $\Gamma$ point, unlike the DFT bands. A weak satellite feature is also apparent near the $\Gamma$ point with $U=1$ and $2$ eV in the lower VHB. The spectral weight is isolated from the original band around $\Gamma$ point opening a gap. The satellite features in the spectral function are previously reported for cubic SrVO₃ [78]. These satellites are signature of weak localization.

The OSMT is also evident from the Fermi surface (FS) shown in the bottom panels of Fig. 4(d). As presented in Figs. 4(a-d), the distribution of spectral function is denoted by blue and yellow intensities. The intensities start disappearing towards higher $U$ values and almost vanish for $U=3$ eV, except, a smaller and intense hole pocket around the $\Gamma$ point. The FSs further support OSMT in H1. On the other hand, H2 neither exhibit OSMT nor the satellites. As apparent from Figs. 4(e-h) the HVB (or LHB) always intersect the FL. Moreover, a shifting in UHBs to the higher energy is visible in H2 also with $U$, but the sensitivity with $U$ is smaller as compared to the H1 HS. The LHBs, unlike UHBs, are almost unaltered with $U$ and show no shifting.

V Conclusions

In conclusion, the CDW bilayer phases of 1T-NbSe₂/TaS₂ (H1) and 1T-NbSe₂/TaSe₂ (H2) are investigated for their structural and electronic properties. It is found that after relaxation H1 form a SoD structure, while the SoD is deformed in case of H2, leading to a unstable CDW structure of H2 as compared to its bulk phase. The ground state energy per formula unit ($E_{fu}$) suggests that H1 is stable in its CDW phase with a $0.02$ eV lesser $E_{fu}$. Also, the bilayer HSs are found to be stable than their monolayer counterparts. The band structure calculations presents flat bands in the spectrum indicating strong correlations due to electron localization. We observe metallic (H1) and nearly insulating (H2) states without $U$. With non-zero $U$ band gaps is produced in both the heterostructures resulting into Mott insulators. We also report a possible QSL behavior in H1 which we explain by non-zero and random magnetic moments on TM atoms. A QSL state is previously observed experimentally in monolayer CDW phase of 1T-TaS₂ and 1T-TaS₂-based bilayers systems. Unlike H1, in H2, we do not find magnetic moments on TMs with or without $U$, indicating QSL is uniquely defined for 1T-TaS₂-based HSs. Moreover, our DMFT calculations predict OSMT in H1, and also satellite features in the spectral function. The OSMT and satellites are absent in H2 HS. We find a significant sensitivity to $U$ in the spectral function for H1, unlike H2 that shows a very slow variation by varying $U$. The LHBs in H2 are almost fixed and hardly show any variation at $U\geq 2$ eV.

VI Acknowledgements

AT and JC acknowledge funding from SERB (Government of India) CRG project No: CRG/2021/001078 under which this work was performed. SP thanks Alexander Hampel for technical help and useful discussions regarding TRIQS software. The authors acknowledge National Supercomputing Mission (NSM) for providing computing resources of ’PARAM Shakti’ at IIT Kharagpur, which is implemented by C-DAC and supported by the Ministry of Electronics and Information Technology (MeitY) and Department of Science and Technology (DST), Government of India.