Ultrafast charge ordering by self-amplified exciton-phonon dynamics in   TiSe$_2$

Chao Lian; Sheng-Jie Zhang; Shi-Qi Hu; Meng-Xue Guan; and Sheng Meng

arXiv:1901.00610·cond-mat.mes-hall·January 7, 2020

Ultrafast charge ordering by self-amplified exciton-phonon dynamics in TiSe$_2$

Chao Lian, Sheng-Jie Zhang, Shi-Qi Hu, Meng-Xue Guan, and Sheng Meng

PDF

TL;DR

This study uses real-time simulations to disentangle and clarify the ultrafast interplay between excitonic and lattice dynamics in TiSe$_2$, revealing how they mutually amplify to invert charge density wave order.

Contribution

It provides the first detailed temporal separation of exciton and phonon contributions in TiSe$_2$, demonstrating their sequential and self-amplifying roles in CDW dynamics.

Findings

01

Exciton order is disrupted within 20 fs by laser pulses.

02

Electronic order weakening triggers ionic movements opposite to original lattice distortions.

03

The interplay leads to a complete inversion of charge density wave order.

Abstract

The origin of charge density waves (CDW) in TiSe $_{2}$ has long been debated, mainly due to the difficulties in identifying the timescales of how and when the excitonic pairing and electron-phonon coupling (EPC) come into play. Without a proper time resolution and microscopic mechanism, one has to assume simultaneous appearance of CDW and periodic lattice distortions (PLD). Here, we accomplish a complete separation of exciton and PLD dynamics and unravel their interplay in the ultrafast time domain in our real-time time-dependent density functional theory simulations. We find that laser pulses knock off the exciton order and induce a homogeneous bonding-antibonding transition in the initial 20 fs, then the weakened electronic order triggers ionic movements antiparallel to the original PLD. The EPC comes into play after the initial 20~fs, and the two processes mutually amplify each other…

Equations60

i ℏ \frac{\partial ψ _{γ k} ( G , t )}{\partial t} = H_{k} (t) ψ_{γ k} (G, t)

i ℏ \frac{\partial ψ _{γ k} ( G , t )}{\partial t} = H_{k} (t) ψ_{γ k} (G, t)

H_{k} (G, G^{'}, t) = = T_{k} (G, G^{'}, t) + V (G, G^{'}, t) \frac{ℏ ^{2}}{2 m} ∣ k + G + A (t) ∣^{2} δ_{G, G^{'}} + V (G, G^{'}, t)

H_{k} (G, G^{'}, t) = = T_{k} (G, G^{'}, t) + V (G, G^{'}, t) \frac{ℏ ^{2}}{2 m} ∣ k + G + A (t) ∣^{2} δ_{G, G^{'}} + V (G, G^{'}, t)

A (t) = - \int_{0}^{t} E (t^{'}) d t^{'}

A (t) = - \int_{0}^{t} E (t^{'}) d t^{'}

∣ ψ_{γ k} (G, t) ⟩ = i \sum c_{iγ, k} (t) ∣ ϕ_{i k} (G, t_{1}) ⟩

∣ ψ_{γ k} (G, t) ⟩ = i \sum c_{iγ, k} (t) ∣ ϕ_{i k} (G, t_{1}) ⟩

H_{k} (G, t_{1}) ∣ ϕ_{i k} (G, t_{1}) ⟩ = ϵ_{i k} (t_{1}) ∣ ϕ_{i k} (G, t_{1}) ⟩,

H_{k} (G, t_{1}) ∣ ϕ_{i k} (G, t_{1}) ⟩ = ϵ_{i k} (t_{1}) ∣ ϕ_{i k} (G, t_{1}) ⟩,

∣ ψ_{i k} (G, t = 0) ⟩ := c_{iγ, k} (t = 0) := ∣ ϕ_{i k} (G, t_{1} = 0) ⟩, δ_{iγ},

∣ ψ_{i k} (G, t = 0) ⟩ := c_{iγ, k} (t = 0) := ∣ ϕ_{i k} (G, t_{1} = 0) ⟩, δ_{iγ},

H_{k} (t) C_{k} (t) = i ℏ \frac{\partial}{\partial t} C_{k} (t),

H_{k} (t) C_{k} (t) = i ℏ \frac{\partial}{\partial t} C_{k} (t),

{C}_{\mathbf{k}}(t)=\left(\begin{array}[]{cccc}c_{11,\mathbf{k}}(t)&c_{12,\mathbf{k}}(t)&\cdots&c_{1N_{b},\mathbf{k}}(t)\\ c_{21,\mathbf{k}}(t)&c_{12,\mathbf{k}}(t)&\cdots&c_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ c_{N_{b}1,\mathbf{k}}(t)&c_{N_{b}2,\mathbf{k}}(t)&\cdots&c_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

{C}_{\mathbf{k}}(t)=\left(\begin{array}[]{cccc}c_{11,\mathbf{k}}(t)&c_{12,\mathbf{k}}(t)&\cdots&c_{1N_{b},\mathbf{k}}(t)\\ c_{21,\mathbf{k}}(t)&c_{12,\mathbf{k}}(t)&\cdots&c_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ c_{N_{b}1,\mathbf{k}}(t)&c_{N_{b}2,\mathbf{k}}(t)&\cdots&c_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

{H}_{\mathbf{k}}(t)=\left(\begin{array}[]{cccc}h_{11,\mathbf{k}}(t)&h_{12,\mathbf{k}}(t)&\cdots&h_{1N_{b},\mathbf{k}}(t)\\ h_{21,\mathbf{k}}(t)&h_{12,\mathbf{k}}(t)&\cdots&h_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ h_{N_{b}1,\mathbf{k}}(t)&h_{N_{b}2,\mathbf{k}}(t)&\cdots&h_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

{H}_{\mathbf{k}}(t)=\left(\begin{array}[]{cccc}h_{11,\mathbf{k}}(t)&h_{12,\mathbf{k}}(t)&\cdots&h_{1N_{b},\mathbf{k}}(t)\\ h_{21,\mathbf{k}}(t)&h_{12,\mathbf{k}}(t)&\cdots&h_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ h_{N_{b}1,\mathbf{k}}(t)&h_{N_{b}2,\mathbf{k}}(t)&\cdots&h_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

H_{k} (t) = H_{k} (t_{1}) + \frac{t - t _{1}}{t _{2} - t _{1}} [H_{k} (t_{2}) - H_{k} (t_{1})],

H_{k} (t) = H_{k} (t_{1}) + \frac{t - t _{1}}{t _{2} - t _{1}} [H_{k} (t_{2}) - H_{k} (t_{1})],

h_{ij, k} (t_{1}) = ⟨ ϕ_{i} (t_{1}) ⟩ H_{k} (t_{1}) ϕ_{j} (t_{1}) = δ_{ij} ϵ_{i, k} (t_{1}),

h_{ij, k} (t_{1}) = ⟨ ϕ_{i} (t_{1}) ⟩ H_{k} (t_{1}) ϕ_{j} (t_{1}) = δ_{ij} ϵ_{i, k} (t_{1}),

H_{ij, k} (t_{2}) = = = = = ⟨ ϕ_{i, k} (t_{1}) ⟩ H_{k} (t_{2}) ϕ_{j, k} (t_{1}) ⟨ ϕ_{i, k} (t_{1}) ⟩ l \sum ϕ_{l, k} (t_{2}) ⟨ ϕ_{l, k} (t_{2}) ⟩ H_{k} (t_{2}) m \sum ϕ_{m, k} (t_{2}) ⟨ ϕ_{m, k} (t_{2}) ⟩ ϕ_{j, k} (t_{1}) l \sum m \sum a_{i l, k} (t_{1}, t_{2}) a_{mj, k}^{*} (t_{1}, t_{2}) ⟨ ϕ_{l, k} (t_{2}) ⟩ H_{k} (t_{2}) ϕ_{m, k} (t_{2}) l \sum m \sum a_{i l, k} (t_{1}, t_{2}) δ_{l m} ϵ_{l k} (t_{2}) a_{i l, k} (t_{1}, t_{2}) a_{mj, k}^{*} (t_{1}, t_{2}) l \sum a_{i l, k} (t_{1}, t_{2}) a_{l j, k}^{*} (t_{1}, t_{2}) ϵ_{l k} (t_{2}),

H_{ij, k} (t_{2}) = = = = = ⟨ ϕ_{i, k} (t_{1}) ⟩ H_{k} (t_{2}) ϕ_{j, k} (t_{1}) ⟨ ϕ_{i, k} (t_{1}) ⟩ l \sum ϕ_{l, k} (t_{2}) ⟨ ϕ_{l, k} (t_{2}) ⟩ H_{k} (t_{2}) m \sum ϕ_{m, k} (t_{2}) ⟨ ϕ_{m, k} (t_{2}) ⟩ ϕ_{j, k} (t_{1}) l \sum m \sum a_{i l, k} (t_{1}, t_{2}) a_{mj, k}^{*} (t_{1}, t_{2}) ⟨ ϕ_{l, k} (t_{2}) ⟩ H_{k} (t_{2}) ϕ_{m, k} (t_{2}) l \sum m \sum a_{i l, k} (t_{1}, t_{2}) δ_{l m} ϵ_{l k} (t_{2}) a_{i l, k} (t_{1}, t_{2}) a_{mj, k}^{*} (t_{1}, t_{2}) l \sum a_{i l, k} (t_{1}, t_{2}) a_{l j, k}^{*} (t_{1}, t_{2}) ϵ_{l k} (t_{2}),

a_{i l, k} (t_{1}, t_{2}) = ⟨ ϕ_{i, k} (t_{1}) ⟩ ϕ_{l, k} (t_{2})

a_{i l, k} (t_{1}, t_{2}) = ⟨ ϕ_{i, k} (t_{1}) ⟩ ϕ_{l, k} (t_{2})

U_{k} (t_{2}, t_{1}) = s = 0 \prod N_{t} \frac{exp [ - i ℏ H ( t _{1} + s d t ) d t /2 ]}{exp [ i ℏ H ( t _{1} + s d t ) d t /2 ]} = s = 0 \prod N_{t} \frac{1 - i ℏ H _{k} ( t _{1} + s d t ) d t /2}{1 + i ℏ H _{k} ( t _{1} + s d t ) d t /2},

U_{k} (t_{2}, t_{1}) = s = 0 \prod N_{t} \frac{exp [ - i ℏ H ( t _{1} + s d t ) d t /2 ]}{exp [ i ℏ H ( t _{1} + s d t ) d t /2 ]} = s = 0 \prod N_{t} \frac{1 - i ℏ H _{k} ( t _{1} + s d t ) d t /2}{1 + i ℏ H _{k} ( t _{1} + s d t ) d t /2},

C_{k} (t_{2}) = U_{k} (t_{2}, t_{1}) C_{k} (t_{1}) .

C_{k} (t_{2}) = U_{k} (t_{2}, t_{1}) C_{k} (t_{1}) .

∣ ψ_{γ k} (t_{2}) ⟩ = i \sum c_{iγ, k} (t_{2}) ∣ ϕ_{i, k} (t_{1}) ⟩ = i \sum c_{iγ, k} (t_{2}) l \sum ∣ ϕ_{l, k} (t_{2}) ⟩ ⟨ ϕ_{l, k} (t_{2}) ⟩ ϕ_{i, k} (t_{1}) = l \sum i \sum a_{i l, k}^{*} (t_{1}, t_{2}) c_{iγ, k} (t_{2}) ∣ ϕ_{l, k} (t_{2}) ⟩ = l \sum c_{l γ, k}^{'} (t_{2}) ∣ ϕ_{l, k} (t_{2}) ⟩

∣ ψ_{γ k} (t_{2}) ⟩ = i \sum c_{iγ, k} (t_{2}) ∣ ϕ_{i, k} (t_{1}) ⟩ = i \sum c_{iγ, k} (t_{2}) l \sum ∣ ϕ_{l, k} (t_{2}) ⟩ ⟨ ϕ_{l, k} (t_{2}) ⟩ ϕ_{i, k} (t_{1}) = l \sum i \sum a_{i l, k}^{*} (t_{1}, t_{2}) c_{iγ, k} (t_{2}) ∣ ϕ_{l, k} (t_{2}) ⟩ = l \sum c_{l γ, k}^{'} (t_{2}) ∣ ϕ_{l, k} (t_{2}) ⟩

c_{l γ, k}^{'} (t_{2}) = i \sum a_{i l, k}^{*} (t_{1}, t_{2}) c_{iγ, k} (t_{2}) .

c_{l γ, k}^{'} (t_{2}) = i \sum a_{i l, k}^{*} (t_{1}, t_{2}) c_{iγ, k} (t_{2}) .

ρ (G, t_{2}) = k \sum γ \sum ∣ ψ_{γ, k} (G, t_{2}) ∣^{2} = k \sum γ \sum i \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2} ∣ ϕ_{i, k} (G, t_{2}) ∣^{2} = k \sum i \sum [γ \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2}] ∣ ϕ_{i k} (G, t_{2}) ∣^{2} = k \sum i \sum q_{i k} (t_{2}) ∣ ϕ_{i k} (G) ∣^{2},

ρ (G, t_{2}) = k \sum γ \sum ∣ ψ_{γ, k} (G, t_{2}) ∣^{2} = k \sum γ \sum i \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2} ∣ ϕ_{i, k} (G, t_{2}) ∣^{2} = k \sum i \sum [γ \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2}] ∣ ϕ_{i k} (G, t_{2}) ∣^{2} = k \sum i \sum q_{i k} (t_{2}) ∣ ϕ_{i k} (G) ∣^{2},

q_{i k} (t_{2}) = γ \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2}

q_{i k} (t_{2}) = γ \sum ∣ c_{iγ, k}^{'} (t_{2}) ∣^{2}

F_{R_{I}} = i k \sum ⟨ ψ_{i k} ⟩ \nabla_{R_{I}} H ψ_{i k},

F_{R_{I}} = i k \sum ⟨ ψ_{i k} ⟩ \nabla_{R_{I}} H ψ_{i k},

M_{I} \frac{d ^{2} R _{I}}{d t ^{2}} = F_{R_{I}},

M_{I} \frac{d ^{2} R _{I}}{d t ^{2}} = F_{R_{I}},

v_{I}^{'} (t) = v_{I} (t) T^{'} (t) / T (t),

v_{I}^{'} (t) = v_{I} (t) T^{'} (t) / T (t),

∣ ϕ_{i k} (G, t_{1}) ⟩ = \hat{T} ∣ \tilde{ϕ}_{i k} (G, t_{1}) ⟩

∣ ϕ_{i k} (G, t_{1}) ⟩ = \hat{T} ∣ \tilde{ϕ}_{i k} (G, t_{1}) ⟩

a_{i l, k} (t_{1}, t_{2}) = = ⟨ \tilde{ϕ}_{i k} (t_{1}) ⟩ \hat{T}^{†} \hat{T} \tilde{ϕ}_{l k} (t_{2}) ⟨ \tilde{ϕ}_{i k} (t_{1}) ⟩ \hat{S} \tilde{ϕ}_{l k} (t_{2})

a_{i l, k} (t_{1}, t_{2}) = = ⟨ \tilde{ϕ}_{i k} (t_{1}) ⟩ \hat{T}^{†} \hat{T} \tilde{ϕ}_{l k} (t_{2}) ⟨ \tilde{ϕ}_{i k} (t_{1}) ⟩ \hat{S} \tilde{ϕ}_{l k} (t_{2})

\hat{S}=\left(\begin{array}[]{cccc}s_{11,\mathbf{k}}(t)&s_{12,\mathbf{k}}(t)&\cdots&s_{1N_{b},\mathbf{k}}(t)\\ s_{21,\mathbf{k}}(t)&s_{12,\mathbf{k}}(t)&\cdots&s_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ s_{N_{b}1,\mathbf{k}}(t)&s_{N_{b}2,\mathbf{k}}(t)&\cdots&s_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

\hat{S}=\left(\begin{array}[]{cccc}s_{11,\mathbf{k}}(t)&s_{12,\mathbf{k}}(t)&\cdots&s_{1N_{b},\mathbf{k}}(t)\\ s_{21,\mathbf{k}}(t)&s_{12,\mathbf{k}}(t)&\cdots&s_{1N_{b},\mathbf{k}}(t)\\ \vdots&\vdots&\cdots&\vdots\\ s_{N_{b}1,\mathbf{k}}(t)&s_{N_{b}2,\mathbf{k}}(t)&\cdots&s_{N_{b}N_{b},\mathbf{k}}(t)\\ \end{array}\right),

∣ ϕ_{i, k} (G, t) ⟩ = I, K \sum a (I, K; i, k; t) ∣ Φ_{I, K} G, t) ⟩,

∣ ϕ_{i, k} (G, t) ⟩ = I, K \sum a (I, K; i, k; t) ∣ Φ_{I, K} G, t) ⟩,

A (K, E, t) = i \sum P (K; k, i; t) δ (E - ϵ_{i, k} (t)),

A (K, E, t) = i \sum P (K; k, i; t) δ (E - ϵ_{i, k} (t)),

P (K; k, i) = I \sum a^{*} (K, I; k, i; t) a (K, I; k, i; t) = G \sum ∣ ϕ_{i, k} (G + K - k, t) ∣^{2},

P (K; k, i) = I \sum a^{*} (K, I; k, i; t) a (K, I; k, i; t) = G \sum ∣ ϕ_{i, k} (G + K - k, t) ∣^{2},

A (K, E, t) = i \sum P (K; k, i; t) w_{i, k} (t) δ (E - ϵ_{i, k} (t)) .

A (K, E, t) = i \sum P (K; k, i; t) w_{i, k} (t) δ (E - ϵ_{i, k} (t)) .

w_{i, k} (t) = q_{i, k} (t)

w_{i, k} (t) = q_{i, k} (t)

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Ultrafast charge ordering by self-amplified exciton-phonon dynamics in TiSe2

Chao Lian

Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing, 100190, P. R. China

Sheng-Jie Zhang

Shi-Qi Hu

Meng-Xue Guan

Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing, 100190, P. R. China

Sheng Meng

Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing, 100190, P. R. China

School of Physical Sciences, University of Chinese Academy of Sciences, Beijing, 100190, P. R. China

Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, P. R. China

Abstract

The origin of charge density waves (CDW) in TiSe2 has long been debated, mainly due to the difficulties in identifying the timescales of how and when the excitonic pairing and electron-phonon coupling (EPC) come into play. Without a proper time resolution and microscopic mechanism, one has to assume simultaneous appearance of CDW and periodic lattice distortions (PLD). Here, we accomplish a complete separation of exciton and PLD dynamics and unravel their interplay in the ultrafast time domain in our real-time time-dependent density functional theory simulations. We find that laser pulses knock off the exciton order and induce a homogeneous bonding-antibonding transition in the initial 20 fs, then the weakened electronic order triggers ionic movements antiparallel to the original PLD. The EPC comes into play after the initial 20 fs, and the two processes mutually amplify each other leading to a complete inversion of CDW ordering. The self-amplified dynamics reproduces the evolution of band structures in excellent agreement with ultrafast photoemission experiment. Hence we resolve the key processes in the initial dynamics of CDW that help elucidate the mechanism underlying the long debated problem.

Charge density wave (CDW) in 1T-TiSe2 has been one of the persistent eye-drawing topics over decades. It is not only an excellent playground to study the interplay between CDW and superconductivity Morosan2006 ; Li2007 ; Kusmartseva2009 ; Yao2018 ; Wei2017b ; Kogar2017a ; Yan2017 ; Medvecka2016 ; Das2015 ; Luna2015 ; Ganesh2014 ; Joe2014 ; Kacmarcik2013 ; Husanikova2013 ; Iavarone2012 ; Jeong2010 ; Zaberchik2010 ; Giang2010 ; Morosan2010 ; Hillier2010 ; Kusmartseva2009 ; Barath2008 ; Li2007b ; Li2007 ; Morosan2006 ; Cui2006 , but also an evidenced excitonic insulator Kohn1967 ; Jerome1967 ; Halperin1968 ; Rossnagel2002 ; Cercellier2007 . It was heavily debated which mechanism – the electron-phonon coupling (EPC) Hughes1977 ; Wakabayashi1978 ; Gaby1981 ; Motizuki1981 ; Lopez-Castillo1987 ; Holt2001 ; Bussmann-Holder2009 ; VanWezel2010a ; Rossnagel2010 ; Calandra2011 ; Zhu2012 ; Zenker2013 or the excitonic pairing Kidd2002 ; Monney2012a ; VanWezel2010 ; Stoffel1982 ; Monney2015 ; Cazzaniga2012 ; Monney2010 ; Anderson1985a ; Zenker2013 ; Peng2015b ; Koley2014 ; Novello2017 ; Monney2012 ; Monney2010a ; Monney2011 ; Monney2009 ; Sugawara2016 ; Monney2012a ; Hildebrand2016a ; Watanabe2015 ; Pillo2000 ; May2011 ; Chen2018c – is the major driving force for the formation of CDW in TiSe2. These two mechanisms disagree on the role of periodic lattice distortions (PLD) in CDW: PLD is essential in forming the CDW according to the EPC mechanism, while it is only a passive consequence of the CDW if the excitonic pairing dominates.

Isolating the PLD from the CDW can solve the debate, but it is not achievable in the ground state. Inspiringly, the ultrafast measurements can distinguish PLD and CDW based on their different time scales Sundaram2002 . However, conclusions from previous ultrafast measurements are highly controversial. Utilizing time-resolved (tr) and angle-resolved photoemission spectroscopy (ARPES) measurements, Rohwer et al. observed very fast ( $<30$ fs) collapses of the long-range order Rohwer2011a , which was interpreted by Mathias et al. as the signals of the mutually-amplified carrier multiplication and gap quenching Mathias2016 . However, this mechanism completely neglects the possible involvements of ionic dynamics and the effects of geometry relaxation. Möhr-Vorobeva et al. observed the nonthermal melting within 250 fs in the ultrafast X-ray measurements and suggested that excitonic pairing generates the CDW Mohr-Vorobeva2011 . Hellmann et al. also supported the excitonic mechanism Hellmann2012 . They observed that the $100$ fs change in the tr-ARPES signals is comparable to the buildup time of the electron-hole screening in the exciton formation. On the other hand, Porer et al. indicated that excitonic pairing was not the sole driving force of CDW Porer2014 . They separately studied the electronic and structural orders via monitoring the characteristic peaks of the two in the transient energy loss spectra. They found that the PLD can persist with the quenched excitonic order. Despite extensive efforts in the past, the conclusions in all these previous studies were mainly derived from indirect mappings between the spectra and the CDW/PLD orders. Notwithstanding the cutting-edge techniques used in these studies, the time-resolved spectra can only provide the averaged response of the material without atomic resolution. With recent significant progresses in real-time (rt) time dependent (TD) density functional theory (DFT) algorithms and computing power Runge1984 ; Bertsch2000 ; Yabana2006 ; Otobe2008 ; Otobe2009 ; Otobe2016b ; Yabana2012 ; Shinohara2010 ; Shinohara2010a ; Shinohara2012 ; Sato2015 ; Sato2015a ; Ren2013 ; Wang2015a , it become now possible to perform ultrafast quantum dynamics simulations fully from first principles, to provide a unified atomic-level picture of ultrafast CDW dynamics.

Here we take advantages of recently developed rt-TDDFT methods for ab initio simulations of ultrafast dynamics in complex materials such as TiSe2. We demonstrate that laser pulses knock down the CDW order without disturbing PLD through inducing a homogeneous bonding-antibonding electronic transition. The reduced CDW order then triggers the ionic movements exactly antiparallel to the original PLD, but cannot solely drives the observed inversion in CDW/PLD. Instead, assisted by electron-phonon couplings, a self-amplification mechanism between electron dynamics and lattice distortion emerges after the initial excitation, reproducing well experimental features observed in tr-ARPES measurements. We propose that both exciton pairing and electron-phonon couplings contribute to the CDW formation, albeit in a different timescale (see Fig. 1): CDW is predominantly initiated by exciton binding ( $<20$ fs) and subsequently enhanced by EPC ( $>20$ fs). These insights hint for a complete microscopic understanding on the nature of charge ordering in quantum materials.

**Results

Ground state properties.** We calculate the atomic geometries of the normal phase and the CDW 2 $\times$ 2 $\times$ 2 phase. The bulk TiSe2 is used with the interlayer separation of 6.69 Å. The optimized PLD displacements { $\mathbf{d}_{i}$ } are shown in Fig. 2(a), where the displacements $\delta_{Ti}=0.091$ Å, $\delta_{Se}=0.030$ Å and $\delta_{Ti}/\delta_{Se}=3.03:1$ . These results excellently reproduce the experimental measurements $\delta_{Ti}=0.085\pm 0.014$ Å and $\delta_{Ti}/\delta_{Se}\sim 3:1$ DiSalvo1976 . To directly compare with ARPES measurements, we calculate the effective band structure (EBS) along K-M- $\mathrm{\Gamma}$ -K by unfolding energy bands from the $2\times 2$ BZ to the $1\times 1$ BZ, as shown in Fig. 2. The calculated EBS agrees well with the experimental spectra. In the normal state, we reproduce the electron and hole pockets at M and $\Gamma$ points, respectively. The CDW opens an indirect gap of 0.18 eV and creates folded valence bands (VB) at M and folded conduction bands (CB) at $\Gamma$ . We note that PBE+U with $U_{\mathrm{Ti}}=3.5$ eV incorrectly indicates that the CDW state is unstable Bianco2015 . Thus, we use $U_{\mathrm{Ti}}=0$ eV for all the dynamics simulations with TDDFT.

Based on the consistency between DFT results and experimental measurements, we briefly discuss the accuracy of DFT in describing the excitonic paring and EPC in TiSe2. It is well known that the semilocal exchange-correlation (XC) functionals (e.g. PBE) poorly describes the long-range Coulomb screening $2\pi/|\mathbf{k}^{\prime}-\mathbf{k}|\sim\infty$ in vertical excitations $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}=0$ Sharma2011 . Computationally expensive corrections such as Bethe-Salpeter equation (BSE) Sottile2003 ; Marini2003 can considerably improve the accuracy. The experimentally-observed superlinear feature is absent in our simulations [Supplementary Fig. 2]. However, the long-range attractions produce spatially uniform forces on the ions. Based on the concept of the excitonic insulator Kohn1967 ; Jerome1967 ; Halperin1968 , the PLD stability is only affected by the inter-valley excitons formed by an attractive interaction $V(\mathbf{w})$ between the electron pocket at the M point and the hole pocket at the $\Gamma$ point, as shown in Fig. 2(d). Here, $\mathbf{w}=\pm\mathbf{b}_{i}/2$ and $\mathbf{b}_{i}$ ( $i=1,2$ ) being the reciprocal lattice vector along the $i$ th direction. Therefore, $V(\mathbf{w})$ is a short-range interaction with a characteristic length scale $1/|\mathbf{w}|=a$ , where $a=|\mathbf{a}_{i}|$ and $\mathbf{a}_{i}$ is the lattice vector. Note $V(\mathbf{w})$ is different from typical long-range interactions in Wannier excitons. The semilocal XC already includes exciton binding between the electron and hole pockets at M and $\Gamma$ , respectively, albeit slightly underestimating the screening effect.

To quantitatively demonstrate the validity of the semilocal functional in describing the intervalley exciton with momentum $\mathbf{q}=\mathbf{w}$ , we compare the linear-response TDDFT results obtained from adiabatic PBE (APBE) and BSE kernels. At $\mathbf{q}=\mathbf{w}$ , the APBE kernel yields similar absorption spectra with those from BSE kernel [Supplementary Fig. 3]. Since the BSE kernel is a well-accepted accurate description of excitons Sottile2003 ; Marini2003 , this indicates that the electron-hole exchange effect has been well described in the semilocal XC. Therefore, due to the unique band structures of TiSe2, the semilocal XC yields acceptable excitonic interactions.

Furthermore, the density functional perturbation theory and molecular dynamics calculations reproduce well experimental phonon spectra Weber2011 ; Singh2017a and thermal conductivity Bianco2015 ; Hellgren2017 , confirming that DFT can accurately describe the EPC in TiSe2. Thus, the semilocal XC and electron-ion dynamics simulations are suitable for tracking photoexcitation physics of TiSe2.

CDW dynamics. To identify the dominant mechanism, we try to supply direct proofs by separating the CDW dynamics from the PLD in the time domain. Instead of directly creating CDW at the normal state, we aim at decreasing the CDW order with the fixed PLD. We analyze the differential charge density $\rho_{\mathrm{tot}}(\mathbf{r},t)=\rho_{\mathrm{chg}}(\mathbf{r},t)-\rho_{\mathrm{atom}}(\mathbf{r},t)$ as a function of time, where $\rho_{\mathrm{chg}}$ is the charge density and $\rho_{\mathrm{atom}}$ is the superposition of the atomic charge densities. Thus, $\rho_{\mathrm{tot}}(\mathbf{r},t)$ features the spacial distribution of the bonding (+) and antibonding (–) densities. We divide $\rho_{\mathrm{tot}}(\mathbf{r},t)$ into two parts $\rho_{1\times 1}(\mathbf{r})$ and $\rho_{\mathrm{CDW}}(\mathbf{r})$ , where $\rho_{1\times 1}(\mathbf{r})=[\rho_{\mathrm{tot}}(\mathbf{r}+\mathbf{a}_{i})+\rho_{\mathrm{tot}}(\mathbf{r})]/2$ , and $\rho_{\mathrm{CDW}}(\mathbf{r})=\rho_{\mathrm{tot}}(\mathbf{r})-\rho_{1\times 1}(\mathbf{r})$ . Obviously, $\rho_{1\times 1}(\mathbf{r}+\mathbf{a}_{i})=\rho_{1\times 1}(\mathbf{r})$ characterizes the original $1\times 1$ order, while $\rho_{\mathrm{CDW}}(\mathbf{r})$ characterizes the strength of the CDW order. A laser induced charge difference is characterized by $\Delta\rho_{i}(\mathbf{r})=\rho_{i}(\mathbf{r},t_{f})-\rho_{i}(\mathbf{r},t_{0})$ ( $i=\mathrm{tot}$ , $1\times 1$ , $\mathrm{CDW}$ ).

Figure 3 shows the two-dimensional contour of the charge density $\rho_{i}(x,y,t)$ ( $i=\mathrm{tot}$ , $1\times 1$ , $\mathrm{CDW}$ ), which is $\rho_{i}(\mathbf{r},t)$ averaged over $z$ direction. Comparing Fig. 3(a)-(c), we find that $\rho_{1\times 1}(x,y,t_{0})$ is the main ingredient of $\rho_{\mathrm{tot}}(x,y,t_{0})$ , even with the presence of PLD, while $\rho_{\mathrm{CDW}}(x,y,t_{0})$ is localized around the Ti positions. Starting from the CDW ground state with the fixed PLD, the CDW order $\rho_{\mathrm{CDW}}(x,y,t)$ decreases after the laser illumination: (i) The laser induces electron transfer from the bonding area to the antibonding area. The induced charge $\Delta\rho_{\mathrm{tot}}(x,y,t)$ [Fig. 3(g)] is opposite to ground state charge $\rho_{\mathrm{tot}}(x,y,t_{0})$ [Fig. 3(a)]. (ii) The majority of the induced charge has the $1\times 1$ periodicity, i.e. $\Delta\rho_{\mathrm{1\times 1}}(x,y,t)$ [Fig. 3(h)] dominates $\Delta\rho_{\mathrm{tot}}(x,y,t)$ [Fig. 3(g)]. (iii) The most important feature is that the induced CDW charge density $\Delta\rho_{\mathrm{CDW}}(x,y)$ [Fig. 3(i)] is opposite to the original $\rho_{\mathrm{CDW}}(x,y,t_{0})$ [Fig. 3(c)], indicating a $20\%$ decrease in the excitonic order.

We analyze the integrated charges $Q_{i}(t)=\int|\rho_{i}(\mathbf{r},t)|d\mathbf{r}$ and $C_{i}(t)=\int|\rho_{i}(\mathbf{r},t)-\rho_{i}(\mathbf{r},t_{0})|d\mathbf{r}$ ( $i=\mathrm{tot}$ , $1\times 1$ , $\mathrm{CDW}$ ) for quantitative comparisons. The former characterizes the strength of bonding, while the latter denotes the weakening of bonding states, i.e. the strength of anti-bonding. As shown in Fig. 3(j), the percentage of the decrease $[Q_{i}(t_{f})-Q_{i}(t_{0})]/Q_{i}(t_{0})$ is 23.6%, 23.8% and 19.7% for $i=\mathrm{tot}$ , $1\times 1$ , $\mathrm{CDW}$ , respectively. Meanwhile, $C_{i}(t)$ increases, with a ratio $C_{\mathrm{tot}}(t_{f}):C_{\mathrm{1\times 1}}(t_{f}):C_{\mathrm{CDW}}(t_{f})=1:0.98:0.20$ , slightly different from the ratio of the initial bonding charges $Q_{\mathrm{tot}}(t_{0}):Q_{\mathrm{1\times 1}}(t_{0}):Q_{\mathrm{CDW}}(t_{0})=1:0.98:0.27$ . Thus, laser-induced bonding-antibonding transfer is nearly homogeneous, lowering both the $1\times 1$ and CDW order proportionally. Since the decrease in $|Q_{\mathrm{1\times 1},\mathrm{bond}}|$ affects all chemical bonds homogeneously, as an overall effect, the decrease in $|Q_{\mathrm{tot},\mathrm{bond}}|$ lowers the stability of CDW.

PLD dynamics. The laser-induced CDW instability would further trigger dynamic changes in PLD, which are simulated from first principles. Among the atoms in the layer of $2\times 2$ cell, only one Ti atom and one Se atom are nontrivial symmetry-inequivalent atoms, as shown in Fig. 4(a). We present the time-dependent displacement $\mathbf{d}_{i}(t)$ ( $i=\mathrm{Ti},\mathrm{Se}$ ) by the projection of $\mathbf{d}_{i}(t)$ on the PLD direction $\delta_{i}(t)=\mathbf{d}_{i}(t)\cdot\mathbf{d}_{i}(t_{0})/|\mathbf{d}_{i}(t_{0})|$ and the angle $\theta(t)$ between $\mathbf{d}_{i}(t)$ and $\mathbf{d}_{i}(t_{0})$ [Fig. 4(b)]. We find that laser triggers a set of movements $\Delta\mathbf{d}_{i}(t)=\mathbf{d}_{i}(t)-\mathbf{d}_{i}(t_{0})$ which are exactly antiparallel to $\mathbf{d}_{i}(t_{0})$ , noticing that $\delta_{i}$ decreases with all $\cos\theta\approx\pm 1$ as shown in Fig. 4(d-e). Besides the in-plane movements, the overall bonding-antibonding charge transition introduces an out-of-plane breathing mode between the Se layer and the Ti layer $z_{\mathrm{Se}}(t)$ [Fig. 4(f)]. This is the established out-of-plane $A_{1g}$ mode with the periodicity of 175 fs (5.7 THz) when $I=I_{0}/4$ , which is slightly smaller than the experimental value 5.9 THz Snow2003 ; Wang2018g . Higher laser fluence $I=I_{0}$ further decreases the frequency to 5.1 THz, due to the weakened Ti-Se bonds by photocarriers.

We analyze the dynamical potential energy surface (PES) related to the PLD dynamics. To focus on the in-plane PLD movements other than the out-of-plane $A_{1g}$ mode, we plot the minus of the kinetic energies of these two movements $-E_{\mathrm{kin}}$ and $-E_{\mathrm{kout}}$ , respectively. Since the total energy $E_{\mathrm{tot}}(t)=E_{\mathrm{p}}(t)+E_{\mathrm{kin}}(t)+E_{\mathrm{kout}}(t)$ is conserved, where $E_{\mathrm{p}}(t)$ is the potential energy, we use $-E_{\mathrm{kin}}(t)$ to characterize the dynamical PES of the PLD. As shown in Fig. 4(g), the dynamical potential energy decreases due to the in-plane PLD movements, together with the energy oscillations attributed to the out-of-plane $A_{1g}$ mode. This indicates that the laser induced PLD movements stabilize the system. Within the simulation time, $\mathbf{d}_{i}(t)$ has been completely inverted $\mathbf{d}_{i}(t)=-\mathbf{d}_{i}(t_{0})$ at around $300$ fs, lowering the dynamical PES and creating a new quasi-equilibrium state. We note that the inversion of PLD produces a degenerate CDW state which is symmetry-equivalent to the original state. As shown in Fig. 4(c), we generate the dynamical PES by plotting $-E_{\mathrm{kin}}(t)$ as a function of $\delta(t)$ . This implies a scenario of ultrafast PLD switch, which includes three consecutive steps: i) In the first 20 fs, the laser pulses pump the electrons from the bonding state to the anti-bonding state, which leads to the change of the potential energy surface.

The calculated double-well PES yields a vibration frequency of 4.3 THz, which is consistent with the oscillation period of the amplitude mode of about 220-240 fs observed in dynamic simulations (Fig. 4d). This frequency is comparable to the experimental value of 3.4 THz observed for the $A_{1g}$ CDW amplitude mode Snow2003 ; Wang2018g . The difference between theory and experiment is attributed to possible mixing of other phonon modes and inadequate accuracy of semilocal functionals to treat low-energy phonons. The softening of the double-well potential is evidenced by the fact that the CDW mode oscillations disappear when the laser fluence increases from $I_{0}/4$ to $I_{0}$ .

ii) The lattice transforms into a structure with the opposite PLD along the non-equilibrium TD-PES during 20–300 fs. iii) The photoexcited system might relax into the ground state after the recombination of photocarriers. This process is beyond our simulations and requires multichannel and large-scale modeling to properly account for the decoherence Floss2017 ; Yamada2018 ; Floss2019 and dissipation effects Tully1990 ; Gossel2018 . Here, we complete the story by briefly discussing the long-time behaviors: in the experiments, the ionic movements in different unit cells have different phases due to the finite thickness of the material, the inhomogeneous spatial distribution of the laser spot, as well as thermal fluctuations. Via phonon-phonon scattering, the equilibrium temperature gradually forms at the timescale of picoseconds. The thermally driven phase transitions may occur through an increasing population of phonons and changes in the atomic potential via anharmonic phonon interactions. Thus, the nonthermal CDW dynamics observed here is separated from thermal transitions at different time scales.

Dynamic interplay. The highly directional PLD movements are clearly not the reason but the consequence of the reduced excitonic order. In contrast to the complete PLD inversion for the $I=I_{0}$ case, $\mathbf{d}_{i}(t)$ only oscillates around the original value when laser intensity decreases to $I_{0}/4$ . The formation/decay of the CDW/PLD is electronically initialized, which supports the excitonic mechanism. This intensity dependence is also experimentally observed Mohr-Vorobeva2011 ; Porer2014 : the $A_{1g}$ CDW mode is softened as the laser fluence increases, which serves as a precursor of the PLD melting. Meanwhile, there are extra factors in the following dynamics. Assuming that the excitonic pairing is the only driving force, an inversion of electronic order $\rho_{\mathrm{CDW}}(\mathbf{r},t_{f})=-\rho_{\mathrm{CDW}}(\mathbf{r},t_{0})$ is required to reverse the PLD $\mathbf{d}_{i}(t)=-\mathbf{d}_{i}(t_{0})$ . Considering $\rho_{\mathrm{CDW}}$ only decreases by about 20%, the PLD should be only perturbed.

We design another numerical experiment by introducing a thermostat and removing a small fraction of kinetic energy at the rate of 0.01 eV atom*-1* ps*-1*. We note that the kinetic energy $E_{\mathrm{k}}(t)$ increases by 0.66 eV atom*-1* ps*-1* within 60 fs, as shown in Fig. 4(g). Thus the thermostat, affecting only 1.4% of the kinetic energy change, is expected to slightly delay the dynamics of PLD $\delta(t)$ . However, despite that the dissipation in the kinetic energy is small ( $1.4\%$ ), the PLD dynamics is completely quenched (Fig. 4, dotted lines). It strongly implies that the increase in $E_{\mathrm{k}}(t)$ is nonlinear, with a lower rate at the beginning. The fragile initial movement of $\delta(t)$ will be completely quenched with a small dissipation rate. Since the PLD dynamics is sensitively dependent on its own trajectory, we infer that there exists a self-amplified electron-phonon mechanism in the PLD dynamics initiated by exciton binding.

Accordingly, we demonstrate the impact of ionic movements on electronic structure through analyzing TD-EBS along $\Gamma$ -M in Fig. 5(a-f). The laser pulse first shifts the EBS downwards and enlarges the bandgap; then VBs are raised by $\sim$ 0.18 eV while CBs barely change during 20-80 fs, in excellent agreement with tr-ARPES data. It leads to the vanishing band gap and forms hole and electron pockets, boosting the relaxation of hole carriers towards the maximum of VBs (VBM). The carrier relaxation lowers the total energy of the system. In comparison, in the quenched case where ions barely move [Fig. 5(h-l)], the dynamic electronic evolution – the gap closing and carrier relaxation – are completely absent, implying that the measured changes in EBS are highly correlated with ionic movements. We thus identify a self-amplified process from these observations: the decrease in the PLD $\rightarrow$ the decrease in the bandgap $\rightarrow$ the decrease in the energy of carriers $\rightarrow$ the further decrease in PLD. We expect such a dynamic interplay between electronic order and PLD, ignored in previous studies, would sustain in a variety of CDW materials such as NbSe2, TaS2, LaTe3, etc Hellgren2017 ; Zhang2018CDWTaS2 . The results highlight the entangled dynamics between different degree of freedoms in quantum materials.

**Discussion

To complete our story, we note the following facts towards building a unified picture for experimental observations. i) The effective lattice temperature is relatively low during the whole laser-induced dynamics. The maximum kinetic energy of ions is 0.06 eV atom-1 in a temporary period (corresponding to a transient effective temperature of $450$ K). Although the low kinetic energy of ions may not be directly connected to nonthermal processes, we believe the processes discussed above are primarily nonthermal since they are short-lived ( $<$ 300 fs) and the excited-state carriers are clearly present (Fig. 5). In this case the structural soft-mode potential is relaxed or even driven to a single well-potential by excitation of the electronic degree of freedom. This is also consistent with the experimental observations Mohr-Vorobeva2011 and previous studies Lian2019Silicon ; Zhang2018CDWTaS2 ; Lian2016NTM . ii) The CDW order is disturbed instantly by the laser pulse within the timescale of 10 fs, since it is originated from homogeneous electronic transitions. No buildup time in electron-hole screening is needed. It shows that, the $\sim$ 100 fs shift in tr-ARPES signals Hellmann2012 is not caused by the buildup of electron-hole screening, but by the EPC-assisted PLD dynamics $\mathbf{d}_{i}(t)$ . Apparently, the phonon response within 100 fs is contradictory to the well accepted timescale of phonon dynamics $\tau_{\mathrm{ph}}>1$ ps. However, we note that the photo-induced PLD dynamics are coherently driven not by the phonon-phonon scattering, but directly by the electronic transitions. Thus, the response shares the same timescale as the electron-phonon scattering of $\sim$ 100 fs. iii) With the laser intensity used in our simulations, the laser pulses induce the PLD dynamics but barely change the intrinsic phonon modes (e.g. out-of-plane $A_{1g}$ ). This explains the sustainable characteristic peaks of PLD with the melted exciton order as previously observed in Ref. Porer2014 . iv) There are open discussions in literature on whether excitonic interaction or the Jahn-Teller distortion is dominant for the formation of CDW in TiSe2. Our simulations seem to indicate the two effects could coexist in a short period perturbed by photoexcitation. From Fig. 3 and 4 it is shown that while there is no significant lattice distortion (thus no changes in Jahn-Teller interaction), the electronic order is suppressed by 20% in the first 20 fs, beaconing excitonic interactions. On the other hand, the following self-amplifying process suggests dynamic Jahn-Teller lattice distortion is crucial for the subsequent decay/formation of CDW state.

We demonstrate an entangled exciton-phonon mechanism for charge ordering in 1T-TiSe2 by ab initio ultrafast rt-TDDFT simulations of photoexcitation dynamics. We show that laser pulses knock down the CDW order via inducing a homogeneous bonding-antibonding electronic transition. The weakened CDW order then triggers ionic movements exactly anti-parallel to the original PLD, but cannot solely drive the observed PLD inversion. Instead, assisted by the EPC, a self-amplification mechanism drives the PLD dynamics initiated by laser excitation, in excellent agreement with the tr-ARPES measurements. We propose that both the excitonic pairing and the electron-phonon coupling contribute to the CDW/PLD formation, but in a different timescale: the excitonic pairing initializes the formation of CDW, while the electron-phonon coupling promotes the following dynamics through a self-amplification process.

**Method

Following our previous efforts in developing time dependent ab initio package (TDAP) Meng2008a ; Lian2018MultiK ; Lian2018AdvTheo , we implemented TDDFT algorithm in the plane-wave code Quantum Espresso Giannozzi2009 . The details are described in Supplementary Note 1. We used the projector augmented-waves method (PAW) Blochl1994 and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation (XC) functional Perdew1996 in both DFT and TDDFT calculations. Pseudopotentials were generated using pslibrary DalCorso2014 . The plane-wave energy cutoff was set to 55 Ry. Brillouin zone (BZ) was sampled using Monkhorst-Pack scheme Monkhorst1976 with a $6\times 6\times 3$ $k$ -point mesh. Band unfolding techniques were utilized to generate the effective band structure (EBS) Ku2010 ; Popescu2012 with a modified version of BandUP code Medeiros2014 ; Medeiros2015 ; Lian2017SiUnfold . Onsite Coulomb repulsion $U=3.5$ eV was added to the Ti atom to reproduce the experimental band structure, while we used $U=0$ in the dynamic TDDFT calculations and structural optimization. The electron timestep $\delta t$ is $4\times 10^{-4}$ a.u.= $0.1935$ attosecond and the ion timestep $\Delta t$ is $0.194$ fs. Laser pulses with a wavelength $\lambda=800$ nm, width $\sigma=12$ fs and fluence $I_{0}=2.1$ mJ cm-2 centered at $t=10$ fs is utilized. The linear-response TDDFT calculations were carried out with the YAMBO package Marini2009 . This setup, which is similar to that used in experiment, reproduces the measured number of excited carriers in TiSe2(Supplementary Fig. 2).

**Acknowledgment

**We acknowledge partial financial supports from MOST (grants 2016YFA0300902 and 2015CB921001), NSFC (grants 91850120, 11774396, 11934003), and CAS (XDB07030100).

**Author Contributions

**S.M. conceived and directed the research. C.L. performed the calculations and analyzed the data. C.L., S.-J.Z., S.-Q.H., M.-X.G. and S.M. participated in the discussion. C.L. and S.M. wrote the paper with contributions from all the other authors.

Competing interests: The authors declare no competing interests.

I Algorithms of TDDFT in Planewave basis

I.1 Planewave and Adiabatic Basis

Following our previous scenario in the Time dependent ab initio package (TDAP) Meng2008a ; Lian2018MultiK ; Lian2018AdvTheo , we implement the TDDFT algorithm in PW basis. The time-dependent Kohn-Sham (TDKS) equation at time $t$ in PW basis $\{\mathbf{G}\}$ reads Runge1984 :

[TABLE]

where $\psi_{\gamma\mathbf{k}}(\mathbf{G},t)$ is TDKS orbital, $\gamma$ denotes the band index, $\mathbf{k}$ is the reciprocal momentum index. $\mathcal{H}_{\mathbf{k}}(t)$ is the Hamiltonian expanded with plane-wave basis, with matrix element

[TABLE]

where $T_{\mathbf{k}}(\mathbf{G},\mathbf{G^{\prime}})=\frac{\hbar^{2}}{2m}|\mathbf{k}+\mathbf{G}+\mathbf{A}(t)|^{2}\delta_{\mathbf{G},\mathbf{G^{\prime}}}$ is the kinetic term, $\mathbf{A}$ is the velocity gauge potential Bertsch2000 ; Yabana2006 :

[TABLE]

where $\mathbf{E}$ is electric field. $V(\mathbf{G},\mathbf{G^{\prime}})$ is the potential term calculated within the corresponding module Quantum Espresso, including such as ion-electron potential, Hartree potential and exchange-correlation potential.

There are $N_{b}$ independent TDKS equations for each $\mathbf{k}$ index, and thus $N_{b}N_{\mathbf{k}}$ equations in total, where $N_{b}$ is the number of bands, typically, $N_{\mathbf{k}}$ is the number of $\mathbf{k}$ points. Typically, the orders of $N_{\mathbf{k}}$ and $N_{b}$ are both $\sim 10^{2}$ . The number of $\{\mathbf{G}\}$ , $N_{G}$ , is usually $\sim 10^{4}$ . Thus, if the number of time steps is $N_{t}$ , the complexity of solving Eq. (1) is $O(N_{t}N_{b}N_{\mathbf{k}}N^{2}_{G})$ .

We express the TDKS orbitals $\psi_{\gamma\mathbf{k}}(\mathbf{G},t)$ using adiabatic basis $\{\phi_{\gamma\mathbf{k}}(\mathbf{G},t_{1})\}$

[TABLE]

where $c_{i\gamma,\mathbf{k}}(t)$ is the TD coefficients and the adiabatic basis $\{\phi_{i\mathbf{k}}\}(t_{1})$ are solved by diagonalizing the Hamiltonian

[TABLE]

where the $\epsilon_{i\mathbf{k}}$ is the eigenvalue. The initial condition is chosen as

[TABLE]

where $\delta_{i\gamma}$ is Kronecker delta. The TDKS equations are interpreted with the evolution of coefficient matrix

[TABLE]

where

[TABLE]

and $h_{ij,\mathbf{k}}=\braket{\phi_{i\mathbf{k}}(t_{1})}{\mathcal{H}_{\mathbf{k}}(t)}{\phi_{j\mathbf{k}}(t_{1})}$ . To distinguish, we use $\mathcal{A}$ and $A$ to represent the matrix $A$ (e.g. Hamiltonian) in the PW basis and the adiabatic basis, respectively.

Note that, no approximation has yet been introduced in the derivation. The complexity of solving the Eq. 5 is the same as solving the Eq. 1, while it is less time consuming to solve Eq. 7 noticing that the dimension of ${C}_{\mathbf{k}}(t)$ and ${H}_{\mathbf{k}}(t)$ is $N_{b}\times N_{b}$ and $N_{b}\ll N_{G}$ . The computation will be accelerated if the time of solving of Eq. 5 is minimized.

I.2 Efficient Evolution in Adiabatic Basis

An efficient scheme is raised by Wang et. al. Wang2015a : Hamiltonian on adiabatic basis $h_{ij,\mathbf{k}}(t)=\braket{\phi_{i\mathbf{k}}(t_{1})}{\mathcal{H}_{\mathbf{k}}(t)}{\phi_{j\mathbf{k}}(t_{1})}$ changes approximately linearly within $[t_{1},t_{2}]$ with $t_{2}=t_{1}+\Delta t$ fs and $\Delta t\sim 0.2$ fs Ren2013 :

[TABLE]

The matrix element $h_{ij,\mathbf{k}}$ is evaluated as

[TABLE]

and

[TABLE]

where

[TABLE]

$\ket{\phi_{l,\mathbf{k}}(t_{2})}$ and $\epsilon_{l\mathbf{k}}(t_{2})$ are the adiabatic basis and the eigenvalue at time $t_{2}$ , which are solved in a self-consistent process, see detail at Sec. I.3 .

Propagator operator $U_{\mathbf{k}}(t_{2},t_{1})$ is calculated with the knowledge of $H_{\mathbf{k}}(t)$ in the Crank-Nicholson scheme

[TABLE]

where $dt=\Delta t/N_{t}$ is the integration time step. Since $dtH\ll 1$ is needed to satisfy the condition of $\exp(-i\hbar dtH)=1-{i\hbar}Hdt$ , we choose $dt\sim 0.1$ attosecond, considering H is about $10^{2}$ eV. The coefficients $c_{i\gamma,\mathbf{k}}$ are thus propagated as

[TABLE]

Thus, we finish the evolution of TDKS orbital from $t_{1}$ to $t_{2}$ as

[TABLE]

where

[TABLE]

Charge density $\rho(t_{2})$ can be calculated with $c^{\prime}_{l\gamma,\mathbf{k}}(t_{2})$ and $\ket{\phi_{l,\mathbf{k}}(t_{2})}$ as

[TABLE]

where

[TABLE]

is the population of the adiabatic states.

I.3 Computational Flowchart

We organized the equations in a flowchart, as shown in Fig. 1. Step (A) (B) (C) (H) (I) (K) are computed with the Quantum Espresso modules. As mentioned above, the algorithm requires foreknowledge of $\rho(\mathbf{G},t_{2})$ in the calculation of Hamiltonian $H_{\mathbf{k}}(t)$ . It is satisfied with a self-consistent process:

An initial guess of $\rho(\mathbf{G},t_{2})$ is built using extrapolation from previous step. 2. 2.

Hamiltonian $H_{\mathbf{k}}(t)$ and propagator $U_{\mathbf{k}}(t_{2},t_{1})$ are built from $\rho(\mathbf{G},t_{2})$ using Eq. 10 and Eq. 14, respectively. 3. 3.

An new $\rho(\mathbf{G},t_{2})$ is calculated with new population propagated from Eq. 15. 4. 4.

The first three steps are repeated until the new and old $\rho(\mathbf{G},t_{2})$ are the same.

I.4 Forces and dynamics

Once the self-consistency in charge density evolution is satisfied, post-processing including the calculation of total energy, Hellmann-Feynman forces, and the ionic trajectory are invoked. For instance, the forces acting on the ions can be calculated through

[TABLE]

where $\mathbf{R}_{I}$ and $\mathbf{F}_{\mathbf{R}_{I}}$ are the position and force of $I$ th ion.

With $\mathbf{R}_{I}$ and $\mathbf{F}_{\mathbf{R}_{I}}$ , we utilize the Ehrenfest theorem for evolving ions according to the equation of motion

[TABLE]

where $M_{I}$ is the mass of $I$ th ion. The velocity $v_{I}(t)=d\mathbf{R}_{I}/dt$ and the temperature $T(t)=\sum_{I}^{N_{I}}M_{I}v^{2}_{I}(t)/2N_{I}$ are also calculated, where $N_{I}$ is the total number of ions.

Besides the conventional NVE ensemble, additional thermostats, such as Nosé-Hoover Nose1984 ; Hoover1985 and Berendsen Berendsen1984 is considered to simulate different environmental conditions. In the damped MD simulations, we utilize a simple velocity-rescaling thermostat. The ionic velocities are rescaled at each time step as

[TABLE]

where $\mathbf{v}_{I}^{\prime}(t)$ and $T^{\prime}(t)=T(t)-\Delta T$ are the rescaled velocity and temperature, respectively. The decreasing rate $\Delta T=0.01$ eV atom*-1* ps*-1* is used in the simulations.

Combining the TDKS equation and the Ehrenfest theorem, the many-body electron-electron interaction and the ionic movement under the excited-state TDKS wavefunction evolution are described in an ab initio way. We expect that the electron-electron interactions at the adiabatic XC level and electron-phonon scatterings within the mean-field average trajectory are present in these simulations. The excess electronic energy could dissipate into available phonon modes via electron-phonon coupling or to low-energy electrons via electron-electron scattering, resulting in carrier thermalization and cooling effect.

I.5 Projector Augmented Wave Method

To expand the TDKS orbital, the adiabatic basis $\ket{\phi_{i\mathbf{k}}(\mathbf{G},t_{1})}$ should be orthonormal. Eigenstates in all electron method and norm-conserving pseudopotential satisfy the orthonormal requirement naturally. However, using projector augmented-waves method (PAW) Blochl1994 bring additional core functions, which lead to nonorthogonal eigenstates. The formalism of DFT quantities such as total energy and forces using PAW method are described in Blöchl’s original paper Blochl1994 . Here, we only list the related changes in TDDFT evolution.

Using PAW method, all the coefficients $C^{\prime}_{\mathbf{k}}(t)$ are related with the pseudo-eigenstates $\ket{\tilde{\phi}_{i\mathbf{k}}(\mathbf{G},t_{1})}$

[TABLE]

where $\hat{T}$ is the transform operator. The only change caused by PAW method is Eq (13). Using PAW, eq (13) becomes

[TABLE]

where

[TABLE]

$s_{ij,\mathbf{k}}=\braket{\beta_{i\mathbf{k}}(\mathbf{G})}{\phi_{j\mathbf{k}}(\mathbf{G})},$ and $\beta_{i\mathbf{k}}(\mathbf{G})$ is the Kleinman-Bylander projectors Giannozzi2009 .

I.6 Time Dependent Band Unfolding

The CDW phase of TiSe2 is a $2\times 2$ cell of the normal phase. The energy bands are folded from $1\times 1$ Brillouin zone (BZ) to the $2\times 2$ BZ. In contrast, ARPES measurements still span over the $1\times 1$ BZ. To bridge the gap between DFT bands and measured ARPES spectra, the band unfolding technique is used to calculate the effective band structure (EBS) of the supercell (SC). Expanding the adiabatic basis $\ket{\phi_{i,\mathbf{k}}(t)}$ of $2\times 2$ SC in the adiabatic basis $\ket{\Phi_{I,\mathbf{K}}(t)}$ of primitive $1\times 1$ cell (PC), we get

[TABLE]

where $\mathbf{K}=\mathbf{k}+\mathbf{B}$ and $\mathbf{B}$ is the reciprocal basis vector of SC. The spectral function is the EBS along the $\mathbf{K}$ path in PCBZ Popescu2012 ; Medeiros2014 :

[TABLE]

where $E$ is the energy and

[TABLE]

We can introduce an extra weight function $w(i,\mathbf{k})$

[TABLE]

The choice of $w_{i,\mathbf{k}}$ is arbitrary Lian2017SiUnfold . Here, we use the population of the adiabatic states as the

[TABLE]

to reproduce the intensity in ARPES spectra.

II Fluence Dependence of Photocarrier Density

We investigate the photo-carrier density as a function of time and laser fluence $n(t,I)$ . Here, the carrier density is calculated as $n(t)=\frac{1}{2}\sum_{i,\mathbf{k}}|q_{i\mathbf{k}}(t)-q_{i\mathbf{k}}(t=0)|$ . We note that $n(t)$ accurately describes the number of excited carriers after the laser field ends $A(t)=0$ . Otherwise, a gauge independent projection on $\psi_{k^{\prime}}(t)$ with $k^{\prime}=k-A(t)/c$ instead of $k^{\prime}=k$ can be used Otobe2008 ; Otobe2009 ; Otobe2016b ; Yabana2012 ; Shinohara2010 ; Shinohara2010a ; Shinohara2012 ; Sato2015 ; Sato2015a . These two projections are identical when the laser field ends $A(t)=0$ .

As shown in Fig. 2, the $n(t=20~{}\mathrm{fs},I)$ are approximately proportional to the laser fluences. After comparing with the experimental measurements Porer2014 , we find the $n(I)$ relations are consistent at both the low fluence area $I<0.04$ mJ cm*-2* and the high fluence area $I>0.22$ mJ cm*-2*, while a clear derivation from linearity is observed in experiments. Exciton correlation drives the superlinear feature but is underestimated in DFT/TDDFT calculations (Fig. 3). Nevertheless, since the difference is minimized when $I>0.22$ mJ cm*-2*, we focus on the experimental phenomena in high laser fluence region, to reproduce the similar excitation states for direct comparisons.

Consistent with experimental observations Mohr-Vorobeva2011 , we note that the PLD dynamics are not sensitive to the photon energies. This is because that the band gap of CDW 1T-TiSe2 (0.18 eV) is smaller than the photon energy in most commonly used laser sources ( $\sim$ 1 eV). Besides, the sub-picosecond laser pulses utilized in the experiments and our simulations would bring up significant broadening in photon energy as well as multi-photon absorption processes.

Bibliography116

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) Morosan, E. et al. Superconductivity in Cu x Ti Se 2 . Nat. Phys. 2 , 544–550 (2006).
2(2) Li, G. et al. Semimetal-to-Semimetal Charge Density Wave Transition in 1T-Ti Se 2 . Phys. Rev. Lett. 99 , 027404 (2007).
3(3) Kusmartseva, A. F., Sipos, B., Berger, H., Forró, L. & Tutiš, E. Pressure Induced Superconductivity in Pristine 1T-Ti Se 2 . Phys. Rev. Lett. 103 , 236401 (2009).
4(4) Yao, Q. et al. Charge transfer effects in naturally occurring van der waals heterostructures (Pb Se) 1.16 (Ti Se 2 ) m (m=1,2). Phys. Rev. Lett. 120 , 106401 (2018).
5(5) Wei, M. J. et al. Manipulating charge density wave order in monolayer 1T-Ti Se 2 by strain and charge doping: A first-principles investigation. Phys. Rev. B 96 , 165404 (2017).
6(6) Kogar, A. et al. Observation of a charge density wave incommensuration near the superconducting dome in Cu x Ti Se 2 . Phys. Rev. Lett. 118 , 027002 (2017).
7(7) Yan, S. et al. Influence of domain walls in the incommensurate charge density wave state of Cu intercalated 1T-Ti Se 2 . Phys. Rev. Lett. 118 , 106405 (2017).
8(8) Medvecká, Z. et al. Observation of a transverse Meissner effect in Cu x Ti Se 2 single crystals. Phys. Rev. B 93 , 100501 (2016).