Charge-density-wave melting in the one-dimensional Holstein model

We study the Holstein model of spinless fermions, which at half filling exhibits a quantum phase transition from a metallic Tomonaga-Luttinger liquid phase to an insulating charge-density-wave (CDW) phase at a critical electron-phonon coupling strength. In our work, we focus on the real-time evolution starting from two different types of initial states that are CDW ordered: (i) ideal CDW states with and without additional phonons in the system and (ii) correlated ground states in the CDW phase. We identify the mechanism for CDW melting in the ensuing real-time dynamics and show that it strongly depends on the type of initial state. We focus on the far-from-equilibrium regime and emphasize the role of electron-phonon coupling rather than dominant electronic correlations, thus complementing a previous study of photoinduced CDW melting [H. Hashimoto and S. Ishihara, Phys. Rev. B 96, 035154 (2017)]. The numerical simulations are performed by means of matrix-product-state based methods with a local basis optimization (LBO). Within these techniques, one rotates the local (bosonic) Hilbert spaces adaptively into an optimized basis that can then be truncated while still maintaining a high precision. In this work, we extend the time-evolving block decimation (TEBD) algorithm with LBO, previously applied to single-polaron dynamics, to a half-filled system. We demonstrate that in some parameter regimes, a conventional TEBD method without LBO would fail. Furthermore, we introduce and use a ground-state density-matrix renormalization group method for electron-phonon systems using local basis optimization. In our examples, we account for up to $M_{\mathrm{ph}}=40$ bare phonons per site by working with $O\left(10\right)$ optimal phonon modes.

Figure 1

(a) Sketch of the different terms in the Holstein model Eq. (1). The fermions can hop from site to site with an amplitude $t_0$. If a fermion is on a particular site it can create or destroy phonon excitations at that site with a coupling strength $\gamma$. Every phononic excitation costs an energy ${\omega }_0$. (b) Sketch of the phase diagram of the half-filled Holstein model [38, 39]. As the electron-phonon coupling $\gamma$ increases there is a phase transition from a metallic Tomonaga-Luttinger liquid phase (TLL) to an insulating charge-density-wave phase (CDW). The arrows represent the different quenches that we will investigate in Sec. 4c, i.e., frequency and coupling quench (FQ and CQ, respectively).

Figure 2

Different steps of the DMRG with the subspace-expansion and local basis optimization update. (a) Shift of the focus to the basis transformation tensor and optimization. (b) Shift of the focus back to the site tensor and truncation. (c) Transformation of the local part of the Hamiltonian matrix product operators into the optimized basis (see also [80]).

Figure 3

Relative error $\mathrm{\Delta }\epsilon$ of the ground-state energy obtained with the $\mathrm{DMRG}3\mathrm{S}+\mathrm{LBO}$ algorithm calculated for $L=4,{\omega }_0/t_0=2$, and $\gamma /t_0=4$ (crosses: Data calculated with $d_o=5$, diamonds: $d_o=8$). Red symbols were calculated with $M_{\mathrm{ph}}=20$ and blue symbols with $M_{\mathrm{ph}}=30$. The discarded weight for the bond dimension is ${10}^{-10}$. The thin black dotted line marks $\mathrm{\Delta }\epsilon ={10}^{-10}$.

Figure 4

(a) General structure of a TEBD algorithm [89, 95]. (b) Different steps in the application of a single local time-evolution operator in the TEBD with local basis optimization algorithm [77].

Figure 5

Comparison between TEBD-LBO data (open black symbols) and Lanczos time-evolution data (small blue symbols) of the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). Calculations are done for system size $L=4$, phonon frequency ${\omega }_0/t_0=2$, and different coupling strengths $\gamma /t_0=1,3,4$ (squares, diamonds, and circles, respectively). In the TEBD-LBO time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=10,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fourth data point that was computed in TEBD-LBO and every twentieth data point from the Lanczos time evolution.

Figure 6

Sketch of the initial states: (a) $|\mathrm{BCDW}\rangle$ Eq. (32) and (b) $|\mathrm{DCDW}\rangle$ Eq. (34).

Figure 7

Time evolution of (a) the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ and (b) the phonon number per fermion $N_{\mathrm{ph}}/N$ when starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). The small black dots in panel (a) are exact analytical results for $\gamma =0$ in the thermodynamic limit [97]. The dashed horizontal lines in panel (b) represent the phonon number in the ground state at the respective parameters. Simulations are performed for $L=13, {\omega }_0/t_0=2$ and different coupling strengths $\gamma /t_0=1,3,4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=10,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 8

Time evolution of (a) the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ and (b) the phonon number per fermion $N_{\mathrm{ph}}/N$ when starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). Simulations are performed for $L=13, {\omega }_0/t_0=10$ and different coupling strengths $\gamma /t_0=5,15,20$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=10,30,40$, respectively, and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every twenty-fifth data point that was computed.

Figure 9

Time evolution of the kinetic energy per fermion $E_{\mathrm{kin}}/N$ when starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). The dashed horizontal lines represent the kinetic energy in the ground state at the respective parameters. Simulations are performed for $L=13, {\omega }_0/t_0=2$ and different coupling strengths $\gamma /t_0=1,3,4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=10,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 10

Time evolution of (a) the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ and (b) the phonon number per fermion $N_{\mathrm{ph}}/N$ when starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). Simulations are performed for $L=13, {\omega }_0/t_0=2$ and coupling strengths $\gamma /t_0=4$. In the time evolution, we use different local phonon cutoffs $M_{\mathrm{ph}}=10,20,40$ to illustrate convergence in this parameter. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-7}$. For clarity, we only show every fifth data point that was computed.

Figure 11

Time evolution of (a) the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ and (b) the phonon number per fermion $N_{\mathrm{ph}}/N$ when starting from the dressed CDW state $|\mathrm{DCDW}\rangle$ Eq. (34). The small black dots in panel (a) are exact analytical results for $\gamma =0$ when starting from the $|\mathrm{BCDW}\rangle$ state in the thermodynamic limit [97]. The dashed horizontal lines in panel (b) represent the number of phonons per fermion in the ground states at the respective parameters. Simulations are performed for $L=13, {\omega }_0/t_0=2$ and different coupling strengths $\gamma /t_0=1,3,4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=20,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 12

Time evolution of the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ when starting from the dressed CDW state $|\mathrm{DCDW}\rangle$ Eq. (34). Here, the time axis is in units of the effective hopping matrix element ${\tilde{t}}_0$ Eq. (8) (system size $L=13$, phonon frequency ${\omega }_0/t_0=2$, and data for different coupling strengths $\gamma /t_0=1,3,4$). In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=20,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth [every twentieth] data point that was computed for $\gamma /t_0=3 [\gamma /t_0=4]$.

Figure 13

Time evolution of the kinetic energy per fermion $E_{\mathrm{kin}}/N$ when starting from the dressed CDW state $|\mathrm{DCDW}\rangle$ Eq. (34). The dashed horizontal lines represent the kinetic energy in the ground states at the respective parameters. Simulations are performed for $L=13, {\omega }_0/t_0=2$ and different coupling strengths $\gamma /t_0=1,3,4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=20,30,40$, respectively. The local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 14

Time evolution of (a) the decay of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ and (b) the staggered displacement ${\mathcal{O}}_{\mathrm{disp}}$ Eq. (35) in a quench from the CDW phase to the TLL phase. Circles: Quench from ${\omega }_{0,\mathrm{init}}/t_0=2$ and ${\gamma }_{\mathrm{init}}/t_0=4$ to ${\omega }_0/t_0=0.1$ and $\gamma /t_0=0.2$ (frequency quench, FQ). Diamonds: Quench from ${\gamma }_{\mathrm{init}}/t_0=4$ to $\gamma /t_0=1$ while ${\omega }_0/t_0=2$ is kept fixed (coupling quench, CQ). The small black dots in panel (a) are exact analytical results in the thermodynamic limit for $|\mathrm{BCDW}\rangle$ as initial state and no coupling to phonons [97]. The dashed horizontal lines in panel (b) represent the value of ${\mathcal{O}}_{\mathrm{disp}}$ in the respective ground states. The system size is $L=13$, local phonon cutoff $M_{\mathrm{ph}}=40$, and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only plot every second [fourth] data point that was computed in the FQ [CQ].

Figure 15

Time evolution of the kinetic energy per fermion $E_{\mathrm{kin}}/N$ in a quench from the CDW phase to the TLL phase. Circles: Quench from ${\omega }_{0,\mathrm{init}}/t_0=2$ and ${\gamma }_{\mathrm{init}}/t_0=4$ to ${\omega }_0/t_0=0.1$ and $\gamma /t_0=0.2$ (frequency quench, FQ). Diamonds: Quench from ${\gamma }_{\mathrm{init}}/t_0=4$ to $\gamma /t_0=1$ while ${\omega }_0/t_0=2$ is kept fixed (coupling quench, CQ). The dashed horizontal lines represent the kinetic energy per fermion $E_{\mathrm{kin}}/N$ in the respective ground states. The system size is $L=13$, local phonon cutoff $M_{\mathrm{ph}}=40$, and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only plot every second [fourth] data point that was computed in the FQ [CQ].

Figure 16

Time evolution of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ when starting from the bare CDW state $|\mathrm{BCDW}\rangle$ Eq. (32). Comparison between different system sizes $L=5$ (open black symbols), $L=9$ (open colored symbols) and $L=13$ (filled symbols). The phonon frequency is ${\omega }_0/t_0=2$ while in panel (a) $\gamma /t_0=1$, in panel (b) $\gamma /t_0=3$ and in panel (c) $\gamma /t_0=4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=10,30,40$, respectively, and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 17

Time evolution of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ when starting from the dressed CDW state $|\mathrm{DCDW}\rangle$ Eq. (34). Comparison between different system sizes $L=5$ (open black symbols), $L=9$ (open colored symbols) and $L=13$ (filled symbols). The phonon frequency is ${\omega }_0/t_0=2$ while in panel (a) $\gamma /t_0=1$, in panel (b) $\gamma /t_0=3$ and in panel (c) $\gamma /t_0=4$. In the time evolution, we use a local phonon cutoff $M_{\mathrm{ph}}=20,30,40$, respectively, and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only show every fifth data point that was computed.

Figure 18

Time evolution of the charge-density-wave order parameter ${\mathcal{O}}_{\mathrm{CDW}}$ for different system sizes $L=5$ (open black symbols), $L=9$ (open colored symbols) and $L=13$ (filled symbols). Panel (a): Quench from ${\omega }_{0,\mathrm{init}}/t_0=2$ and ${\gamma }_{\mathrm{init}}/t_0=4$ to ${\omega }_0/t_0=0.1$ and $\gamma /t_0=0.2$ (frequency quench). Panel (b): Quench from ${\gamma }_{\mathrm{init}}/t_0=4$ to $\gamma /t_0=1$ while ${\omega }_0/t_0=2$ is kept fixed (coupling quench). The local phonon cutoff is $M_{\mathrm{ph}}=40$ and the local discarded weight is set to ${\mathrm{\Delta }}_{\mathrm{loc}}={10}^{-8}$. For clarity, we only plot every second data point that was computed in panel (a) and every fourth data point that was computed in panel (b).