Photoinduced charge-order melting dynamics in a one-dimensional interacting Holstein model

TL;DR

This study investigates the ultrafast dynamics of charge order melting in a one-dimensional fermion-phonon system after optical excitation, revealing different mechanisms depending on phonon frequency and confirming the soliton-based interpretation.

Contribution

It provides a detailed numerical analysis of photoinduced charge order melting in a one-dimensional Holstein model, highlighting the roles of quantum phonons and soliton dynamics under different conditions.

Findings

Charge order melting is explained by soliton propagation without phonons.

Large phonon frequency induces polaronic soliton-driven melting.

Small phonon frequency involves energy transfer and soliton band motion.

Abstract

Transient quantum dynamics in an interacting fermion-phonon system are investigated. In particular, a charge order (CO) melting after a short optical-pulse irradiation and roles of the quantum phonons on the transient dynamics are focused on. A spinless-fermion model in a one-dimensional chain coupled with local phonons is analyzed numerically. The infinite time-evolving block decimation algorithm is adopted as a reliable numerical method for one-dimensional quantum many-body systems. Numerical results for the photoinduced CO melting dynamics without phonons are well interpreted by the soliton picture for the CO domains. This interpretation is confirmed by the numerical simulation for an artificial local excitation and the classical soliton model. In the case of the large phonon frequency corresponding to the antiadiabatic condition, the CO melting is induced by propagations of the polaronic solitons with the renormalized soliton velocity. On the other hand, in the case of the small phonon frequency corresponding to the adiabatic condition, the first stage of the CO melting dynamics occurs due to the energy transfer from the fermionic to phononic systems, and the second stage is brought about by the soliton motions around the bottom of the soliton band. Present analyses provide a standard reference for the photoinduced CO melting dynamics in low-dimensional many-body quantum systems.

Equations

$$\mathcal{H}=$$

$$+V\sum _i\left(n_i-\frac{1}{2}\right)\left(n_{i+1}-\frac{1}{2}\right)$$

$$+{\omega }_0\sum _i{a}_i^{†}{a}_i-g\sum _i(a_i+a_i^{†})\left(n_i-\frac{1}{2}\right),$$

$$t_{i,i+1}\to t{e}^{\mathrm{i}A(\tau )},$$

$$A(\tau )=\frac{A_p}{\sqrt{2\pi }{\tau }_p}{e}^{{\tau }^2/(2{\tau }_p^2)}\cos ({\omega }_p\tau ),$$

$$|\Psi ⟩=\sum _{{\sigma }_1,{\sigma }_2,\dots }\mathrm{Tr}(B^{{\sigma }_1}{B}^{{\sigma }_2}\cdots )|{\sigma }_1,{\sigma }_2,\dots ⟩,$$

$$|{\Psi }_{\mathrm{GS}}⟩\propto \prod _{n=1}^N\exp (-\mathcal{H}\delta \tau )|{\Psi }_0⟩,$$

$$O=\frac{1}{N}\left|\sum _i(-1{)}^i\left(⟨n_i⟩-\frac{1}{2}\right)\right|,$$

$$q=\frac{1}{N}\left|\sum _i(-1{)}^i⟨a_i^{†}+a_i⟩\right|.$$

$$\Delta E_X(\tau )=⟨\Psi (\tau )|{\mathcal{H}}_X|\Psi (\tau )⟩-⟨{\Psi }_{\mathrm{GS}}|{\mathcal{H}}_X|{\Psi }_{\mathrm{GS}}⟩,$$

$${\chi }_{\mathcal{AA}}(\omega )=⟨{\Psi }_{\mathrm{GS}}\left|{\mathcal{A}}^{†}\delta (\omega -\mathcal{H}+E_0)\mathcal{A}\right|{\Psi }_{\mathrm{GS}}⟩.$$

$${\chi }_{\mathcal{AA}}(\omega )=\frac{1}{\pi \sqrt{1-{\omega }^2}}\left[g_0{\mu }_0+2\sum _{n=1}^{N_{\mathrm{ch}}-1}{g}_n{\mu }_n{T}_n(\omega )\right],$$

$$g_n=\frac{(N_{\mathrm{ch}}-n+1)\cos \frac{\pi n}{N_{\mathrm{ch}}+1}+\sin \frac{\pi n}{N_{\mathrm{ch}}+1}\cot \frac{\pi }{N_{\mathrm{ch}}+1}}{N_{\mathrm{ch}}+1}.$$

$${\chi }_{jj}(\omega )=⟨{\Psi }_{\mathrm{GS}}\left|j\delta (\omega -\mathcal{H}+E_0)j\right|{\Psi }_{\mathrm{GS}}⟩,$$

$$E(k,k^{\prime })={\epsilon }_k+{\epsilon }_{k^{\prime }},$$

$$N_p(\tau )=\frac{\Delta E(\tau )}{N{\omega }_p},$$

$$O(\tau )=O_0\exp (-\gamma \tau ),$$

$$t_{i,i+1}\to \{\begin{array}{ll}t{e}^{\mathrm{i}A(\tau )}& (i=\cdots -L,0,L,\cdots )\\ t& \text{other sites},\end{array}$$

$$O_i(\tau )=(-1{)}^i\left(⟨n_i⟩-\frac{1}{2}\right).$$

$$\gamma =-\frac{2v\ln (1-2L{n}_p)}{L},$$

$${\mathcal{H}}_{\mathrm{eff}}=-t\sum _{i=\mathrm{odd}}\left(b_i^{†}{b}_{i+2}+d_i^{†}{d}_{i+2}+H.c.\right)+\frac{V}{2},$$

$$v_k=|4t\sin 2k|.$$

$$v_{{\omega }_p}=\sqrt{(4t{)}^2-({\omega }_p-V{)}^2}.$$

$$v^{\ast }(\tau )=-\frac{1}{4N_p(\tau )}\frac{dO(\tau )}{d\tau },$$

$$\gamma /(4n_{p}t)\propto e^{-g^2/{\omega }_0^2},$$

Full text

I Introduction

Ultrafast dynamics in correlated electron systems induced by the intensive optical laser pulse have been widely accepted as an attractive theme in recent solid-state physics. nasu ; PIPT_JPSJ This issue is motivated not only from the fundamental research in highly nonequilibrium many-body systems but also from the potential applications to the efficient energy conversion and ultrafast optical switching. One of the prototypical targets for the nonequilibrium many-body systems is the Mott insulator and the carrier-doped Mott insulator. A variety of exotic nonequlibrium phenomena emerges owing to the spin and orbital degrees of freedom, han ; matsuda ; iwai ; okamoto ; novelli ; stojchevska ; cavalleri ; kubler ; perfetti and a number of theoretical studies has been also performed in the optical induced nonequilibrium state in the Mott insulators. takahashi ; oka ; freericks ; moritz ; meisner ; werner ; vidmar ; mierzejewski ; lu ; iyoda

Another promising target for the photoinduced nonequilibrium dynamics in correlated electron systems is the charge ordered (CO) state, in which the averaged electron density per site is 0.5 in most of the cases. The CO states, in which the translational symmetry of the electronic charge distribution is broken, are ubiquitously observed in a wide class of solids. On the verge of CO, several rich phenomena emerge, such as superconductivity, metal-insulator transition, ferroelectricity, colossal magneto-resistance and so on. Since the CO system is susceptible against the external stimuli due to its small energy scale in comparison with the Mott insulator, this system is a plausible candidate for the nonequilibrium electron dynamics research. kiryukhin ; poli ; beaud ; iwai2 ; ichikawa ; onda ; iwano09 ; rincon ; hashimoto14 ; hashimoto15 Actually, a variety of CO compounds have been investigated by the ultrafast optical pump-probe measurements as well as the several time-resolved experimental techniques. In contrast to the Mott insulators, the spin degree of freedom hides behind. Instead, the lattice degree of freedom plays crucial roles for the equilibrium and nonequilibrium states as follows: i) The non-uniform electronic charge distribution and the alternate lattice distortion are stabilized cooperatively. ii) Different energy and time scales in the electronic and phononic systems provide a rich variety of the transient dynamics. iii) The energy transfer possibly occurs from the highly excited electronic to phononic systems.

During the several CO systems, one dimensional CO insulator is a simple example for study of the nonequilibrium dynamics induced by the optical excitation. The organic molecular solid (TMTTF)${}_{2}X$ ($X$: monovalent anion) is a prototypical one-dimensional CO compound, in which the number of holes per the molecular orbital is 0.5. ishiguro ; jacobsen ; balicas ; dressel ; monceau ; naitoh The quantum fluctuation in a low dimensional lattice promotes a strong competition between the long-range CO insulator and the correlated metal, and governs the photoinduced transient dynamics between the two states. Quantum nature of lattice vibrations is also crucial on the energy transfer and the relaxation dynamics in the nonequilibrium transient states. werner2 ; sayyad ; murakami ; golez ; lee ; dorfner ; filippis From the theoretical viewpoint, several numerical algorithms have been applied for one-dimensional quantum systems. maeshima05 ; maeshima06 ; maeshima07 ; iwano11 ; hassanieh Reliable calculations in one-dimensional systems may provide good references for the highly nonequilibrium transient dynamics in correlated electrons, in which the quantum effect, the many body effect and the transient characters are treated properly on an equal footing.

In this paper, we address the photoinduced CO melting dynamics in a one-dimensional correlated electron system. In particular, we focus on roles of the quantum fluctuations of electrons and phonons on the transient dynamics. We analyze the one-dimensional spinless fermion model at half filling by using the infinite time-evolving block decimation (iTEBD) method. This method is well accepted as a reliable numerical algorithm for the one-dimensional interacting quantum model with high accuracy. In the model without the phonon degree of freedom, the CO melting is caused by soliton propagations, corresponding to the domain-wall motions in the CO state. In the case of the finite fermion-phonon coupling, the transient dynamics depend qualitatively on the adiabatic parameter $\alpha\equiv\omega_{0}/t$, where $\omega_{0}$ and $t$ are the phonon frequency and the fermion hopping integral, respectively. In the case of large $\alpha$, the CO melting occurs in the same way to the case without the phonon degree of freedom and is induced by the soliton motions in which the velocity is renormalized by the polaronic effect. On the other hand, in the case of small $\alpha$, the CO melting occurs by the energy transfer to the phononic system at the first stage, and then by the soliton motion around the bottom of the soliton band. The present comprehensive analyses provide a standard reference for the photoinduced CO melting dynamics in low dimensional systems.

In Sect. II, the fermion-phonon coupled model is introduced, and the numerical method based on the iTEBD algorithm is briefly explained. In Sect. IIIA, the numerical results for the ground state before photoexcitation are presented. The transient dynamics without the fermion-phonon interaction are introduced in Sect IIIB, and the results with the phonon degree of freedom are presented in Sect. IIIC. Section IV is devoted to discussion and summary. Detailed formulation for the classical soliton model is given in Appendix A.

II Model and Method

The model in which we analyze the photoinduced CO dynamics is the interacting spinless-fermion model on a one-dimensional chain coupled with the local phonons. The Hamiltonian is defined as

[TABLE]

where $c_{i}(c_{i}^{\dagger})$ is the annihilation (creation) operator of a spinless-fermion at site $i$, $n_{i}(\equiv c_{i}^{\dagger}c_{i})$ is the fermion number operator, and $a_{i}(a_{i}^{\dagger})$ is the annihilation (creation) operator of a local phonon at site $i$. The first term $({\cal H}_{t})$ and the second term $({\cal H}_{V})$ represent the fermion hopping and the Coulomb interaction between the fermions in the nearest neighboring sites, respectively. The third term $({\cal H}_{\rm ph})$ and the last term $({\cal H}_{\rm ep})$ represent the phonon energy and the Holstein-type fermion-phonon interaction with a coupling constant $g$, respectively. As already introduced in Sect. I, $\alpha=\omega_{0}/t$ is the adiabatic parameter, and $\alpha\ll 1$ and $\alpha\gg 1$ are termed the adiabatic and antiadiabatic limits, respectively. We also introduce the dimension-less fermion-phonon coupling constant $\lambda\equiv g/\omega_{0}$.

The light pulse applied along a chain direction is introduced as the Peierls phase in the hopping integral as

[TABLE]

where the lattice constant, the light velocity and the Planck constant are taken to be one. The vector potential for the pump pulse at time $\tau$ is given as

[TABLE]

where $A_{p}$, $\omega_{p}$, and $\tau_{p}$ are the amplitude, frequency and pulse width of the vector potential, respectively. The center of the vector potential is set to be $\tau=0$.

The wave functions for the fermions and phonons in the ground state and the time-evolved states are calculated using the iTEBD method. vidal07 ; dorfner ; takayoshi ; ono1 ; ono2 The wave function is represented as the matrix-product state form as

[TABLE]

where $\sigma_{i}$ describes the local states for the fermion and phonon at site $i$, and $B^{\sigma_{i}}$ is a matrix with the dimension $\chi$. The ground state wave function is obtained by the imaginary-time evolution as

[TABLE]

where $|\Psi_{0}\rangle$ is the initial wave function, and $\delta\tau$ is small difference of the imaginary time. In the calculation for the time-evolved states, $|\Psi(\tau)\rangle$, we use Eq. (5), in which the imaginary time is replaced by the real time as $\delta\tau\rightarrow i\delta\tau$. To calculate the exponential factor in Eq. (5), we adopt the second and fourth order Suzuki-Torotter decompositions for the imaginary-time and real-time equations, respectively. hatano It is assumed that the wave functions are invariant under shifts of the two sites in most of the calculations, and have a periodicity under shifts of $64$ sites in the calculations for the local excitations introduced later. In order to monitor the CO state, we calculate the staggered CO parameter defined as

[TABLE]

and the alternate ion displacement defined as

[TABLE]

Here, $\langle\cdots\rangle$ represents the expectation value with respect to $|\Psi_{\rm GS}\rangle$ in the ground state, and that with respect to $|\Psi(\tau)\rangle$ for the time-evolved state. For each term of the Hamiltonian in Eq. (1), we introduce the transient energy at $\tau$ measured from the initial energy defined as

[TABLE]

where $X$ identifies the term of the Hamiltonian.

Accuracies in the numerical calculations are checked by calculating the CO parameter and the ground state energies for various parameter sets. The present model at $g=0$, i.e. the interacting spinless fermion model without phonon, is equivalent to the $S=1/2$ XXZ model in a one-dimensional lattice, and $V=2t$ corresponds to the Heisenberg model. mila ; ejima The exact ground state is obtained by the Bethe ansatz: samaj a Luttinger-liquid state as the unique ground state in $V/t<2$ and the two-fold degenerated CO states in $V/t>2$. Differences of the calculated ground-state energy per site $(E_{0})$ from the exact results $(E_{\rm exa})$, i.e. $|(E_{0}-E_{\rm exa})/E_{\rm exa}|$, are of the order of $10^{-14}$ ($10^{-9}$) at $V/t=6$ ($V/t=3$) for $\chi=64$. Time dependences of the CO parameter at $g=0$ and $V/t=6$ are shown in Fig. 1(a) for various values of $\chi$. It is confirmed that the data for $\chi\geq 256$ converge up to $\tau t\simeq 30$. In the calculations for $g\neq 0$, the maximum phonon number per site $(N_{\rm ph})$ is truncated. In Fig. 1(b), $\Delta E_{\rm ph}$ are shown for several values of $N_{\rm ph}$ and $\chi$. The calculated results are almost converged up to around $\tau t=25$ at $N_{p}=5$ and $\chi=256$. By taking into account of the above results, in most of the numerical calculations for the transient states, we adopt $\chi=512$ for $g/\omega_{0}<0.5$, $\chi=384$ for $0.5\leq g/\omega_{0}<1$, and $\chi=256$ for $1\leq g/\omega_{0}$, and $N_{\rm ph}=3-8$.

The excitation spectra are calculated based on the matrix-product state formalism, in which finite size clusters of $L_{c}$ sites with the open boundary condition are adopted. Here, the spectral functions are expanded using the Chebyshev polynomial. holzner We define the excitation spectra induced by the operator ${\cal A}$ as

[TABLE]

The spectral function is expanded as

[TABLE]

where $T_{n}(\omega)$ is the $n$-th Chebyshev polynomial, $\mu_{n}=\langle\Psi_{\rm GS}|{\cal A}^{\dagger}T_{n}({\cal H}){\cal A}|\Psi_{\rm GS}\rangle$ is the Chebyshev moment, and $g_{n}$ is the correction factor defined by

[TABLE]

Based on the recursion formula for the Chebyshev vectors, $|t_{n+1}\rangle=2{\cal H}|t_{n}\rangle-|t_{n-1}\rangle$ with $|t_{0}\rangle={\cal A}|\Psi_{\rm GS}\rangle$ and $|t_{1}\rangle={\cal H}|t_{0}\rangle$, these are calculated in the matrix-product state formalism.

III Results

III.1 Ground state

Before showing the photoinduced dynamics in the CO systems, we introduce the ground state properties as the initial state. In Figs. 2(a) and (b), the CO parameter defined in Eq. (6) and the lattice distortions in Eq. (7) for several values of $g$ are presented, respectively. Without the fermion-phonon coupling ($g=0$), the CO parameter tends to vanish at $V/t=2$ with decreasing $V/t$ as predicted by the exact solution. With increasing $g$, the CO phase is strengthened. The alternate ion distortion occurs in the CO phase, and the saturated values of $q$ increase with $g$.

The optical absorption spectra are shown in Fig. 2(c) for several values of $\alpha$ as well as $g=0$, and in Fig. 2(d) for several values of $\lambda$. We define the optical absorption spectra as

[TABLE]

where $j=\sum_{i}({\rm i}tc_{i}^{\dagger}c_{i+1}+{\rm H.c.})$ is the current operator along the chain. The elementary excitations in the ordered phase of the one-dimensional XXZ model is known as the domain wall motions. The optical absorption spectra, in which the net momentum transfer is zero, are ascribed to the kink and antikink pair with the opposite momenta. ishimura ; kanamori The energy for a non-interacting pair of the kink with momentum $k$ and the antikink with momentum $k^{\prime}$ is given by

[TABLE]

with $\varepsilon_{k}=-2t\cos(2k)+V/2$. This formula indicates that the optical absorption spectra distribute in $V-4t\leq\omega\leq V+4t$, which is consistent with the main part of the optical spectra shown in Fig. 2(c) where $V/t=6$. We have confirmed that a weak hump structure below $V-4t$, such as a small peak at around $\omega/t=1.9$, reduces with increasing the cluster size, and this is an artifact due to the size effect. On the other hand, the main structure does not depend sensitively on the cluster size. The spectra in the antiadiabatic case ($\alpha=3.2$) are similar to the spectra without the fermion-phonon interaction, and some fine structures appear in the adiabatic case ($\alpha=0.4$). With increasing $\lambda$, multiple peak structures appear above the main absorption edge at around $V+n\omega_{0}$ with an integer number $n$ as shown in Fig. 2(d). In the calculations of the photoinduced dynamics shown in the following section, the pump-photon frequency is chosen around the energy in which $\chi_{jj}(\omega)$ is finite, implying that the pump photon yields resonant excitations.

III.2 Transient state without the fermion-phonon coupling

In this subsection, the photoinduced CO melting in the case without the fermion-phonon interaction are introduced. In Fig. 3(a), the time profiles of $\Delta E_{V}(\tau)$ after the excitation are shown for several values of $\omega_{p}$. We introduce the effective absorbed photon number per site defined as $n_{p}=N_{p}(\tau\gg\tau_{p})$ with

[TABLE]

where $\Delta E(\tau)=\langle\Psi(\tau)|{\cal H}|\Psi(\tau)\rangle-E_{0}$ is a change in the total energy. In several cases of $\omega_{p}$ adopted in Fig 3(a), $n_{p}$ are unified to be approximately 0.0005 by adjusting $A_{p}$. After introducing the pump pulse, $\Delta E_{V}(\tau)$ rapidly increases and settles at almost time-independent values; the photoexcited state is stable after turning off the pump pulse. In the case of small $\omega_{p}$, a oscillatory behavior appears and is damped at around $\tau=10/t$. As shown in Fig. 3(b), $n_{p}$ increases quadratically as a function of $A_{p}$. This is interpreted that the number of the domain walls in the CO states, i.e. the kinks and antikinks, increases with $A_{p}$.

The time profiles of the CO parameter are shown in Fig. 4(a). Reduction in $O(\tau)$ occurs monotonically for all values of $A_{p}$. In contrast to the time dependence of $\Delta E_{V}(\tau)$ shown in Fig. 3(a), the CO melting occurs in a long time scale after turning off the pump pulse. In the case of large $A_{p}$, $O(\tau)$ reaches zero, implying almost complete melting of the initial CO. As shown in the semi-logarithmic plot in Fig. 4(b), the time dependences of $O(\tau)$ follows an exponential function as

[TABLE]

where $O_{0}$ is the CO parameter in the initial state, and $\gamma$ is the damping factor. We have confirmed that the exponential dependence of $O(\tau)$ is commonly observed in wide parameter ranges of $\omega_{p}$ and $V$. The damping factor is calculated for several values of $A_{p}$ at fixed $\omega_{p}(=4t)$, and is plotted as a function of $n_{p}$ in Fig. 5(a). It increases linearly for small $n_{p}$, and deviates upward with increasing $n_{p}$. As shown in Fig. 5(b), the $\omega_{p}$ dependence of $\gamma$ shows a dome shape. These two behaviors of $\gamma$ are well explained by the results of the soliton picture plotted by bold lines in Figs. 5(a) and (b), as explained later.

The CO melting dynamics are analyzed by the following artificial electronic excitation. The vector potential for the pump pulse is applied at a specific site instead of the uniform vector potential introduced in Eq. (2). This is termed “the local excitation” from now on and is defined as

[TABLE]

where $L$ is a periodicity of the sites in which the excitations are induced, and the time dependence of $A(\tau)$ is given in Eq. (3). Responses to the local excitation are monitored by the CO parameter defined at each site termed “the local CO parameter” defined as

[TABLE]

An intensity plot of $O_{i}(\tau)$ induced by the local excitations in a space-time plane is plotted in Fig. 6(a). A periodicity of the local excitations is chosen to be $L=64$. After introducing the pump pulse, a CO domain characterized by the reduced $O_{i}(\tau)$ is generated, and the boundaries between the domains propagate almost linearly with time. At $\tau\approx 9/t$, the domain wall collides at $i=\pm 32$ with other domain walls which are generated at $i=\pm L$. After the collision, $O_{i}$ reduces furthermore. Next collisions occur at approximately $\tau=18/t$ and $i=0$. A velocity of the domain wall is about 4 sites per the unit time of $1/t$. After the collisions, both the widths and velocities of the domain walls are almost unchanged. From these characteristics, the domain wall dynamics in the photoinduced CO state are identified as the soliton motions. As plotted in Fig. 6(b), the CO parameter $O(\tau)$ calculated from the data of $O_{i}(\tau)$ shows the exponential decay in the same way to the results in Fig. 4. We conclude that the soliton motions are responsible for the photoinduced CO melting.

Next, we derive the damping factor based on the soliton picture. We set up the model in which the photogeneraged solitons are treated as classical particles. A photoinduced fermionic excitation is simulated by a stochastic generation of solitons. We assume that the soliton velocity is constant $v$, and the collisions between the solitons are elastic. Details are presented in Appendix A. We obtain an expression of $\gamma$ as

[TABLE]

which is approximately $\gamma=4vn_{p}$ in the limit of $n_{p}\ll 1$. The soliton velocity in Eq. (18) is evaluated from the effective Hamiltonian derived from Eq. (1) in the limit of $V\gg t$. The number of the soliton-antisoliton pair is assumed to be one. The effective Hamiltonian is given as kanamori

[TABLE]

where $b_{i}^{\dagger}$ ($b_{i}$) and $d_{i}^{\dagger}$ ($d_{i}$) are the creation (annihilation) operators of the hard-core bosons for the soliton and antisoliton, respectively. A summation with respect to $i$ takes for odd integers. The interaction between the soliton and antisoliton given by the second oder of $t$ is not included. The eigenstates are classified by the momenta of the soliton $k$ and antisoliton $k^{\prime}$, and the eigenenergy coincides to Eq. (13). The group velocity of soliton is obtained as

[TABLE]

A factor 2 in the sine function is attributed to the fact that the fermion hopping between the nearest neighbor site corresponds to the soliton hopping between the next-nearest neighbor sites. In the case where the photon momenta is zero and the resonant excitation is induced, i.e. $\omega_{p}=E(k,-k)$, the velocity is represented as a function of $\omega_{p}$ as

[TABLE]

The maximum group velocity is given as $4t$, being consistent with the results in Fig. 6(a).

The analytical expression of the damping factor in Eq. (18) with Eq. (21) is compared with the numerical data. Here, $L$ in Eq. (18) is treated as a fitting parameter and is determined to be approximately $L=6.5$. It is shown that the $n_{p}$ dependence of $\gamma$ plotted in Fig. 5(a) well reproduces the calculated data. As shown in Fig. 5(b), $\gamma$ as a function of $\omega_{p}$ shows a dome shape structure which takes its maximum at approximately $\omega_{p}=V$. This behavior is confirmed by the numerical data for several values of $V$. The analytical formulation little deviates from the numerical data in small $V$, in which the higher-order quantum fluctuation neglected in Eq. (18) might be important.

III.3 Effect of the fermion-phonon coupling

In this subsection, the numerical results in the case of $g\neq 0$ are presented. Figure 7(a) shows the time profiles of $O(\tau)$ for several values of $\alpha$ at $\lambda=0.25$. For comparison, we also present the results at $g=0$, in which $O(\tau)$ decreases exponentially as a function of time as explained previously. With introducing the fermion-phonon coupling, the reduction of $O(\tau)$ is suppressed. The exponential decay of $O(\tau)$ remains in the case of large phonon frequency $(\alpha=3.2)$, but the time profiles of $O(\tau)$ deviates from the exponential decay with decreasing $\alpha$. To clarify the change in the CO melting dynamics by the fermion-phonon coupling, the concept of the soliton velocity introduced previously is extended to the transient state. From Eqs. (15) and (18) under the conditions of $\gamma\tau\ll 1$ and $n_{p}\ll 1$, we have $O(\tau)\sim e^{-\gamma\tau}\sim 1-4vn_{p}\tau$. Then, the transient soliton velocity at time $\tau$ is introduced as

[TABLE]

where $N_{p}(\tau)$ is defined in Eq. (14). The calculated results of $v^{\ast}(\tau)$ are shown in Fig. 7(c). In the case of the large $\alpha$ ($\alpha=3.2$), $v^{\ast}(\tau)$ is almost independent of time after turning off the pump pulse in the same way of $g=0$. The oscillatory behaviors seen in $\tau\geq 10$ are suppressed with increasing $\chi$, and are attributed to the numerical artifact. On the other hand, for the small $\alpha$ $(\alpha=0.4,{\rm and\ }0.8)$, $v^{\ast}(\tau)$ changes remarkably in a short time scale and tends to be saturated in a long time scale. This characteristic time scale in the adiabatic case ($\alpha=0.4$ and 0.8) will be discussed later. The time profiles of $O(\tau)$ for several values of $\lambda$ are presented in Figs. 7(b) and (d), where $\alpha$ is fixed to be 3.2, i.e. the antiadiabatic condition. The exponential time dependence of $O(\tau)$ and the almost constant $v^{\ast}(\tau)$ are observed for $0\leq\lambda\leq 3.2$, and $v^{\ast}(\tau)$ monotonically decreases with increasing $\lambda$. Detailed CO melting dynamics in the cases of the adiabatic and antiadiabatic conditions are explained in the following two subsections.

III.3.1 Antiadiabatic case (\alpha\ \raise 1.50696pt\hbox{>}\kern-7.5pt\lower 2.15277pt\hbox{\sim}\ 1)

First, we focus on the dynamics in the antiadiabatic case at $\alpha=3.2(>1)$, where the CO parameter decays exponentially. The time profiles of $\Delta E_{t}(\tau)$ and $\Delta E_{\rm ph}(\tau)+\Delta E_{\rm ep}(\tau)$ are shown for several values of $\lambda$ in Figs. 8(a) and (b), respectively. Except for small dip and hump structure just after turning off the pump pulse at $\lambda=2.4\ {\rm and}\ 3.2$, both of the energies do not show remarkable time dependences. With increase $\lambda$, the saturated value of $\Delta E_{t}(\tau)$ decreases monotonically. The intensity plots in the space-time plane for $O_{i}(\tau)$ induced by the local excitation are presented in Fig. 9. For all cases of $\lambda$, the soliton motions are stable and their velocities are almost unchanged. With increasing $\lambda$, the velocity decreases, being consistent with the results in Fig. 7(d).

From the numerical data shown in Fig. 7, the normalized damping rate $\gamma/(4n_{p}t)$, corresponding to the soliton velocity, is plotted as a function of $\lambda$ in Fig. 10. The velocity decreases monotonically with increasing $\lambda$. The bold line in the figure represents the equation given by

[TABLE]

which is expected from the renormalized carrier hopping in the polaronic state based on the Lang-Firsov theory. lang The analytical formula fits well the numerical data even in the region of small $\lambda$, implying that the coherent soliton motion associated with the ion displacement occurs from weak to strong coupling regimes. This is termed that polaronic-soliton picture from now on, and is consistent with the fact that $\Delta E_{\rm ph}(\tau)$ does not show remarkable change after turning off the pump pulse (not shown). We conclude that, in the antiadiabatic case, the CO melting dynamics is caused by the soliton motion renormalized by the polaron formation.

III.3.2 Adiabatic case (\alpha\ \raise 1.50696pt\hbox{<}\kern-7.5pt\lower 2.15277pt\hbox{\sim}\ 1)

Next, we focus on the CO melting dynamics in the adiabatic case in $\alpha<1$, where the time dependence of $O(\tau)$ deviates from the exponential decay as shown in Fig. 7(a). The time profiles of $\Delta E_{t}(\tau)$ and $\Delta E_{\rm ph}(\tau)+\Delta E_{\rm ep}(\tau)$ are shown for several values of $\alpha$ in Figs. 8(c) and (d), respectively. It is confirmed that $\Delta E_{V}(\tau)$ suddenly increases by pulse irradiation, and does not show remarkable time dependence after turning off the pump pulse (not shown). The time profiles of $\Delta E_{t}(\tau)$ are clearly different from those in the antiadiabatic case presented in Fig. 8(a); $\Delta E_{t}(\tau)$ gradually decreases after turning off the pump pulse and are saturated at approximately $\tau t=10$ and 15 for $\alpha=0.8$ and 0.4, respectively. Counter behaviors are observed in $\Delta E_{\rm ph}(\tau)+\Delta E_{\rm ep}(\tau)$, which increases gradually with time. We also confirm that $\Delta E_{\rm ep}(\tau)$ is one order smaller than $\Delta E_{\rm ph}(\tau)$. The above results are interpreted that the excess energy by the pump excitation is transfered from the fermionic to phononic systems, and the charge distributions are compatible to the ion displacements.

In order to analyze the energy transfer, the time profiles of $\Delta E_{\rm ph}(\tau)+\Delta E_{\rm ep}(\tau)$ are fitted by a function $C(e^{-\tau/\tau_{r}}-1)$ where $\tau_{r}$ is the energy relaxation time, and $C$ is a numerical constant. In Fig. 11, the energy relaxation times calculated for several values of $\alpha(<1)$ and $\lambda$ are plotted as a functions of $g^{2}/(\omega_{0}t)$. Almost linear dependence is observed until $g^{2}/(\omega_{0}t)\sim 0.16$, implying the energy relaxation in accordance with the Fermi golden rule through the fermion-phonon coupling. werner2 ; sayyad ; murakami ; golez Above $g^{2}/(\omega_{0}t)=0.16$, the width of the pump pulse chosen to be $\tau_{p}=1/t$ is comparable to the energy relaxation time and the relaxation occurs before turning off the pump pulse. In Figs. 11(b) and (c), respectively, intensity maps of $O_{i}(\tau)$ in the space-time plane after the local excitation are presented for $\alpha=0.4$ and 0.8 at $\lambda=0.25$. It is shown that the wave fronts of the domain walls become broad and the averaged velocities decrease with time. This is clearly in contrast to the case of large $\alpha$ (see Fig. 9), in which the soliton velocities are almost independent of time. The reduction of the soliton velocity is understood in the soliton model in Eq. (19): as shown in Eq. (20), the soliton velocity is minimum (maximum) at the bottom (middle) in the soliton band. Thus, the energy relaxation from the fermionic to phononic systems leads to the reduction of the velocity, implying a suppression of the CO melting.

Finally, characteristics in the damping factor is summarized; $\gamma/(4n_{p}t)$ are plotted as a function of the inverse of $\alpha$ in Fig. 12. The damping factors are obtained from the time profiles of $O(\tau)$ at a time when $v^{\ast}(\tau)$ is fully saturated. With increasing $\alpha$, the numerical data approaches to a red broken line, which represents the damping factor without the fermion-phonon coupling multiplied by the polaronic factor $e^{-(g/\omega)^{2}}$. Below $\alpha\sim 1$, reduction of the calculated results from the red line is remarkable and is attributable to the energy relaxation from the fermionic to phononic systems.

IV Discussion and Summary

We present in this paper a comprehensive analyses of the photoinduced melting dynamics of the one-dimensional CO state interacting with the local phonons. By using the iTEBD method, we obtain the reliable numerical data for the time dependences of the CO parameter, the electronic and phononic energy and so on. We also analyze the CO melting dynamics in the system in which the artificial local photoexcitation is introduced. As well known, the excited state in the one-dimensional CO state without the electron-phonon coupling is described by the soliton and antisoliton picture. This is directly confirmed by the calculations of the local CO parameter induced by the local excitation shown in Fig. 6; the domain-wall motions collapse the CO structure and give rise to the exponential time dependence of the CO parameter as $O(\tau)\propto e^{-\gamma\tau}$. The characteristic CO melting induced by the local excitation is reproduced by the classical soliton picture presented in Appendix A, in which the solitons are treated as classical particles. We have the expression for the damping constant $\gamma$ as shown in Eq. (18), being proportional to the soliton velocity in the limit of the weak photoexcitation.

The roles of the quantum phonons on the CO melting dynamics are characterized by the dimension-less adiabatic parameter $\alpha=\omega_{0}/t$. We perform the numerical analyses in a wide parameter region of $\alpha$. In the case of the antiadiabatic region (\alpha\ \raise 1.50696pt\hbox{>}\kern-7.5pt\lower 2.15277pt\hbox{\sim}\ 1), the CO melting dynamics occurs by the polaronic soliton picture. The stable propagations of the domain walls are survived even by taking account of the fermion-phonon coupling, but the soliton velocity is renormalized. The damping constant $\gamma$ decreases with increasing the fermion-phonon coupling. This reduction of $\gamma$ in a wide region of the coupling constant ($0<g/\omega_{0}<1.25$) is well fitted by the expression for the polaronic state deduced by Lang and Firsov, lang which is known to be derived in the strong coupling regime. This implies that the solitons and ionic displacements propagate coherently even in the case of small $g$. This might be attributed to the fact that the ions are distorted cooperatively in the CO ground state even for small g, in which incoherent motions of electrons and ions produce a energy loss in macroscopic order.

On the contrary, in the adiabatic case (\alpha\ \raise 1.50696pt\hbox{<}\kern-7.5pt\lower 2.15277pt\hbox{\sim}\ 1), the CO melting in the early stage of the photoexcited transient state is not induced by the propagations of the solitions but their relaxations. The photogenerated solitons have a band energy of $\varepsilon_{k}=-2t\cos 2k$ and a velocity $v=4t\sin 2k$, where the interaction with phonons and that between solitons are neglected. Trough the fermion-phonon coupling, the kinetic energy of the solitons are dissipated to the phononic system, and the soliton velocity is reduced. The solitons are settled down around the bottom of the band, at which these cannot be relaxed furthermore by emitting the quantum phonons. After the relaxation, the CO melting is induced by the soliton propagation, in the same way with the antiadiabatic case. This scenario is confirmed by the calculation of the transient electron energy and phonon energy, in which the relaxation time is scaled by $g^{2}/(\omega_{0}t)$. The change of the time dependence of the CO parameter, that is, the non-exponential decay in the early stage and the exponential decay after the energy relaxation, is understood in this scenario.

We thank M. Naka, H. Seo and S. Iwai for their helpful discussions. This work was supported by MEXT KAKENHI Grant No. 26287070, 15H02100 and 17H02916. Some of the numerical calculations were performed using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo.

Appendix A A classical soliton model

In this Appendix, we introduce a classical model for the photoinduced soliton and antisoliton in order to derive the analytical expression of the damping factor given in Eq. (18). We introduce a classical particle model in a one-dimensional chain, in which hard core particles (antiparticles) correspond to the solitons (anti-solitons) in the CO state. We assume that the (anti)particle velocity is constant $v$, and the elastic collision occurs between the (anti)particles. The pump photon irradiation is simulated by a stochastic generation of the particle-antiparticle pairs in the bonds connecting the nearest neighboring sites with a periodicity $L$ and probability $n$. We assume the $\pi$ phase shift of CO and a step-function like CO parameter across the (anti)soliton as shown in Fig. 13(a). The local CO parameter is calculated in this model and the intensity map of $O_{i}(\tau)$ in the space-time map is shown in Fig. 13(b). The result well reproduces the numerical results shown in Fig. 6(a). At time $\tau_{m}=Lm/(2v)$ when the $m$-time collision occurs as indicated by broken lines in Fig. 13(b), the local order parameters are spatially uniform. We confirm at least up to $m=5$ that the amplitude of the CO parameter at the time $\tau_{m}$ are given approximately by $O(\tau_{m})=(1-2n)^{m}$. This result is interpreted as follows. Let us consider the state at $\tau_{1}=L/(2v)$, until when the (anti)solitons propagate independently. The local CO parameters at the sites, where the (anti)soliton pass and does not pass over, are $1$ and $-1$, respectively. Since the probability for the photogeneration of solitons are $n$, the CO parameter at this time is obtained by the average calculation as $O(\tau_{1})=1\times(1-n)+(-1)\times n=1-2n$. At the time $\tau_{2}=2L/(2v)$, the second collision occurs between the soliton and antisoliton which are generated at $i=0$ and $L$. The four photoinduced states are possible: the (anti)soliton is generated at both the sites of $i=0$ and $L$, at either $i=0$ site or $i=L$ site, and neither the sites. Since these probabilities are $n^{2}$, $n(1-n)$, $n(1-n)$ and $(1-n)^{2}$, respectively, the CO parameter is obtained by the average as $O(\tau_{2})=1\times n^{2}+(-1)\times n(1-n)+(-1)\times n(1-n)+1\times(1-n)^{2}=(1-2n)^{2}$. The order parameters at $\tau_{m}$ with $m\geq 3$ are obtained in the same way. By using the relation $n=n_{p}L$ where $n_{p}$ is the effective photon number introduced in Sec. III.2, the CO parameter is given as $O(\tau_{m})=(1-2n)^{m}=\exp[2v\tau_{m}\ln(1-2Ln_{p})/L]$. This expression can be extended to the continuous time. By comparing this result with Eq. (15), the damping factor for the CO parameter is identified as $\gamma=-2vL^{-1}\ln(1-2Ln_{p})$ which is Eq. (18).