Pump-probe spectroscopy of the one-dimensional extended Hubbard model
at half filling

Recent years have seen significant progress in studying novel phenomena in materials driven out of equilibrium by high-intensity laser irradiation de la Torre et al. (2021); Dong et al. (2022). Emergent phenomena, such as photoinduced (photoenhanced) superconductivity Kaiser (2017); Cavalleri (2018); Fausti et al. (2011); Suzuki et al. (2019); Isoyama et al. (2021); Buzzi et al. (2020); Werner et al. (2018); Wang et al. (2018); Bittner et al. (2019); Schlawin et al. (2019); Werner et al. (2019) and optical control of magnetic order Kirilyuk et al. (2010); Ishihara (2019); Mikhaylovskiy et al. (2015); Baierl et al. (2016); Afanasiev et al. (2019); Schlauderer et al. (2019); Mentink et al. (2015); Li et al. (2018) in strongly correlated electron systems, have also started to attract both experimental and theoretical attention. In this context, pump-probe spectroscopy has been a crucial tool for exploring nonequilibrium phenomena in these materials. This technique involves exciting a system with a high-intensity pump pulse, followed by examining the dynamical properties of the nonequilibrium state through the linear response of a low-intensity probe pulse.

Analyzing the nonequilibrium state of strongly correlated systems within a one-dimensional (1D) model provides a beneficial starting point since the theoretical treatment of a 1D system is simpler compared to two- or three-dimensional systems. Despite being a special case, essential characteristics of correlated systems can be elucidated from 1D systems. Furthermore, the 1D Mott insulator, which manifests in actual materials and has been thoroughly studied in the literature, is of particular interest. Ultrafast phenomena of 1D Mott insulators, such as organic salt ET-F₂TCNQ Okamoto et al. (2007); Uemura et al. (2008); Wall et al. (2011); Miyamoto et al. (2019); Takamura et al. (2023) and halogen-bridged transition-metal compounds Iwai et al. (2003); Matsuzaki et al. (2006, 2014), have been intensively explored in preceding experiments because their optical response is significantly influenced by the formation of doublon-holon bound states due to nonlocal interactions resulting from photoexcitation.

From a theoretical perspective, photoinduced nonequilibrium states in the model with nonlocal interaction have been discussed in the literature Gomi et al. (2005); Al-Hassanieh et al. (2008); Takahashi et al. (2008); Lu et al. (2012); Bittner et al. (2019); Shao et al. (2019); Murakami et al. (2022, 2023). Numerical simulations of pump-probe spectroscopy have reported the emergence of light-induced in-gap states in the half-filled 1D extended Hubbard model (1DEHM) including nearest-neighbor Coulomb interaction. This was achieved by examining the nonequilibrium optical conductivity using time-dependent exact diagonalization Lu et al. (2015) and density-matrix renormalization group Rincón and Feiguin (2021) techniques. However, previous studies have shown system-size dependencies in their results since their calculations were conducted using finite clusters.

In this paper, by employing time-dependent tensor-network algorithms in the infinite matrix-product state (iMPS) representation, we simulate the nonequilibrium dynamics of the 1DEHM directly in the thermodynamic limit to investigate nonequilibrium phenomena induced by an intense, short-time pulse. We elucidate the dynamical properties of this system by examining the linear response of a subsequent weak probe pulse, corresponding to pump-probe spectroscopy experiments. Our focus lies primarily on optical conductivity and time- and angle-resolved photoemission spectroscopy (TARPES) in nonequilibrium situations.

Calculating with an infinite system enables distinct differentiation between excitation spectra arising from continuous and discrete levels. In finite-system analyses, spectra may appear discretized even if the excitation derives from continuous levels. By directly simulating an infinite system, this issue can be circumvented. Additionally, thanks to the high-resolution calculation concerning momentum, we can easily obtain detailed peak structure and dispersion relations, facilitating unambiguous comparisons between theory and experiments. Furthermore, we demonstrate that the nonequilibrium dynamics can be comprehensively understood by studying both optical conductivity and TARPES in a complementary fashion. Our study aims to deepen the understanding of nonequilibrium phenomena in strongly correlated electron systems and provide valuable insights for experimental investigations.

The remainder of this paper is organized as follows. In Sec. II, we introduce the 1DEHM and describe the numerical method used in this paper. In Sec. III, we present the numerical results of nonequilibrium optical conductivity induced by an intense pump pulse. In Sec. IV, we demonstrate the simulated single-particle excitation spectra, which are expected to be observed in TARPES experiments. Finally, we provide conclusions and future outlook in Sec. V.

In this section, we introduce the Hamiltonian of the 1DEHM and provide a brief explanation of the numerical calculations that we performed.

We estimate the optical conductivity by examining the response of an electric current to a weak electric field. The current operator in a vector potential $A(t)$ is written as

$\displaystyle\hat{J}_{A(t)}=-\partialderivative{\hat{H}}{A}=t_{\mathrm{h}}\sum_{j,\sigma}\quantity(ie^{iA(t)}\hat{c}^{\dagger}_{j,\sigma}\hat{c}_{j+1,\sigma}+\mathrm{H.c.}).$

(3)

Upon applying a weak electric field $E_{\mathrm{pr}}(t)=-\partialderivative{A_{\mathrm{pr}}(t)}{t}$ as a probe pulse in addition to the pump pulse, the induced deviation in the current per site satisfies

$$j_{\mathrm{pr}}(t)=\frac{1}{L}\quantity(\langle\hat{J}_{A+A_{\mathrm{pr}}}\rangle_{t}-\langle\hat{J}_{A}\rangle_{t})=\int^{t}_{-\infty}\sigma(t,t^{\prime})E_{\mathrm{pr}}(t^{\prime})\differential{t^{\prime}},$$

(4)

where $L$ is the system size and $\sigma(t,t^{\prime})$ represents the linear-response function for the electric field. Since the response function should satisfy causality, i.e., $\sigma(t,t^{\prime})=\theta(t-t^{\prime})\sigma(t,t^{\prime})$, Eq. (4) can be expressed as $j_{\mathrm{pr}}(t)=\int^{\infty}_{-\infty}\sigma(t,t^{\prime})E_{\mathrm{pr}}(t^{\prime})\differential{t^{\prime}}$. Taking the Fourier transform of both sides with respect to $t$, we obtain $j_{\mathrm{pr}}(\omega)=\int^{\infty}_{-\infty}\sigma(\omega,t^{\prime})E_{\mathrm{pr}}(t^{\prime})e^{i\omega t^{\prime}}\differential{t^{\prime}}$, where $\sigma(\omega,t^{\prime})\equiv\int^{\infty}_{-\infty}\sigma(t,t^{\prime})e^{i\omega(t-t^{\prime})}\differential{t}$. Assuming that the probe pulse $E_{\mathrm{pr}}(t)$ is nonzero only within the period $t_{\mathrm{pr}}\pm\tau/2$, where $\tau$ is much smaller than the characteristic time scale of the system, and that $\sigma(\omega,t^{\prime})$ remains constant during this period, we approximate $j_{\mathrm{pr}}(\omega)\simeq\sigma(\omega,t_{\mathrm{pr}})\int^{\infty}_{-\infty}E_{\mathrm{pr}}(t^{\prime})e^{i\omega t^{\prime}}\differential{t^{\prime}}=\sigma(\omega,t_{\mathrm{pr}})E_{\mathrm{pr}}(\omega)$ Shao et al. (2016). In this case, the optical conductivity of the frequency $\omega$ at probe time $t_{\mathrm{pr}}$ can be evaluated from

$\displaystyle\sigma(\omega,t_{\mathrm{pr}})=\frac{j_{\mathrm{pr}}(\omega)}{i(\omega+i\eta)A_{\mathrm{pr}}(\omega)},$

(5)

where $A_{\mathrm{pr}}(\omega)=E_{\mathrm{pr}}(\omega)/i(\omega+i\eta)$ with a damping factor $\eta$. Even though this factor is introduced for the convergence of our numerical Fourier transformation, it is associated with the lifetime of the quasiparticles due to, for example, impurity scattering in actual materials. In this way, we can directly determine equilibrium and nonequilibrium optical conductivities in the thermodynamic limit by means of the iTEBD method. We adopt a weak, narrow probe pulse $A_{\mathrm{pr}}(t)=A_{0,\mathrm{pr}}e^{-(t-t_{\mathrm{pr}})^{2}/2\sigma_{\rm pr}^{2}}$, where $\omega A_{0,\mathrm{pr}}\sigma_{\mathrm{pr}}\ll 1$ is satisfied.

The advantage of this method is that it allows simultaneous calculation of the response to an external field across all frequencies. This is because the probe pulse can be regarded as a delta function when $\sigma_{\mathrm{pr}}$ is very small compared to the characteristic time scale of the system. In other words, this probe pulse is a superposition of waves at all frequencies. In particular, when assuming $\sigma_{\mathrm{pr}}\to 0$, Eq. (5) yields the same result as the optical conductivity obtained by applying the Kubo formula, originally formulated for thermal systems, to a nonequilibrium state Shao et al. (2016); Rincón and Feiguin (2021). Note that the same method was used to simulate the nonequilibrium optical conductivity in previous studies Lu et al. (2015); Shao et al. (2016); Shinjo and Tohyama (2018); Rincón and Feiguin (2021); Shinjo et al. (2022). The ultrashort probe pulses used in this paper are idealized, and in actual pump-probe spectroscopy experiments, the probe pulse is a wave packet with a finite width. The interpretation of nonequilibrium optical conductivity is complicated by the uncertainty relation between energy and time. There is an ongoing debate about the theoretical description of the optical conductivity observed in actual pump-probe spectroscopy Shao et al. (2016); Eckstein and Kollar (2008); Kennes et al. (2017). A detailed quantitative analysis to reconcile the experimental results remains a subject for future work. In the following, we set $A_{0,\mathrm{pr}}=0.05$, $\sigma_{\rm pr}=0.05$, and $\eta=0.1$ for the computations of optical conductivity. We denote $\Delta t_{\mathrm{pr}}=t_{\mathrm{pr}}-t_{0}$ and rewrite Eq. (5) as $\sigma(\omega,\Delta t_{\mathrm{pr}})$.

Figure 1(a) shows the real parts of the optical conductivity of the SDW states for $U=10$ and various $V$ in the absence of the pump pulse [$A(t)=0$], which is consistent with previous dynamical density-matrix renormalization group studies Jeckelmann (2002, 2003). A broad peak appears above a charge gap at $V=0$, originating from the excitations between the continuous levels of the upper-Hubbard band (UHB) and the lower-Hubbard band (LHB). Turning on the intersite interaction $V$, the energy level of a doublon-holon bound state (exciton) emerges, and decreases as $V$ increases. The energy level of the exciton becomes smaller than the bottom of the energy continuum for $V\geq 2$, leading to the emergence of a sharp peak below the charge gap Gallagher and Mazumdar (1997). We should note that in our numerical calculations the amplitude of this excitonic peak is finite due to the finite $\eta$ and diverges as $\eta\to 0$ Jeckelmann (2003). The excitonic energy level becomes the lowest value at the SDW-CDW transition point $V\simeq U/2$.

The optical conductivities of the CDW state $V>U/2$ are shown in Fig. 1(b). Here, the peak position of $\real\sigma(\omega)$ increases as $V$ increases. This peak position corresponds to the energy required to dissociate a doublon.

The peak positions of the optical conductivities can be readily estimated in the strong-coupling limit ($U,V\gg t_{\mathrm{h}}$). In this limit, the SDW state comprises singly occupied sites. Thus, the energy of the first-excited state, characterized by the presence of an adjacent doublon and holon, is approximately $U-V$. On the other hand, doublons and holons align alternately in the CDW state. Therefore, the first-excited state, where two adjacent sites become singly occupied, requires an excitation energy of approximately $3V-U$.

We now discuss the time-dependent single-particle excitation spectra in the 1DEHM, taking into account TARPES experiments. The intensity of single-particle excitation spectra with momentum $k$ and energy $\omega$ is given by Freericks et al. (2009, 2015)

	$\displaystyle A^{-}(k,\omega,\Delta t_{\mathrm{pr}})$
	$\displaystyle=\frac{1}{L}\sum_{j,\ell,\sigma}e^{-ik(r_{j}-r_{\ell})}\int^{\infty}_{-\infty}\differential{t_{1}}\int^{\infty}_{-\infty}\differential{t_{2}}e^{i\omega(t_{1}-t_{2})}$
	$\displaystyle\quad\times s(t_{1}-t_{\mathrm{pr}})s(t_{2}-t_{\mathrm{pr}})\expectationvalue{\hat{c}^{\dagger}_{\ell,\sigma}(t_{2},-\infty)\hat{c}_{j,\sigma}(t_{1},-\infty)}_{-\infty},$		(6)

where $s(t-t_{\mathrm{pr}})$ is an envelope function of the wave packet of the probe pulse at central time $t_{\mathrm{pr}}$ and $\hat{c}_{j,\sigma}(t,t^{\prime})=\hat{U}^{\dagger}(t,t^{\prime})\hat{c}_{j,\sigma}\hat{U}(t,t^{\prime})$ is the Heisenberg representation. The envelope function utilized in this paper is a Gaussian function written as

$$s(t-t_{\mathrm{pr}})=\frac{1}{\sqrt{2\pi}\sigma_{\mathrm{pr}}}\exp(-\frac{\quantity(t-t_{\mathrm{pr}})^{2}}{2\sigma_{\mathrm{pr}}^{2}}),$$

(7)

where $\sigma_{\mathrm{pr}}$ represents the width of the wave packet. To simulate single-particle excitation spectra in nonequilibrium, we construct the window state from the iMPS obtained from the iTEBD method and employ the infinite-boundary conditions with a uniform update scheme Zauner et al. (2015). For more detail, refer to Ref. Ejima et al. (2022). In addition, we define the integrated photoemission spectra as

$$A^{-}(\omega,\Delta t_{\mathrm{pr}})=\int^{\pi}_{-\pi}\frac{\differential{k}}{2\pi}A^{-}(k,\omega,\Delta t_{\mathrm{pr}}),$$

(8)

which in the following will be denoted as the time-resolved density of states (TDOS).

Increasing the width of the probe-pulse wave packet improves the energy resolution, but decreases the time resolution due to the uncertainty relation. In this section, we set the width of the probe pulse to $\sigma_{\mathrm{pr}}=3$. We find that a window-state size of $L_{\mathrm{w}}=64$ is sufficient for this case. We present calculations for the single-particle excitation spectra of nonequilibrium states at $\Delta t_{\mathrm{pr}}=0$ and $8$. However, we have confirmed that the spectral shape for $\Delta t_{\mathrm{pr}}>8$ is almost unchanged from that for $\Delta t_{\mathrm{pr}}=8$, implying that the single-particle excitation spectra are essentially stationary over time after photoexcitation.

We explored the pump-probe spectroscopy of the 1DEHM at half filling in an infinite system. In the strong-coupling regime, the GS of this model resides in the SDW phase when $U\gtrsim 2V$ and in the CDW phase when $U\lesssim 2V$. In the SDW phase, doublons and holons, which are generated by the intense pump pulse, form bound states known as excitons through nonlocal interactions. In the CDW phase, the charge order dissolves due to photoexcitation. The dynamical response in the nonequilibrium state was revealed using the iTEBD method.

We detected an in-gap state at small $\omega$ in the nonequilibrium optical conductivity for the model with $U=10$ and $V=3$, which resides in the SDW phase. This state can be interpreted as the transition between the odd-parity exciton and the even-parity exciton. Furthermore, we investigated a stronger interaction model with $U=20$ and $V=6$. In addition to this in-gap state originating from the different parity excitons, we discovered that additional peaks appear below and above the excitonic energy. The origin of these additional peaks can be understood by examining the single-particle excitation spectra. Specifically, the new dispersions appearing above and below the LHB after the pump pulse irradiation are associated with these new peaks in the optical conductivity. We also observed that the LHB is replicated in a higher-energy region, where the energy difference is equal to the excitonic energy.

Moreover, we examined the nonequilibrium optical conductivity for the model with $U=10$ and $V=6$, which resides in the CDW phase. In this case, we also found an in-gap state after the pump pulse irradiation. Additionally, we discovered that the Drude weight becomes finite, suggesting the metallization of the system. The origin of this in-gap state can be understood from the single-particle excitation spectra.

It would also be interesting to explore nonequilibrium physics close to SDW-CDW phase boundaries, although this paper focused on pump-probe spectroscopy deep in the SDW and CDW phases. For instance, the time-resolved single-particle spectral function has been studied around these phase boundaries by means of the exact-diagonalization technique Shao et al. (2020). Furthermore, a recent study indicates peculiar optical responses in high-harmonic generation near the phase boundary Shao et al. (2022). The issue with performing iTEBD simulations near the quantum phase transition point is the increasing bond dimensions required. It is thus highly desirable to improve the accuracy of numerical techniques and simultaneously reduce the computational time.

Lastly, we address the correspondence between our theoretical findings and experimental observations. The 1D Mott insulator ET-F₂TCNQ is well described by the 1DEHM with interaction parameters $U=10$ and $V=3$ Yamaguchi et al. (2021). In fact, the emergence of the in-gap state in optical conductivity after pump pulse irradiation has been observed. If TARPES becomes feasible in this material, we expect that the single-particle excitation spectra obtained in this paper would also be observable. Another approach to realizing our results is by employing a cold atomic system Bohrdt et al. (2018). With an artificial gauge field mimicking the pump pulse, the observation of single-particle excitation spectra may be feasible. It is worth emphasizing that our calculations, performed on an infinite system, allow for a direct comparison of the spectra observed in future experiments with our results.

In this paper, our focus was on nonequilibrium optical conductivity and single-particle excitation spectra. Recently, time-resolved resonant inelastic x-ray scattering (RIXS) spectra have become observable through pump-probe spectroscopy Cao et al. (2019); Mitrano and Wang (2020); Gel’mukhanov et al. (2021). RIXS allows us to examine the dynamical correlations of charge and spin, including their momentum dependence, thereby enabling us to obtain more detailed information on strongly correlated materials. We anticipate further developments in theoretical studies of pump-probe spectroscopy using tensor-network algorithms in the future.

K.S. was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grants No. JP20H01849, No. JP21K03439, and No. JP23K03286 and by Japan Science and Technology Agency (JST) COI-NEXT Program Grant No. JPMJPF2221. S.E. was supported by the DLR Quantum Computing Initiative and the Federal Ministry for Economic Affairs and Climate Action dlr . The iTEBD calculations were performed using the ITensor library Fishman et al. (2022).