Real-time decay of a highly excited charge carrier in the one-dimensional Holstein model

We study the real-time dynamics of a highly excited charge carrier coupled to quantum phonons via a Holstein-type electron-phonon coupling. This is a prototypical example for the nonequilibrium dynamics in an interacting many-body system where excess energy is transferred from electronic to phononic degrees of freedom. We use diagonalization in a limited functional space (LFS) to study the nonequilibrium dynamics on a finite one-dimensional chain. This method agrees with exact diagonalization and the time-evolving block-decimation method, in both the relaxation regime and the long-time stationary state, and among these three methods it is the most efficient and versatile one for this problem. We perform a comprehensive analysis of the time evolution by calculating the electron, phonon and electron-phonon coupling energies, and the electronic momentum distribution function. The numerical results are compared to analytical solutions for short times, for a small hopping amplitude and for a weak electron-phonon coupling. In the latter case, the relaxation dynamics obtained from the Boltzmann equation agrees very well with the LFS data. We also study the time dependence of the eigenstates of the single-site reduced density matrix, which defines the so-called optimal phonon modes. We discuss their structure in nonequilibrium and the distribution of their weights. Our analysis shows that the structure of optimal phonon modes contains very useful information for the interpretation of the numerical data.

Figure 1

Sketch of the initial condition, Eq. (2), and the time evolution: we start from the state with one electron at momentum $k=\pi$ and no phonon. The total energy is thus equal to the initial electronic kinetic energy $E_{\mathrm{kin}}(t=0)={\epsilon }_{\pi }=2t_0$, where ${\epsilon }_k=-2t_{0}cosk$. Due to the coupling to phonons, the electron loses energy by exciting phonons while moving through the lattice, which also results in the redistribution of its momentum occupations.

Figure 2

Sketch of the states in the limited functional space that are generated by Eq. (16) during the first three iterations (we take $M=1$ for simplicity). The ring represents a lattice with $L=6$ sites and periodic boundary conditions. Each state in the figure represents the parent state of a translationally invariant state. The site of the electron is represented by an arrow for clarity, while the local phonons are represented by harmonic oscillators. Starting from the free electron state $(n_{\mathrm{h}}=0)$, a single phonon (since $M=1)$ is created on the site of the electron. In the next generation, the $n_{\mathrm{h}}=1$ state is taken as the starting state: the electron hops to the left or right of the phonon (leading to two new states) and a state with an additional phonon is created. Next, the second generation $(n_{\mathrm{h}}=2)$ states serve as starting states and the new unique states of the third generation $(n_{\mathrm{h}}=3)$ are added to the basis. This procedure is repeated until the desired number of generations is set up (when a state is generated that already exists in an earlier generation, it is omitted from the generation).

Figure 3

Time evolution of the phonon energy $E_{\mathrm{ph}}\left(t\right)$ for $\lambda =0.5$ and ${\omega }_0=t_0$. (a) Comparison of the numerical methods for $L=8$. The solid line represents results obtained by exact diagonalization (ED) in the complete Hilbert space, using the global phonon cutoff $N_{\max }=12$. The dashed line represents results obtained by diagonalization in a limited functional space (LFS) by using the parameters $N_{\mathrm{h}}=17$ and $M=1$ [see Eq. (16)]. (b) Convergence of $E_{\mathrm{ph}}\left(t\right)$ by using LFS for $L=12$. For the parameters of the basis generation, we keep $M=1$ while $N_{\mathrm{h}}$ is varied. The legends display the corresponding basis dimension $D$.

Figure 4

Comparison of the phonon energy $E_{\mathrm{ph}}\left(t\right)$ obtained by diagonalization in LFS and the TEBD method. We set $\lambda =0.5$, ${\omega }_0=t_0$, and $L=12$. Open boundary conditions are used in this figure (we use the free-electron eigenstate with the highest energy as the initial state).

Figure 5

Short-time dynamics of the phonon energy $E_{\mathrm{ph}}\left(t\right)$. The bold-dashed line represents the results from perturbation theory, given by Eq. (20). All other curves are numerical results using LFS at ${\omega }_0=t_0$ and different $\lambda$ (we use $L=12$ for $\lambda =0.5$ and 1, and $L=8$ for $\lambda =2$ and $4.5)$. Time is measured in units of $1/\gamma$.

Figure 6

Time evolution of the phonon energy $E_{\mathrm{ph}}\left(t\right)$ in the weak-coupling antiadiabatic limit for $L=12$. We set $t_0/{\omega }_0=0.1$ and $\gamma /{\omega }_0=0.01$, which corresponds to $\lambda =0.5\times {10}^{-3}$. The solid line is the numerical result using LFS, while the dashed line is the result from perturbation theory, given by Eq. (37).

Figure 7

Relaxation of the electronic kinetic energy $E_{\mathrm{kin}}\left(t\right)$ at ${\omega }_0=t_0$ in the weak-coupling regime ${(\gamma /t_0)}^2=0.4$, which corresponds to $\lambda =0.2$. The solid line is the numerical result using LFS (for $L=12)$, the dashed line is the result from the numerically integrated Boltzmann equations (50). Time is measured in units of $\tau$, defined in Eq. (53).

Figure 8

Relaxation of the electronic kinetic energy $E_{\mathrm{kin}}\left(t\right)$ in the weak-coupling regime. Results are shown for a fixed $\gamma$ and $t_0$ [setting ${(\gamma /t_0)}^2=0.4$] and different values of ${\omega }_0/t_0=0.5,0.8,1.0,1.2,1.5$, which correspond to $\lambda =0.4,0.25,0.2,0.17,0.13$, respectively. We use LFS for $L=12$. The dashed line is the function $E_{\mathrm{kin}}\left(t\right)/t_0=2[1-2(t/\tau )]$ with the relaxation time $\tau$ defined in Eq. (53).

Figure 9

The electronic kinetic energy $E_{\mathrm{kin}}\left(t\right)$ in the strong-coupling antiadiabatic regime at $\gamma /{\omega }_0=1$ and $L=12$. The dashed line is the result from perturbation theory, given by Eq. (68), while the solid and thin dashed line represent the TEBD result for $t_0/{\omega }_0={10}^{-3}$ and ${10}^{-1}$ (corresponding to $\lambda =500$ and 5, respectively).

Figure 10

Time evolution for the two parts of the Hamiltonian: (a) the phonon energy divided by the bandwidth $E_{\mathrm{ph}}\left(t\right)/\left(4t_0\right)$ and (b) the electron kinetic energy $E_{\mathrm{kin}}\left(t\right)$, plotted for several values of the adiabaticity ratio ${\omega }_0/t_0$. We use LFS at $\lambda =0.2$ and $L=12$.

Figure 11

Density plot showing the electronic momentum distribution function $n_k\left(t\right)$ using LFS at ${\omega }_0=t_0$ and $L=12$. (a) $\lambda =0.5$ and (b) $\lambda =2$.

Figure 12

Time evolution for the three competing parts of the Hamiltonian: (a) phonon energy $E_{\mathrm{ph}}\left(t\right)$, (b) electron-phonon coupling energy $E_{\mathrm{coup}}\left(t\right)$, and (c) electronic kinetic energy $E_{\mathrm{kin}}\left(t\right)$. We set ${\omega }_0=t_0$ and show results using LFS for the three parameter values $\lambda =0.5$, 1 and 2 (we use $L=12$ for the first two systems and $L=8$ for the last one).

Figure 13

Standard deviation ${\sigma }_{\Delta E_{\mathrm{ph}}}$ [see Eq. (69)] of the phonon energy $E_{\mathrm{ph}}\left(t\right)$ in the stationary regime, as a function of the inverse lattice size. Circles represent results using LFS, while the dashed lines in (a) and (b) are the linear fitting functions $f\left(L\right)$ with one and two free parameters, respectively. (a) $\lambda =0.5$, where oscillations vanish in the thermodynamic limit (the two points for the smallest $L$ were not included in the fit). (b) $\lambda =2.0$, where oscillations persist for all $L$.

Figure 14

Time evolution of the phonon energy $E_{\mathrm{ph}}\left(t\right)$ in the strong-coupling regime. The solid lines represent numerical results using LFS for $t_0={\omega }_0$ and $L=8$, while the dashed lines represent the single-site dynamics $(t_0=0)$, given by Eq. (64). (a) $\lambda =2$, which corresponds to $g=2$. (b) $\lambda =4.5$, which corresponds to $g=3$.

Figure 15

Time-averaged values (t-AVG) of the three competing parts of the Hamiltonian in the stationary state as a function of electron-phonon coupling $\lambda$ at $L=12$ and ${\omega }_0=t_0$, obtained by using LFS. The values are compared to the polaron ground state (GS) energies $E_{\mathrm{kin}}$, $E_{\mathrm{coup}}$ and $E_{\mathrm{ph}}$ at $K=\pi$. (a) Main panel: $E_{\mathrm{kin}}+E_{\mathrm{coup}}$ vs $\lambda$ in the ground state (diamonds) compared to the corresponding time-averaged values in the stationary state (circles). The insets display the comparison of the individual terms. The data at $\lambda \to 0$ exhibit an $L$ dependence due to the large spatial extent of the polaron. (b) Phonon energy $E_{\mathrm{ph}}$ vs $\lambda$. Since the averages in the stationary state clearly exceed the ground-state values, we also compare to $E_{\mathrm{kin}}(t=0)-(E_{0,\mathrm{coup}}+E_{0,\mathrm{kin}})$, where $E_{\mathrm{kin}}(t=0)=E_{\mathrm{total}}$, and $E_{0,\mathrm{coup}}$, $E_{0,\mathrm{kin}}$ are measured in the ground state.

Figure 16

The electron kinetic energy $E_{\mathrm{kin}}\left(t\right)$ at ${\omega }_0/t_0=10$ for several values of the electron-phonon coupling $g=\gamma /{\omega }_0$. We use LFS and set $L=12$.

Figure 17

Time evolution of the von Neumann entropy $S_{\text{vN}}\left(t\right)$ for ${\omega }_0=t_0$ and different values of the electron-phonon coupling $\lambda =0.5$, 2, and $4.5$. The values at $t=0$ are determined by Eq. (78). The horizontal lines represent the corresponding ground-state (GS) values. We use LFS with $L=8$.

Figure 18

Weights ${\tilde{w}}_{\alpha }$ of the optimal modes in descending order for ${\omega }_0=t_0$ and different values of the electron-phonon coupling $\lambda =0.5$, 2, and $4.5$ (we use LFS with $L=12$ for $\lambda =0.5$ and $L=8$ for the latter two systems). Filled symbols represent weights in the stationary regime while open symbols represent the corresponding ground-state (GS) values. We choose a time $t_{\max }$ in the stationary regime at which the phonon energy $E_{\mathrm{ph}}\left(t\right)$ has a local maximum [c.f. Figs. 12 and 14], i.e., $t_{\max }{t}_0=15.8$ for $\lambda =0.5$ and $t_{\max }{t}_0=9.3$ for $\lambda =2$ and $4.5$.

Figure 19

Components ${\left|\langle n|\alpha \rangle \right|}^2$ of the four highest optimal modes $|\alpha \rangle$ in the phonon occupation number basis $\left\{|n\rangle \right\}$. Results are shown in the stationary regime using LFS at $\lambda =0.5$, ${\omega }_0=t_0$ and $L=12$ for two different times: at $t_{\min }$ (squares) the phonon energy $E_{\mathrm{ph}}\left(t\right)$ has a local minimum, while at $t_{\max }$ (circles) it has a local maximum [see the time evolution of $E_{\mathrm{ph}}\left(t\right)$ in Fig. 12]. ${\tilde{w}}_{\alpha }^{\min }$ and ${\tilde{w}}_{\alpha }^{\max }$ represent the optimal mode weights (79) at the given times.

Figure 20

Components ${\left|\langle n|\alpha \rangle \right|}^2$ of the four highest optimal modes $|\alpha \rangle$ in the phonon occupation number basis $\left\{|n\rangle \right\}$. Results are shown in the stationary regime using LFS at $\lambda =4.5$, ${\omega }_0=t_0$, and $L=8$ for two different times: at $t_{\min }$ (squares) the phonon energy $E_{\mathrm{ph}}\left(t\right)$ has a local minimum, while at $t_{\max }$, (circles) it has a local maximum [see the time evolution of $E_{\mathrm{ph}}\left(t\right)$ in Fig. 14]. ${\tilde{w}}_{\alpha }^{\min }$ and ${\tilde{w}}_{\alpha }^{\max }$ represent the optimal mode weights (79) at the given times. Dashed lines in (a) and (b) represent Poisson distributions (71).