Renormalization of the Mott gap by lattice entropy: The case of 1T-TaS2

In many transition-metal oxides and dichalcogenides, the electronic and lattice degrees of freedom are strongly coupled, giving rise to remarkable phenomena, such as metal-insulator transition (MIT) and charge-density wave (CDW) order. We study this interplay by tracing the instant electronic structure under ab initio molecular dynamics. Applying this method to a 1T-TaS2 layer, we show that the CDW-triggered Mott gap undergoes a continuous reduction as the lattice temperature raises, despite a nearly constant CDW amplitude. Before the CDW order undergoes a sharp first-order transition around the room temperature, the dynamical CDW fluctuation already shrinks the Mott gap size by half. The gap size reduction is one order of magnitude larger than the lattice temperature variation. Our calculation not only provides an important clue to understand the thermodynamics behavior in 1T-TaS2, but also demonstrates a general approach to quantify the lattice entropy effect in MIT.

1T-TaS 2 has perhaps the richest electronic phase diagram of all transition-metal dichalcogenides because of the intertwined lattice, charge, orbital and spin degrees of freedom [1]. While the low-temperature commensurate charge density wave (CCDW) order and the accompanied Mott transition have been investigated for a long time by diffraction [2], transport [3], scanning tunneling microscopy [4][5][6] and angle-resolved photoemission [7][8][9], the absence of magnetic susceptibility of the Mott phase remains puzzling [10]. The possibility of a quantum spin liquid state due to the lattice frustration was proposed recently [11], which aroused revived theoretical interest and stimulated a series of recent experiments [12][13][14][15][16].
The general consensus [17] is that below 200 ± 20 K (T C ), the √ 13 × √ 13 CCDW order is fully established. The Ta atoms are grouped into 13-atoms clusters with a "Star-of-David" (SD) arrangement. It is widely perceived that this phase can be viewed as a cluster Mott insulator -each SD acts effectively as a correlated site with an odd number of electrons and the SDs form a layered triangular lattice [11]. Above T C to around 350 K (T N C ), the so-called nearly CCDW state is phase separated, where the CCDW domains are separated by incommensurate areas. Above T N C , the SD clusters completely disappear, leaving a weak incommensurate CDW order.
The controversy mainly focuses on the spin dynamics of the CCDW phase. In contrast to a band insulator, the spin degree of freedom in a Mott insulator is not inert, which will manifest in thermodynamic measurements. The experimental results are quite unusual. Although it has been long noticed that the low-T magnetic susceptibility displays a Curie-like tail, the magnitude is much smaller than one would expect from one free elec-tron spin per SD [3]. Recent muon spin rotation measurement observed no spin ordering signal down to 70 mK [12]. The new specific heat data below 2 K clearly displays a magnetic-field-dependent linear-T ingredient [14]. Nuclear quadruple resonance measurement showed an anomalous crossover of the 181 Ta nuclear spin relaxation rate around 50 K [12].
These interesting experimental observations jointly point to a highly exotic spin state. However, a unified theoretical explanation is still elusive. Our previous study formulated a multi-orbital Anderson-type Hamiltonian based on the characterization of two-types of molecular orbitals [18]. It was also shown that the Hamiltonian parameters sensitively depended on the CCDW order, and the latter could be tuned by various mechanisms, such as pressure [1] and doping [19][20][21][22][23]. We note that it is not known a priori whether this multiorbital Hamiltonian has significant temperature dependence. This question has important implications on the effective spin model, which in principle is obtained by projecting out the lattice, charge and orbital degrees of freedom, and thus will inherit the temperature dependence.
The CCDW order, i.e. the SD superstructure amplitude φ SD (see definition below), as the primary order paramter is essentially a periodic lattice distortion, which can be nicely described by ab initio molecular dynamics (MD). Within the Born-Oppenheimer (BO) approximation, the ensemble-averaged electronic properties can be computed as a MD time average.
Our calculations are performed using the Vienna ab initio Simulation Package (VASP) [24][25][26][27]. Nucluei are subject to the Newton's equation of motion on the BO potential surface using a time step of 2 fs. To simulate a canonical ensemble with constant T, the Nose thermostat is used [28][29][30]. For each T, the MD simulation lasts for 20 ps, and the last 4 ps is used to calculate thermodynamic properties. We note that the nuclear quantum effects are not taken into account in the present calculation, which in principle can be evaluated by more sophisticated path-integral MD [31].
At any instant of time, the ground state of electrons is calculated self-consistently within DFT. We employ the projector augmented wave method [32] and the exchangecorrelation functional due to Perdew, Burke and Ernzerhof [33]. The +U correction is employed to capture the Coulomb interaction of Ta 5d orbitals on the Hartree-Fock level, following the simplified (rotational invariant) approach introduced by Dudarev [34]. We employ an effective U = 2.27 eV as previously derived from the linearresponse [35] calculation, which was found to nicely reproduce the experimental Mott gap. We use a plane-wave cutoff of 300 eV, and the Brillouin zone was sampled with the Γ point only.
The simulation cell contains a single layer of 52 Ta atoms sandwiched by 104 S atoms (in total N atom = 156), which can accommodate up to 4 SDs (Fig. 1). A 15Å vaccuum layer is included in the z-direction. The complexities of the stacking order of the single layer and the interlayer coupling are beyond the scope of the current calculation.
Numerically, it is important to guarantee that the last 4 ps has already achieved thermal equilibrium. For a better convergence, we start from low T, which is closest to the DFT relaxed static structure. Then, the structure of the final MD step is used as the initial structure of the next temperature, which is elevated progressively. When all the SDs melt, we reversely reduce the temperature progressively, using the equilibrium structure at the higher temperature as input. Finally, the simulation forms a complete heating-cooling cycle. Our criterion for thermal equilibrium are that (i) the temperature fluctuation has already converged to 2/(3N atom ) = 6.5% as expected from a Boltzmann distribution, (ii) clear periodicity with constant amplitude can be observed from the atomic displacement, and (iii) the observables from the heating and cooling processes coincide when the temperature is away from the transition point.
We define the CCDW order parameter as φ SD = d inter −d intra , whered inter (d intra ) is the time-averaged Ta-Ta bond length between SDs (within a SD). The definition of inter-(intra-)SD bonds is ambiguous in the high-T phase, so we always refer to the SD positions in the CCDW phase. Bothd inter andd intra are calculated from the time-averaged lattice structure as a function of T. In general when the SDs melt, it is expected that φ SD vanishes. Figure 2(a) plots φ SD versus temperature. A sharp first-order transition can be observed. The loop shows hysteresis, indicating that there are two competing phases and the atoms have not fully relaxed to the low free-energy phase within the simulation duration. The transition temperature T C is 250 K ∼ 300 K.
It is usually challenging to reliably simulate the firstorder phase transition, because the equilibrium time might be too long for MD. The good agreement obtained here should be attributed to a relatively shallow energy barrier between the two phases. Given that the simulation is performed for a 2D layer, one natural question is how to reconcile with the Mermin-Wagner theorem. The answer lies in the imposed periodic boundary condition, which cuts off any thermal fluctuation with a wave vector larger than the cell vectors. On the other hand, the most important thermal fluctuation that melts the SDs has the √ 13 × √ 13 wave vector, which is fully accommodated in our simulation cell. The transition occurs when the kinetic energy of the atoms becomes large enough to escape the √ 13 × √ 13 potential well. This physics is faithfully characterized by MD. Therefore, it is not surprising that our 2D simulation is able to predict a finite transition temperature.
We should also note that our simulation cell is still not large enough to describe phase separation and longwave incommensurate CDW. Nevertheless, in Fig. 1(b), some inhomogeneity can already be seen in the 300 K structure. Experimentally, between T C (the SDs start to melt) and T N C (the SDs completely melt), there is a wide range in which the CCDW domains and the discommensurate regions coexist [17]. It is understandable that our simulated transition temperature falls between the experimental T C and T N C .
Next, the finite-T electronic properties are calculated within the static approximation and the BO approximation respectively. In the former treatment, the time average is performed for the lattice structure first, and then the band structure is calculated based on the timeaveraged lattice, from which we read the Mott gap E static g . In the latter treatment, an instant electronic structure is calculated at each MD step. The instant gap size is time averaged afterwards to obtain E BO g . Within the static approximation, the temperature dependence solely arises from the variation of the CCDW order parameter φ SD (T ). Figure 2(b) plots E static g versus T. A sharp metal-to-insulator transition takes place in accompany with the CCDW transition. Figure 3(a-d) show the band structure at four typical temperatures.
A detailed comparison between the DFT+U results and STM data at 5 K was carefully made in our previous work [18]. The two sets of narrow bands centered around ±0.2 eV was assigned as the upper Hubbard (UH) and the lower Hubbard (LH) bands, respectively, which arise from the central Ta atoms in the SDs, and are separated apart from the dispersive conduction bands (CB) and valence bands (VB) consisting of the surrounding Ta atoms by a small CDW gap (∆ CDW ). Note that our MD cell is four times larger than the √ 13 × √ 13 supercell used previously, and thus all the bands are folded.
The static bands only show marginal temperature dependence below T C . If we carefully examine Figs. 3(ac), we can observe a slight reduction of E static g , as well as an increasingly overlap between UH (LH) and CB (VB), which is consistent with a small decrease of φ SD . Above T C , the Mott gap collapses abruptly [ Fig. 3(d)].
Going beyond the static approximation, Figs. 3(e-h) plot the instant energy levels at the Γ point at each MD step. At 5 K, the gap size agrees well with the static band structure, despite slight temporal fluctuations. At higher temperatures, the temporal fluctuations significantly reduces the static gap size. At 275 K, instant gap closing can be observed, indicating that the system is close to the phase transition. The time-averaged E BO g as a function of temperature is plotted in Fig. 2(c). The gap size decreases from around 0.4 eV at 5 K to around 0.2 eV at T C .
The comparison between E static g and E BO g clearly indicates that the driving force of the Mott gap reduction is not the static φ SD order but the vibrations around the SD superstructure. In other words, the Mott gap is strongly renormalized by the electron-phonon coupling (EPC). In conventional semiconductors, E g typically decreases with increasing temperature at a rate ∼ 0.1 meV/K, i.e. ∆E g /k B ∆T ∼ −1, as a consequence of lattice dilatation and EPC [36]. Here, ∆E BO g /k B ∆T is of the order of −10. The absolute gap size shrinks by half when the temperature approaches T C .
This giant temperature dependence is closely related to the correlated nature of the Mott gap. It is important to notice that the lattice vibration not only perturbs the lattice potential, as commonly considered in the single-electron EPC calculation, but also strongly modifies the effective Coulomb repulsion experienced by the Hubbard electrons. An interesting parallel can be drawn with respect to Se-substituted 1T-TaS 2 , by viewing Sesubstitution induced distortions as some "frozen phonon" modes. Our previous work reported similar gap reduction [18]. The underlying physics was explained by a multi-orbital Hubbard model using two types of molecular orbitals as the basis. It was demonstrated that the energy spacing between these two molecular orbitals could be sensitively tuned by Se substitution, which in turn modified the effective Coulomb repulsion experienced by the Hubbard electrons. Consequently, the gap size varies significantly.
In summary, the key information revealed by the MD simulation is that the Mott gap of 1T-TaS 2 in the CCDW phase is NOT a constant. Considering that the spin superexchange magnitude is correlated with the Mott gap, the gap size change will play an important role in correctly interpreting the unusual spin dynamics as observed in the specific heat and nuclear quadruple resonance measurements. In principle, a temperature-dependent effective spin model should be formulated. We noticed that some signatures of gap size reduction could indeed be observed by comparing the recent scanning tunneling spectroscopy data measured at different temperatures [18,37,38]. A previous angle-resolved photoemission measurement also reported the temperature dependence of the photoelectron spectra, despite a different interpretation in terms of pseudogap [39].
One remaining question is whether we can quantify the temperature dependence of the spin exchange parameters based on the MD data. For conventional magnets, the Heisenberg exchange is typically fitted by comparing the energy difference of various static spin configurations. For 1T-TaS 2 , we find that all the spin configurations we can obtain in the MD cell result in nearly degenerate energy. A trivial interpretation is that the in-plane spin interaction is indeed weak. Then, one has to resort to interlayer coupling in order to explain the missing of the magnetic moment. A nontrivial interpretation is that the degeneracy arises from frustration, and the effective spin model may contain additional terms beyond the Heisenberg exchange. One possible solution is to derive the spin interactions from the multi-orbital Hubbard model by perturbation calculation. We will leave this for future investigation.