Photoinduced Dynamics of Commensurate Charge Density Wave in 1T-TaS$_{2}$ Based on Three-Orbital Hubbard Model

We study the coupled charge-lattice dynamics in the commensurate charge density wave (CDW) phase of the layered compound 1T-TaS$_{2}$ driven by an ultrashort laser pulse. For describing its electronic structure, we employ a tight-binding model of previous studies including the effects of lattice distortion associated with the CDW order. We further add on-site Coulomb interactions and reproduce an energy gap at the Fermi level within a mean-field analysis. On the basis of coupled equations of motion for electrons and the lattice distortion, we numerically study their dynamics driven by an ultrashort laser pulse. We find that the CDW order decreases and even disappears during the laser irradiation while the lattice distortion is almost frozen. We also find that the lattice motion sets in on a longer time scale and causes a further decrease in the CDW order even after the laser irradiation.


Introduction
Transition-metal dichalcogenides (TMDs) are layered compounds that are a testbed to study electron correlations and electron-lattice interactions in two dimensions [1,2]. In particular, 1T-TaS 2 shows a variety of phases including charge density waves (CDWs) accompanied by periodic lattice distortions (PLDs) at low temperatures [3][4][5][6][7]. In the lowest-temperature phase below 180 K, electrons experience a lock-in to the commensurate CDW (CCDW) that forms a √ 13 × √ 13 hexagram structure (see Figure 1), which has been also referred to as the Star of David. The CCDW phase is insulating as revealed by early resistivity and susceptibility measurements [1], and this is confirmed by recent experiments of angle-resolved photo-emission spectroscopy (ARPES) [8]. This is believed to be a Mott insulator [9] (see also [10,11] for alternative scenarios), and the possibility of spin liquid is under debate [12,13].
Theoretical studies for these intriguing nonequilibrium phenomena are quite limited except phenomenological Ginzburg-Landau theories [30][31][32]. This is partly because modeling the electronic structure has not been settled yet for the CCDW phase in 1T-TaS 2 . Constructing an empirical tight-binding model dates back to Mattheiss [33], who proposed a three-orbital Hamiltonian for describing the normal state without PLD. Smith et. al. [34] studied the case in the presence of the PLD and discussed the band structure well below the Fermi energy. Rossnagel and Smith [35] included the spin-orbit interaction and showed that this isolates one narrow band near the Fermi energy. They conjectured the Mott insulating phase of this band, but it remains open to examine its validity. arXiv:1812.00297v1 [cond-mat.str-el] 2 Dec 2018 In this paper, we extend the empirical tight-binding model of Rossnagel and Smith [35] by adding the on-site Coulomb interaction, and obtain a band gap at the Fermi level within a mean-field approach. The gap formation is in line with their conjecture [35] and consistent with the electronic property of the CCDW phase of 1T-TaS 2 . We then treat the PLD of the hexagram shape as a dynamical variable and formulate the coupled equations of motion for the electrons and the PLD. On the basis of these equations, we study the electron-lattice cooperative dynamics induced by an ultrashort laser pulse within a time-dependent mean-field approximation. On a short time scale, the lattice distortion is almost frozen, but the laser input induces coherent electron dynamics and the charge density wave decreases and even disappears. On a longer time scale, the lattice motion sets in and causes further decrease in the charge-density-wave order even after the laser irradiation.

Three-Orbital Tight-Binding Model: Previous Studies
We first introduce a model for the phase without PLD, and this is the tight-binding Hamiltonian used in [35]. This model involves three d-orbitals (3z 2 − r 2 , x 2 − y 2 , and xy) of Ta ions on the two-dimensional triangular lattice, and is defined aŝ where the summation ∑ i,j runs over the nearest neighbor site pairs. Hereĉ imσ denotes the electron annihilation operator at site i, in orbital m (= 3z 2 − r 2 , xy, xy), with spin σ (=↑, ↓), andn im ≡ ∑ σĉ † imσĉ imσ is the orbital density operator. The operatorsL i andŜ i represent the orbital and the spin angular momenta, respectively, and their expressions are given in Section 6.1. The parameters d m are the crystal field energies, while the hopping integrals t 0,ij mm are obtained as linear combinations of the Slater-Koster parameters [36,37] shown in Table 1. The spin-orbit interaction energy is set as ξ = 0.3129 eV. The number of electrons is one per site, and the Fermi energy or the chemical potential is determined accordingly.
The effect of the PLD is taken into account by changes in the transfer integrals [34]. Let r eq i denote the position of the i-th Ta atom in the absence of a PLD, while r i is the position when a PLD is present. Upon the change in atomic positions, we set the transfer integral between sites i and j as following the empirical law for the d-electrons. Thus, the Hamiltonian with taking account of a PLD is given byĤ Note that we neglect inhomogeneous changes in local potential ∆d im ({r j }) associated with the PLD.
Following the previous studies [34,35], we consider a specific PLD mode of the hexagram shape shown in Figure 1(a). Here, the unit cell has 6-fold rotation symmetry and contains 13 sites categorized into three types as shown in Figure 1(b); one "A" (red), six "B" (blue), and six "C" (orange) sites. The distances from A to B and C sites are denoted by r B and r C and these are parameterized by one parameter x quantifying the amplitude of the PLD as Here the length unit is the lattice constant 3.36 Å in the absence of a PLD. In [35], the value x = −1 corresponding to the values (r B , r C ) = (0.936, 1.66) is used to be consistent with experimental results. We note that x < 0 (x > 0) corresponds to a PLD with shrinking (expanding) hexagrams. In calculation, we generate all positions {r i } for a given x without shifting any A site, and obtain the transfer integrals t ij mm ({r i }) according to the modified distances |r i − r j | for all pairs of nearest neighbors. The band structure in our tight-binding model is summarized in Figure 2. Panel (a) shows the result calculated fromĤ 0 TB (1) in the absence of a PLD. Panel (b) shows the same result but the bands are folded by taking a hexagram as a unit cell. Panel (c) shows the result calculated fromĤ TB (3) with distortion x = −1.5. This case shows a large band restructuring, and one narrow band separates from the others and extends near the Fermi energy. This isolation becomes clearer for larger |x| with x < 0. Note that we consider the case of x < −1, while the value x = −1 was used in [35]. This is because the value |x| = 1 is not large enough to obtain a splitting of the narrow band within our treatment of the electron-electron interactions discussed below.

Electron-Electron Interactions
We now add toĤ TB (3) the electron-electron interactions of on-site Coulomb typê Here, we use the standard choice U = U − 2J and ignore the pair hopping terms [38]. For simplicity, throughout this paper, we focus on the typical case of U = 3J, where the last term in Equation (5) vanishes. We use a mean-field approach and approximateĤ where n imσ ≡ n imσ is determined self-consistently and we have neglected the orbital off-diagonal contributions ĉ † imσĉ im σ for simplicity. In the mean-field approximation, we use the unit cell made of two hexagrams, which accommodates 26 electrons, and it is important that the electron number is even. To simulate the charge distribution in the Mott insulator phase, we introduce a fictitious magnetic order in the mean-field approximation. Yet, this redundant magnetic order is not very harmful for discussing the charge dynamics driven by strong laser fields. Another drawback is that the translation symmetry is further lowered than a periodic array of hexagrams, but we shall show later that this effect is not very strong.
We compare in Figures 2(d) and (e) the band structure without and with the electron-electron interaction. Panel (d) is plotted for comparison and shows the same energy bands as in (c) (U = 3J = 0), but the bands are now folded corresponding to a two-hexagram unit cell. Figure 2(e) is the result for U = 3J = 0.816 eV calculated at zero temperature T = 0. The interactions have relatively weak effects at energies well below the Fermi energy E F . However, near E F , they cause a splitting of the two bands and make the system insulating. The band gap at the Γ point is 69.2 meV, which is a few times smaller than the experimental value [14]. While the gap size can be reproduced by tuning U, U , and J, we do not go into further detail in this paper.
Let us now examine charge and spin densities. The results are shown in Figure 3 for the same parameters U = 3J = 0.816 eV at T = 0. All the hexagrams have the identical charge density distribution: the density increases toward the center of each hexagram, and the 2-fold rotational symmetry around an A site is preserved. The 6-fold rotational symmetry is broken due to the doubled unit cell, but its asymmetry is so small (about 1%) that this effect is ignorable. In the doubled unit cell, the two hexagrams have the spin density distribution with the opposite signs, and this is due to a fictitious Neel order. The spin density distribution in each hexagram reflects the profile of the wave functions of the narrow band in Figure 2(c). In fact, as was shown in [35], these wave functions have a large amplitude near the hexagram center (A site). We emphasize again that the spin polarization appears due to the technical reason in treating a Mott insulator by the mean-field approximation (6). In fact, the CCDW phase of 1T-TaS 2 is paramagnetic [12,13]. Our aim is not to discuss the spin structure, but to investigate the charge properties, which play a main role upon strong laser drivings. For this purpose, we have now set up a reasonable model by adding the Coulomb interaction (5) to the tight-binding Hamiltonian (3) used in previous studies.

Lattice Degree of Freedom
To discuss cooperative phenomena, we now generalize our model and treat the lattice distortion x as a dynamical variable. We discuss its statics in this section and shall do its dynamics in Section 4.
We parametrize the elastic energy per two-hexagram unit cell accompanied by a PLD as where K 2 , K 3 , and K 4 are positive parameters so that V lat (x) is convex. Equation (7) is a Taylor series in x, which starts from x 2 since x = 0 represents the equilibrium position at high temperatures. We take account of the x 3 term, since the asymmetry between x < 0 and x > 0, or shrinking and expanding hexagrams, become nonnegligible for larger x. We set the parameters K 2 , K 3 , and K 4 so as to satisfy the following three conditions: (i) the CDW phase appears at T ∼ 500 K, (ii) x = −1.5 in equilibrium at T = 0, and (iii) the free energy difference between the CCDW (x = −1.5) and the uniform (x = 0) states is not very large, or typically less than 1 eV at T = 0 per two-hexagram unit cell. Here the free energy means the sum of V lat (x) and the electronic free energy F el (x) calculated forĤ TB (x) +Ĥ int within the mean-field approximation (6). The condition (i) determines K 2 , and then the conditions (ii) and (iii) do K 3 and K 4 . We have chosen K 2 = 12.1 eV, K 3 = 5.44 eV, and K 4 = 2.54 eV. Figure 4(a) shows the electronic energy F el (x) at T = 0 and V lat (x) as functions of x. We note that the free energy coincides with the energy at T = 0. The PLD decreases F el (x) and increases V lat (x) in both directions x < 0 and x > 0. The competition between these gain and cost results in the minimum of the total free energy F el (x) + V lat (x) at x = −1.5 as shown in Figure 4(b). There also exists a very shallow local minimum at x ∼ 0.2 corresponding to an expanding hexagram. At higher temperatures above ∼800 K, the total free energy F el (x) + V lat (x) has only one minimum at x = 0 and the CDW is not present. We note that this transition temperature is probably overestimated since we have performed the mean-field approximation (6) and neglected phonon fluctuations. We make a remark on the simplification that we have made on the PLD. We have treated one specific mode of the hexagram shape in lattice distortion, but the distortion has many more modes. In fact, 1T-TaS 2 shows an incommensurate CDW below 540 K and a nearly commensurate CDW below 350 K, which cannot be described within the present model. Thus we do not further discuss the free energy at intermediate temperatures, and shall focus on the CCDW state thus obtained at T = 0.

Equations of Motion
To study the photoinduced charge dynamics, we use the coupled equations of motion for the electrons and the PLD. We treat the electrons quantum-mechanically and describe their state by a mean-field wave function |Ψ(t) at each time t. For the PLD, we classically treat it and describe its coordinate and momentum, x(t) and p(t), by Hamilton's equations of motion.
The total Hamiltonian is given bŷ where M is the effective mass of the PLD and its value will be discussed later, andĤ int is treated within the mean-field approximation (6). The explicit time dependence ofĤ TB (x(t), t) comes from the coupling to the laser electric field, which oscillates in time. We assume that the electric field E(t) is spatially uniform and treat it in terms of the vector potential A(t) satisfying E(t) = −dA(t)/dt. Then H TB (x(t), t) is obtained from Equation (3) with the Peierls substitution where the elementary charge is set to unity.

of 13
The equations of motion are given by In the Schrödinger equation (12), we have ignored the c-number term H lat (x, p) on the right-hand side since it affects no physical observables. Equations (13) and (14) are Hamilton's equations of motion for the PLD. The right-hand side of Equation (14) represents the force acting on the PLD and consists of the elasticity and the reaction from the electrons. We note that, when E(t) = 0,Ĥ el (x(t), t) does not depend explicitly on t and the total energy Ψ(t)|Ĥ el (x(t))|Ψ(t) + H lat (x(t), p(t)) is a conserved quantity.

Cooperative Dynamics of Charge and Lattice
Motivated by recent experiments, we study the electronic excitation driven by an ultrashort laser pulse. We take the following form for the vector potential: Here the real parameter A 0 denotes the peak amplitude of the vector potential, and is the polarization vector and assumed to be = (1, 0) in this paper. The results are not sensitive to the direction of . Following the experiment [19], we set the central frequency ω pump and the pulse width T pump so that hω pump = 1.55 eV and 2T pump = 35 fs. The profile A x (t) is shown in Figure 5(a). We note that the peak amplitude of the electric field is approximately given by E 0 = ω pump A 0 since T pump 1/ω pump = 2.4 fs. In the following, we shall use E 0 rather than A 0 making it easier to compare our results with experimental results. average of (8,11) (9,12) (10,13)   (3,6), and (4, 7)], and the "C" sites [ (8,11), (9,12), and (10, 13)], respectively. Panels (d) and (e) show those of charge-density difference in the pairs of the "B" sites [ (2,5), (3,6), and (4, 7)], and the "C" sites [ (8,11), (9,12), and (10, 13)], respectively. Panel (f) shows the time evolution of the lattice distortion x(t).
The laser electric field induces charge dynamics on short time scale of order T pump . In Figure 5, we show the numerical results of local charge densities at 13 sites in a hexagram unit for a typical peak amplitude E 0 = 0.92 MV/cm, which is comparable with the experimental condition [19] (see Figure 1 for the labeling of sites). Since we are not interested in the spin density, the charge density averaged between the two hexagrams is plotted. At our initial time of calculation t = t init = −52.5 fs, the charge-density distribution is the one shown in Figure 3(a). As noted before, the 6-fold rotational symmetry is slightly broken down to the 2-fold symmetry. As time advances, the electric field envelope gradually increases, and the charge density starts to oscillate in time at each site.
To analyze the charge-density oscillations in detail, we pair up the six "B" sites as (2,5), (3,6), and (4, 7), each of which locates symmetrically about the central A site. We do the similar pairing for the "C" sites as well. We look into the average and the difference of the charge densities for each site pair, which are even and odd, respectively, under the 2-fold rotation. We note that the charge density at site 1 (A) itself is even under this operation. These even and odd quantities show different kinds of dynamics. The even quantities shown in panels (a), (b), and (c) oscillate with frequency ∼ 2ω pump , whereas the odd ones (d) and (e) do with frequency ∼ ω pump . This is because the vector potential A(t) has an odd parity with respect to the 2-fold rotation, and these two quantities are coupled to A(t) in the second and the first orders, respectively.
In addition to those oscillations, the charge dynamics changes its spatial distribution and the CDW order decreases on the short time scale of order T pump . The charge-density imbalance in panels (a-c) decreases during the laser irradiation, and the local charge density slightly approaches the uniform one.
This CDW order decrease is caused not just by the decrease in the PLD. Figure 5(f) shows that the lattice distortion x(t) remains almost unchanged until t ∼ 0, while the significant changes have already appeared in the charge density. Besides, at 0 < t < 50 fs, |x(t)| gradually decreases, while the charge dynamics rather slows down. Therefore, the short-time dynamics is dominated by the electrons driven by the laser field, and the CDW is rapidly suppressed.
On a longer time scale, the lattice distortion x(t) plays a main role in dynamics. We show x(t) up to t = 1000 fs in Fig 6. Its time evolution is well approximated by a harmonic oscillation starting at t ∼ 0 fs, and its frequency is determined as 2.14 THz by a sinusoidal fitting. This agrees well with the experimentally observed value 2.1 THz of the coherent phonon [15]. Of course, the frequency depends on the effective mass M of the PLD. We have set M = 8.00 × 10 3 u in the unified atomic mass unit to make the frequency consistent with the experimental result. The lattice dynamics can be interpreted by an instantaneous modulation of the free energy profile as schematically illustrated in Figure 6(b). In the initial state, the lattice distortion x is located at the minimum of the total free energy F eq . Then, on the short time scale, the laser drives the electrons to a high-energy state, and the free energy is changed into a new one F neq , whose minimum shifts toward x = 0. Since the lattice distortion remains almost unchanged during the driving, the free energy is not minimum and the distortion x oscillates harmonically around the new minimum of the free energy.
We remark that the harmonic oscillation is not damped in our model since we incorporate no dissipation processes and the total energy is conserved. In reality, the time scale of the damping is known to be ∼ 4 ps. Thus we do not proceed our analysis after t ∼ 1000 fs, where dissipation would start to play a main role in dynamics.
We briefly summarize the cooperative dynamics revealed by our analysis. The short-time dynamics is dominated by the electronic degrees of freedom driven by the laser. The driving rapidly suppresses the CDW order on the time scale of the pulse width, which is 35 fs in our analysis. In this time scale, the PLD is approximately frozen and does not decrease much. On the long time scale of order 100 fs, the PLD starts to oscillate at frequency ∼ 2 THz and will be damped in general through dissipation. These observations are consistent with the experimental finding [15], which observed a delay of the PLD decrease after the CDW order decrease.

Melting of Hexagrams
Now we further increase the laser amplitude and see how the lattice dynamics changes. Figure 7 shows the dynamics of the charge densities averaged over A, B, and C sites for the larger peak amplitudes E 0 = 4.61 and 23.1 MV/cm in panels (a) and (b), respectively. For E 0 = 4.61 MV/cm, the CDW order decreases by ∼ 20%, which is much enhanced than the case of E 0 = 0.92 MV/cm shown in Figure 5. For E 0 = 23.1 MV/cm, the CDW order disappears and the charge distribution becomes nearly uniform within ∼ 20 fs. We note that the PLD remains present on this time scale as shown in panel (c). Only the CDW order of the hexagram shape melts down on the short time scale. Next we discuss the energy efficiency of the photoinduced CDW melting. We have calculated the energy absorption ∆E abs ≡ E(t = 0 fs) − E init , and obtained ∆E abs = 4.6 and 39 eV per two hexagrams for E 0 = 4.61 and 23.1 MV/cm, respectively. We compare these values with the energy difference ∆E CDW ≡ E uniform − E init . Here E uniform denotes the energy of a virtual uniform state and is defined as the expectation value ofĤ total (8) at x = −1.5 calculated for the electronic ground state ofĤ total for x = 0. We obtain ∆E CDW = 5.8 eV per two hexagrams, and it is an estimate of the energy needed to wash out the CDW order without deforming the PLD. It is consistent with the fact that ∆E abs < ∆E CDW for E 0 = 4.61 MV/cm, where the CDW is suppressed but remains, whereas ∆E abs > ∆E CDW = 0.15∆E abs for E 0 = 23.1 MV/cm, where the CDW melts down. Thus 15% of the absorbed energy is used for the CDW melting, and the remaining 85% is converted into heat.
To melt the CDW order, the laser-electron interaction is not the only way since the lattice motion is involved on the long time scale. We show the lattice dynamics in Figure 7 for some strong laser fields. For E 0 = 4.61 MV/cm, x(t) hits zero, meaning that the PLD of the hexagram shape melts down and the CDW order of electrons thereby disappears. We note that this field amplitude is not strong enough to completely melt the CDW order on the short time scale. However, in the long run, the laser energy absorbed by the electrons is transferred to the lattice motion, which affects the charge density through the electron-lattice coupling. On the long time scale, the electron-lattice cooperative dynamics becomes more important.
These electron-lattice dynamics are quite different from those found in experimental [39] and theoretical studies [40,41] on charge order melting in quasi-two-dimensional organic conductors. In these materials, the charge order is mainly produced by long-range Coulomb repulsion and weakly assisted by an electron-lattice coupling. Thus, the charge order can be melted efficiently by a pulse laser with little changing the lattice distortion. See also [42] for the efficiency of the charge order melting in a quasi-one-dimensional organic conductor.

Discussion and Conclusion
In this work, we have extended the tight-binding model (3) of the previous studies [34,35] and included the on-site Coulomb interactions (5). By this extension and our mean-field analysis, we have succeeded in opening a band gap at the Fermi level as shown in Figure 2(e). Thus we have established a lattice model which reproduces qualitative features of the electronic properties of 1T-TaS 2 . We have further extended the model to treat the lattice distortion as a dynamical variable. Considering a specific PLD of the hexagram shape, we have modeled the lattice potential, reproduced the CCDW as a thermal equilibrium state, and formulated the coupled equations of motion for the electrons and the PLD.
We have investigated the short-time charge dynamics induced by a pulse laser on the basis of our model. We have shown that the short-time dynamics is dominated by electrons and the hexagram CDW order is suppressed by the laser, while the lattice is effectively frozen. This CDW decrease is caused at the second order of the laser electric field, and the CDW order melts down as fast as 20 fs for the strong laser amplitude E 0 = 23.1 MV/cm for the 35 fs pulse. We note that the necessary amplitude E 0 should be smaller for a longer pulse.
The lattice distortion plays an important role on the long time scale 100 fs. Its motion is dominated by the landscape of the free energy instantaneously converted from the equilibrium one [see Figure 6(b)]. The CDW order may be washed out on the long time scale through the coupling to the lattice motion, even when the laser amplitude is not strong enough to melt the CDW on the short time scale.
We comment on a few open questions beyond our present study. We have considered only one mode of the hexagram shape in lattice distortion, but there exist many other modes giving rise to the incommensurate CDW, the nearly-commensurate CDW, and so on. Including those relevant modes may make it possible to describe photoinduced transitions from the CCDW to the other CDW phases. In that case, our picture of the free-energy modulation is extended in multiple dimensions and several local minima may exist. We note that the phonon modes also work as a heat reservoir. Our present analysis does not include a coupling to any reservoir, or dissipation, and our results would be modified especially at long times. We leave these important issues as future works.

Spin-Orbit Interaction
We give the the expressions forL i andŜ i : Here τ α 's are the Pauli matrices and L α mm 's are the 3 × 3 sub-matrices for m = 3z 2 − r 2 , x 2 − y 2 , xy of the L = 2 representation of the angular momentum [43]. We note that onlyL z iŜ z i is active inL i ·Ŝ i within our subspace spanned by the states with L z i = 0, ±2. The nonzero elements of L z are thus L z x 2 −y 2 ,xy = −L z xy,x 2 −y 2 = i.