# Optical response of laser-driven charge-transfer complex described by
Holstein–Hubbard model coupled to heat baths: Hierarchical equations of motion
approach

Kiyoto Nakamura<EMAIL_ADDRESS>Yoshitaka Tanimura
<EMAIL_ADDRESS>Department of Chemistry, Graduate School of
Science, Kyoto University, Sakyoku, Kyoto 606-8502, Japan

###### Abstract

We investigate the optical response of a charge-transfer complex in a
condensed phase driven by an external laser field. Our model includes an
instantaneous short-range Coulomb interaction and a local optical vibrational
mode described by the Holstein–Hubbard (HH) model. Although characterization
of the HH model for a bulk system has typically been conducted using a complex
phase diagram, this approach is not sufficient for investigations of dynamical
behavior at finite temperature, in particular for studies of nonlinear optical
properties, where the time irreversibility of the dynamics that arises from
the environment becomes significant. We therefore include heat baths with
infinite heat capacity in the model to introduce thermal effects characterized
by fluctuation and dissipation to the system dynamics. By reducing the number
of degrees of freedom of the heat baths, we derive numerically “exact”
hierarchical equations of motion (HEOM) for the reduced density matrix of the
HH system. As demonstrations, we calculate the optical response of the system
in two- and four-site cases under external electric fields.The results
indicate that the effective strength of the system–bath coupling becomes large
as the number of sites increases. Excitation of electrons promotes the
conductivity when the Coulomb repulsion is equivalent to or dominates the
electron–phonon coupling, whereas excitation of optical vibrations always
suppresses the conductivity.

##  I INTRODUCTION

In various organic conductors such as bis(ethylenedithio)tetrathiafulvalene
(BEDT-TTF)Saito and Yoshida (2007) and the fullerides (salts of C60
anions),Gunnarsson (1997) the interplay of on-site electron–electron and on-
site phonon–electron interactions plays a significant role. Owing to the
competition between these two interactions, such materials exhibit a variety
of distinct states in the antiferromagnetic (AF), charge-density-wave (CDW),
spin-density-wave (SDW), and superconducting (SC) phases. For a description of
the complex phase diagram of correlated electron systems, the Holstein–Hubbard
(HH) model has been introduced.Holstein (1959a, b); Beni, Pincus, and Kanamori
(1974); Aubry (1995); Takada and Chatterjee (2003); Sykora, Hübsch, and Becker
(2008); Miyao (2018); Berciu (2007); Yamada, Ishizuka, and Ōno (2014) This
model consists of electrons on a discrete lattice, with on-site Coulomb
repulsion between electrons and a local coupling of electrons to longitudinal
optical (LO) phonons. It is characterized by the transfer integral $t_{h}$,
the on-site Coulomb repulsion $U$, the LO phonon frequency $\Omega$, and the
electron–phonon coupling strength $g$.

For the half-filled HH model, when the Coulomb repulsion $U$ is much larger
than the electron–phonon coupling $g$, two electrons with opposite spin
directions cannot occupy the same site, and electrons with up spin and those
with down spin are aligned in an alternating sequence in the Mott AF phase,
whereas when $g$ is much larger than $U$, two electrons with opposite spin
directions occupy the same site with the aid of an LO phonon, and bipolarons
are formed in the Peierls CDW phase. Then, between the AF and CDW phases, the
existence of a SC state is predicted, which is responsible for the high-
temperature SC state of alkali fullerides.Hebard  _et al._ (1991); Rosseinsky
_et al._ (1991); Palstra  _et al._ (1995) The role of phonons in the SC state
of cuprates has also been investigated along the same lines,Lanzara  _et al._
(2001); Huang  _et al._ (2003); Johnston  _et al._ (2012); Ishihara and
Nagaosa (2004); Mishchenko  _et al._ (2008); De Filippis  _et al._ (2009) and
the HH model, despite its simplicity, does seem to represent a versatile
approach for the investigation of these materials.

Many numerical studies, based, for example, on exact diagonalization, Takada
(1996); Chatterjee (2010) the density matrix renormalization group (DMRG),
Tezuka, Arita, and Aoki (2007) and the quantum Monte Carlo (QMC) method, Von
der Linden, Berger, and Valášek (1995); Macridin, Sawatzky, and Jarrell
(2004); Clay and Hardikar (2005); Nowadnick  _et al._ (2012); Johnston  _et
al._ (2013); Ohgoe and Imada (2017); Weber and Hohenadler (2018) have been
conducted to reveal the thermodynamic properties of the HH model at zero and
finite temperatures. Most of these investigations have focused on bulk
systems, typically through the use of a complex phase diagram, and, with a few
exceptions (see, e.g., Ref. Tsuji, Oka, and Aoki, 2009), there has been little
exploration of the dynamical behavior of the HH model. The recent discovery of
SC-like properties initiated by optical pulses and observed at temperatures
higher than the SC critical temperature $T_{c}$ provides considerable
motivation for dynamical investigations of the HH model,Mitrano  _et al._
(2016) as do the current efforts to achieve room-temperature superconductivity
with the aid of optical pulses, in the context of which both theoretical Mazza
and Georges (2017); Sentef  _et al._ (2016); Kennes  _et al._ (2017); Kaneko
_et al._ (2019); Ido, Ohgoe, and Imada (2017); Kim  _et al._ (2016); Werner
and Eckstein (2015); Golež  _et al._ (2012) and experimentalCantaluppi  _et
al._ (2018); Nicoletti  _et al._ (2017); Mankowsky  _et al._ (2017)
investigations have already been conducted.

To investigate dynamical behavior of an HH system, not only the LO phonon
modes but also environmental degrees of freedom have to be taken into account,
because the latter control the system temperature and the time scales of
relaxation processes of excited states through fluctuations and dissipation
related through the fluctuation–dissipation theorem. However, investigations
of dynamics, in particular investigations of nonlinear optical properties at
finite temperature, have not been conducted thoroughly, even for a small
cluster system.

In this paper, we develop an open quantum dynamics theory for the HH model, in
which the HH system is coupled to baths consisting of surrounding atoms or
molecules. Here, by an open quantum dynamics theory we mean a theory
pertaining to a system that is further coupled to bath systems that are
modeled by infinite numbers of harmonic oscillators.Weiss (2012) Such
system–bath models have been employed for investigations of electron
transferGarg, Onuchic, and Ambegaokar (1985) and exciton transfer
processesIshizaki and Fleming (2009); Hein  _et al._ (2012) of a molecular
system, in particular for calculations of the nonlinear optical spectrum.
However, the effects of on-site Coulomb repulsion between electrons, which
plays an essential role in charge transfer complexes, has not been
investigated in any great depth.

We derive hierarchical equations of motion (HEOM) for the reduced density
matrix of the HH system that can treat non-Markovian and nonperturbative
system–bath interactions at finite temperature in a numerically “exact”
manner. Tanimura and Kubo (1989); Tanimura (2006, 2020) The HEOM approach has
been applied to the HolsteinChen, Zhao, and Tanimura (2015); Zhou  _et al._
(2016); Dunn, Tempelaar, and Reichman (2019) and Peierls–HolsteinCainelli and
Tanimura (2021) problems. Photosynthetic systems have also been investigated
using this approach, and it has been shown that the efficiency of the exciton
transfer process exhibits a maximum in the nonperturbative region of
system–bath coupling,Ishizaki and Fleming (2009) owing to the quantum
entanglement between the system and bath (“bathentanglement”).Tanimura (2020)
Using the HEOM approach, we can explore the possibility of finding such a
mechanism in this HH system. We illustrate our formalism by calculating the
optical conductivity of the one-dimensional HH model at finite temperature. In
addition, by applying time-dependent external fields, we calculate the optical
conductivity under highly nonequilibrium conditions. Because the heat bath
facilitates thermal dissipation and excitation in a thermodynamically
consistent manner, we can obtain reliable results under such conditions, even
when the molecular system described by the HH model is very small.

The remainder of the paper is organized as follows. In Sec. II, we introduce
the system–bath Hamiltonian for the HH model. In Sec. III, we present the HEOM
for the HH model. Section IV is devoted to a presentation of the numerical
results for two- and four-site HH-plus-bath models. The optical conductivities
of the equilibrium states and the steady states under external driving fields
are investigated. We conclude the paper in Sec. V.

##  II MODEL

We consider a one-dimensional HH model described by

$\displaystyle\hat{H}_{S}(t)=$
$\displaystyle-\sum_{i,\sigma}(t^{i}_{h}e^{iA(t)}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}+\mathrm{H.c.})+\sum_{i}U^{i}\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}$
$\displaystyle+\sum_{i}g^{i}(\hat{n}_{i,\uparrow}+\hat{n}_{i,\downarrow})(\hat{b}_{i}^{\dagger}+\hat{b}_{i})+\sum_{i}\Omega^{i}\hat{b}_{i}^{\dagger}\hat{b}_{i}$
$\displaystyle-E_{\mu}(t)\sum_{i}(\hat{b}_{i}^{\dagger}+\hat{b}_{i}),$ (1)

where $c_{i,\sigma}$ and $c_{i,\sigma}^{\dagger}$ are the annihilation and
creation operators of an electron at site $i$ with spin state
$\sigma\in\\{\uparrow,\downarrow\\}$, $b_{i}$ and $b_{i}^{\dagger}$ are the
annihilation and creation operators of an LO phonon at site $i$, and
$t^{i}_{h},U^{i},g^{i}$, and $\Omega^{i}$ are the transfer integral, on-site
Coulomb repulsion, electron–phonon coupling strength, and LO phonon frequency,
respectively, at site $i$. The functions $A(t)$ and $E_{\mu}(t)$ are the time-
dependent external perturbations for electrons and LO phonons, described by
the vector potential and electric field; the static case is given by $A(t)=0$
and $E_{\mu}(t)=0$. Because we treat each site explicitly, these parameter
values can be site-dependent. For example, using a machine learning approach,
we can set a system-specific parameter value targeting a specific complex
molecular system.Ueno and Tanimura (2021) The probability distribution of the
electron is expressed as
$\hat{n}_{i,\sigma}=\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$. Here and
in the rest of this paper, we set the reduced Planck constant and the
elementary charge equal to $1$, i.e., $\hbar=e=1$.

In this study, we consider a case in which the LO modes are coupled to
intermolecular vibrations or other molecular degrees of freedom, and we
describe these degrees of freedom using a harmonic heat bath. Because we can
study only a small system with a finite number of phonon modes associated with
the system site, inclusion of the heat bath is important to maintain the
stability of the equations.Dunn, Tempelaar, and Reichman (2019) Thus, we
consider the situation in which each LO oscillator interacts with a heat bath
that gives rise to dissipation and fluctuation in the LO modes. The total
Hamiltonian is then expressed as

$\displaystyle\hat{H}_{\mathrm{tot}}(t)=\hat{H}_{S}(t)+\hat{H}_{I}+\hat{H}_{B}.$
(2)

The bath Hamiltonian $\hat{H}_{B}$ and the system–bath interaction
$\hat{H}_{I}$ are expressed as

$\displaystyle\hat{H}_{B}$
$\displaystyle=\sum_{i,\alpha}\left(\frac{\hat{p}_{i,\alpha}^{2}}{2m_{i,\alpha}}+\frac{1}{2}m_{i,\alpha}\omega_{i,\alpha}^{2}\hat{x}_{i,\alpha}^{2}\right),$
(3) $\displaystyle\hat{H}_{I}$
$\displaystyle=-\sum_{i}\hat{V}_{i}\sum_{\alpha}k_{i,\alpha}\hat{x}_{i,\alpha}+\sum_{i}\hat{V}_{i}^{2}\sum_{\alpha}\frac{k_{i,\alpha}^{2}}{2m_{i,\alpha}\omega_{i,\alpha}^{2}},$
(4)

where $\hat{p}_{i,\alpha},\hat{x}_{i,\alpha},m_{i,\alpha}$, and
$\omega_{i,\alpha}$ are the momentum, position, mass, and frequency of the
$\alpha$th oscillator of the $i$th bath, respectively. The system operator
that represents the effects of the $i$th local heat bath is given by
$\hat{V}_{i}$, and $k_{i,\alpha}$ is the coupling constant for the interaction
between the system and the $\alpha$th oscillator of the $i$th bath. The second
term in $\hat{H}_{I}$ is a counterterm that compensates for the
renormalization of the potential energy caused by the first term in
$\hat{H}_{I}.$Weiss (2012) Depending on the situation, we may choose a local
or nonlocal system–bath interaction by the way in which we set $\hat{V}_{i}$:
for high-frequency intramolecular vibrational modes, we consider the case in
which each LO oscillator is coupled to its own bath, which is expressed as
$\hat{V}_{i}=\hat{b}_{i}^{\dagger}+\hat{b}_{i}$, whereas for low-frequency
intermolecular vibrational modes or phonon modes, we consider the case in
which multiple LO oscillators are coupled to a single global bath, described
as $\hat{V}_{i}\to\hat{V}_{0}=\hat{b}^{\dagger}_{0}+\hat{b}_{0}$, with $i$
being ignored in the summations in Eqs. (3) and (4). If necessary, we can
further treat a nonlocal heat bath described by $\hat{V}_{ij}$ between the
$i$th and $j$th sites, which often plays an important role in determining the
efficiency of an energy transfer process.Cainelli and Tanimura (2021); Ueno
and Tanimura (2021)

The $i$th heat bath can be characterized by the spectral distribution function
(SDF), defined by

$\displaystyle
J_{i}(\omega)=\sum_{\alpha}^{N_{i}}\frac{k^{2}_{i,\alpha}}{2m_{i,\alpha}\omega_{i,\alpha}}\delta(\omega-\omega_{i,\alpha}).$
(5)

For the heat bath to be an unlimited heat source possessing an infinite heat
capacity, the number of heat-bath oscillators $N_{i}$ is effectively made
infinite by replacing $J_{i}(\omega)$ with a continuous distribution. The
counterterm of $\hat{H}_{I}$ can then be rewritten in terms of the
reorganization energy
$\lambda_{i}=\int_{0}^{\infty}d\omega\,J_{i}(\omega)/\omega$ as
$\sum_{i}\lambda_{i}\hat{V}_{i}^{2}$.

Although we can reduce the number of phonon degrees of freedom by
diagonalizing the LO phonon plus bath degrees of freedom, which leads to the
Hubbard model with Brownian-type heat bath, as described in Refs. Garg,
Onuchic, and Ambegaokar, 1985; Ueno and Tanimura, 2021; Tanimura, 1990; Tanaka
and Tanimura, 2009, 2010; Tanimura, 2012; Cainelli and Tanimura, 2021, here we
treat the LO phonon modes expressed in the eigenstate representation
explicitly to reduce the computational cost. In this way, we can investigate
the case in which the LO phonon modes are directly excited by a laser field
(see Sec. IV.2.3). Moreover, if necessary, we can include an anharmonicity of
the phonon modes by setting each excitation energy of the phonon state
separately without incurring any additional cost.

##  III HEOM FOR THE HH-PLUS-BATH SYSTEM

Our open quantum dynamics theory is constructed on the basis of the reduced
density operator, defined as
$\hat{\rho}(t)=\operatorname{tr}_{B}\\{e^{-i\hat{H}_{\mathrm{tot}}(t)}\hat{\rho}_{\mathrm{tot}}(0)e^{i\hat{H}_{\mathrm{tot}}(t)}\\}$,
where $\operatorname{tr}_{B}\\{~{}\\}$ represents the partial trace over the
bath degrees of freedom.

If we adopt the Drude–Lorentz SDF expressed as

$\displaystyle
J_{i}(\omega)=\frac{\eta_{i}\omega}{\pi}\frac{1}{1+(\omega/\gamma_{i})^{2}},$
(6)

then the time evolution of the reduced density operator is described by the
HEOMTanimura (2020) for the HH-plus-bath model as follows (see Appendix A):

$\displaystyle\frac{\partial}{\partial
t}\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{N}}(t)={}$
$\displaystyle-\left(i\hat{H}_{S}^{\times}(t)+\sum_{i=1}^{N}\sum_{k=0}^{K}n_{i,k}\nu_{i,k}\right)\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{N}}(t)$
$\displaystyle-\sum_{i=1}^{N}\left(i\lambda_{i}(\hat{V}_{i}^{2})^{\times}\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{N}}(t)+i\hat{V}_{i}^{\times}\sum_{k=0}^{K}\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{i}+\vec{e}_{k},\ldots,\vec{n}_{N}}(t)\right.$
$\displaystyle\hskip
32.0pt+i\hat{V}_{i}^{\times}\sum_{k=0}^{K}n_{i,k}c_{i,k}\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{i}-\vec{e}_{k},\ldots,\vec{n}_{N}}(t)\left.+\hat{V}_{i}^{\circ}\frac{\eta_{i}\nu_{i,0}^{2}}{2}\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{i}-\vec{e}_{0},\ldots,\vec{n}_{N}}(t)\right).$
(7)

Here, the vector $\vec{n}_{i}=(n_{i,0},n_{i,1},\ldots,n_{i,K})$ consists of
nonnegative integers, and $\hat{\rho}_{\vec{0},\ldots,\vec{0}}(t)$ corresponds
to the reduced density operator. The density operators whose vectors
$\\{\vec{n}_{i}\\}$ involve positive integers are referred to as the auxiliary
density operators (ADOs). The sets $\\{c_{i,k}\\}$ and $\\{\nu_{i,k}\\}$ are
evaluated for the description of the real part of the two-time correlation
function of the heat bath, $C^{\prime}_{i}(t)$, and the quantity $K$ is the
number of expansion coefficients for $C^{\prime}_{i}(t)$. The number of sites
is represented by $N$. The symbols $\times$ and $\circ$ represent the
commutator and anticommutator defined as
$\hat{O}^{\times}\bullet=[\hat{O},\bullet]$ and
$\hat{O}^{\circ}\bullet=\\{\hat{O},\bullet\\}$, respectively. We set
$\nu_{i,0}=\gamma_{i,0}$ to simplify the notation. The vector $\vec{e}_{k}$ is
the unit vector of the $k$th element. To conduct numerical integrations, we
set $\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{N}}(t)=0$, where
$\sum_{i,k}n_{i,k}>N_{\max}$, to truncate the equations. The number $N_{\max}$
is set sufficiently large for the calculations to converge, and it may change
depending on the parameters of the system and bath. For an outline of the
derivation of the HEOM, see Appendix A.

We truncate the eigenstates of the LO phonon in the same manner as the ADOs.
That is, the eigenstates that satisfy $\sum_{i\sigma}\hat{n}_{\sigma}>M_{LO}$
are ignored, where $M_{LO}$ is an integer representing the LO cutoff number.
With a sufficiently large value of $M_{LO}$, we obtain reliable results
because the populations in the higher eigenstates are almost zero, owing to
the rapid removal of the excess energy of the phonon modes to the heat bath.

##  IV NUMERICAL RESULTS

###  IV.1 Equilibrium distribution

We can obtain the reduced thermal equilibrium state of the total system (which
is not the Boltzmann distribution of the main system itself owing to the
contribution of the system–bath interaction) by numerically integrating Eq.
(III) from a temporal initial state to a time sufficiently long that all of
the hierarchical elements have reached a steady state.Tanimura (2014, 2015,
2006, 2020) Although we can set any strength and cutoff frequency of the
system–bath interaction, described as $\eta_{i}$ and $\gamma_{i}$, with any
form of $\hat{V}_{i}$, here we limit our analysis to the weak-coupling case to
justify the description of the present model in comparison with a case without
a heat bath.

We set the transfer integral, on-site Coulomb repulsion, electron–phonon
coupling strength, and LO phonon frequency at each site identically as
$t^{i}_{h}=t_{h}$, $U^{i}=U$, $g^{i}=g$, and $\Omega^{i}=\Omega$,
respectively, for all $i$. We choose the transfer integral $t_{h}$ as the unit
of frequency. Then, for all calculations, we fix the LO phonon frequency, the
coupling strength between electrons and LO phonons, and the cutoff frequency
of the heat bath as $\Omega/t_{h}=2$, $g/t_{h}=\sqrt{3.6}$, and
$\gamma_{i}/t{{}_{h}}=3$. We choose these parameter values because the quantum
phase transition for the SC state has been reported under such
conditions.Takada and Chatterjee (2003); Tezuka, Arita, and Aoki (2007) The
inverse temperature divided by the Boltzmann constant $k_{B}$ is denoted by
$\beta=1/k_{B}T$. We set $\beta t_{h}=1$ (high-temperature case) and $\beta
t_{h}=4$ (low-temperature case).

To examine the description of the HH-plus-bath model, we calculate the
expectation value of the system Hamiltonian, $\braket{\hat{H}_{S}}$, using the
steady solution of the HEOMZhang  _et al._ (2017) and compare the result with
the value calculated from the Boltzmann distribution of the HH system, in
which the eigenenergies are evaluated by exact diagonalization.Chatterjee
(2010) For this purpose, we consider a two-site HH model with parameter values
$U/t_{h}=3$, $M_{LO}=5$, and $\beta t_{h}=4$. For the HH-plus-bath model, we
consider a local heat bath described by the interaction
$\hat{V}_{i}=\hat{b}^{\dagger}_{i}+\hat{b}_{i}$, in the case of weak
system–bath coupling described as $\eta_{i}=0.001$ $(i=1,2)$ with the ADO
parameter values $N_{\max}=2$ and $K=2$. We obtain
$\braket{\hat{H}_{S}}=-5.00$ for the HH-plus-bath model, which agrees very
well with the value $\braket{\hat{H}_{S}}=-4.94$ obtained from diagonalization
of the HH model. It should be noted that in the case of strong system–bath
coupling, the system energy calculated from the steady-state solution of the
HEOM for the HH-plus-bath model does not agree with that calculated from the
Boltzmann distribution of the HH system, owing to the contribution of the
energy from the system–bath interaction.

###  IV.2 Optical response

We now consider the optical response of the HH-plus-bath system. While optical
signals are, in general, formulated on the basis of the optical dipole moment,
here we consider an expectation value of the optical current,
$\braket{\hat{j}(t)}=\operatorname{tr}\\{\hat{j}(t)\hat{\rho}(t)\\}$, where
$\hat{j}(t)=i\sum_{i,\sigma}(t_{h}e^{i(A(t)+A^{\prime}(t))}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}-\mathrm{H.c.})$
is the current operator, under a sufficiently weak probe excitation described
by a time-dependent vector potential $A^{\prime}(t)$. This consideration is
with future extensions to bulk material in mind. The optical conductivity is
then defined as

$\displaystyle\sigma[\omega]$
$\displaystyle=\frac{1}{E^{\prime}[\omega]}\int_{0}^{\infty}dt\,e^{i\omega
t}[\braket{\hat{j}(t)}-\braket{\hat{j}(t)}_{0}],$ (8)

where $\braket{\cdots}_{0}$ is the current without the vector potential, which
implies $A^{\prime}(t)=0$. The quantity $E^{\prime}[\omega]=i\omega
A^{\prime}[\omega]$ is the electric field that arises from the vector
potential $A^{\prime}(t)$, where we denote the Fourier transform of $f(t)$ as
$f[\omega]\equiv\int_{-\infty}^{\infty}dt\,e^{i\omega t}f(t)$. In this paper,
we calculate the optical conductivity of the HH-plus-bath system for two- and
four-site cases with the parameter values listed in Table 1.

Table 1: System and bath parameter values for two- and four-site HH-plus-bath models.  | $\beta t_{h}$ | $U/t_{h}$ | $g/t_{h}$ | $\Omega/t_{h}$ | $\gamma_{i}/t_{h}$ | $\eta_{i}$ | $N_{\max}$ | $K$ | $M_{LO}$
---|---|---|---|---|---|---|---|---|---
Two-site | $0.1$ | $1$ |  |  |  | $1$ | $30$ | $0$ | $9$
|  |  |  | $0.05$ | 3
| $1$ |  $1$, $3$, $5$ | $\sqrt{3.6}$ | $2$ | $3$ | $1$ | $2$
| $4$ |  $1$, $3$, $5$ |  |  |  | $2$
Four-site | $0.1$ | $1$ |  |  |  | $5$ | $0$ | $7$

####  IV.2.1 Equilibrium response

First, we consider linear absorption spectrum evaluated from Eq. (8) with
$A(t)=0$ and $E_{\mu}(t)=0$ in Eq. (1). In this case, the optical conductivity
is defined in terms of the linear response function as Mahan (2000); Lenarčič
_et al._ (2014)

$\displaystyle\sigma[\omega]$
$\displaystyle=\frac{i}{\omega}\\!\left(\braket{\hat{j}^{d}}_{\mathrm{eq}}-\int_{0}^{\infty}dt\,e^{i\omega
t}R^{(1)}(t)\right),$ (9)

where

$\displaystyle R^{(1)}(t)$
$\displaystyle=i\theta(t)\braket{[\hat{j}^{p}(t),\hat{j}^{p}(0)]}_{\mathrm{eq}},$
(10) $\displaystyle\hat{j}^{d}$
$\displaystyle=t_{h}\sum_{i,\sigma}(\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}+\mathrm{H.c.}),$
(11) $\displaystyle\hat{j}^{p}$
$\displaystyle=it_{h}\sum_{i,\sigma}(\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}-\mathrm{H.c.})$
(12)

are the linear response function and the diamagnetic and paramagnetic current
operators, respectively.Bergeron  _et al._ (2011) Here, $\theta(t)$ is the
step function. For details of the derivation of Eq. (9), see Appendix B.

To compare the present results with those obtained previously from the HH
model without a heat bath, we consider a fast modulation case with a weak
system–bath interaction described as $\gamma_{i}/t_{h}=3$ and $\eta_{i}=0.05$
$(i=1,2)$. We then chose the on-site Coulomb repulsion as $U/t_{h}=1$, $3$,
and $5$ and the inverse temperature as $\beta t_{h}=1$ and $4$.

Figure 1: Real part of the optical conductivity
$\operatorname{Re}\\{\sigma[\omega]\\}$ for (a) weak ($U/t_{h}=1$), (b)
intermediate ($U/t_{h}=3$), and (c) strong ($U/t_{h}=5$) Coulomb repulsion in
the low-temperature case $\beta t_{h}=4$ (blue curves) and the high-
temperature case $\beta t_{h}=1$ (red curves) calculated from Eq. (9). The LO
phonon frequency and the electron–phonon coupling strength are chosen as
$\Omega/t_{h}=2$ and $g/t_{h}=\sqrt{3.6}$, respectively. Each line shape is
normalized such that the maximum peak in the region $0\leq\omega\leq 0.3$ is
$1$. The entire peak profiles are presented in Appendix C.

It should be mentioned that, in a case of a bulk HH system without a heat
bath, we observe continuous energy transition bands. If the bulk system is a
conductor, we observe the peak at zero frequency. In the present two-site
model, however, continuous bands do not appear, owing to the finite number of
eigenstates, and the peak representing the conducting state appears at nonzero
frequency.

Figure 1 displays the low-frequency region of the real part of the optical
conductivity, $\operatorname{Re}\\{\sigma[\omega]\\}$, for the two-site HH-
plus-bath model. The peaks over the entire spectral region are presented in
Appendix C. In the low-temperature case $\beta t_{h}=4$ (blue curves), the
electronic transition peaks are observed at $0.235,0.073$, and $0.202$ for (a)
weak ($U/t_{h}=1$), (b) intermediate ($U/t_{h}=3$), and (c) strong
($U/t_{h}=5$) Coulomb repulsion, respectively. In the weak-repulsion case in
Fig. 1(a), we also observe a peak at $0.019$, which is due to the Drude weight
reflecting the sensitivity of the system to an external vector potential. This
peak appears even in the case of an insulator.Castella, Zotos, and Prelovšek
(1995); Zotos and Prelovšek (1996) The appearance of this peak depends
strongly on the system size, and hence the peak disappears as the system size
increases. Therefore, we shall not discuss this low-frequency peak further
here, and we focus our analysis only on higher transition peaks. As Figs.
1(a)–1(c) indicate, the position of the transition peak first decreases then
increases as a function of $U$. This tendency agrees with the previous results
obtained on the basis of the HH model, as presented in Ref. Tezuka, Arita, and
Aoki, 2007 (note that the present parameter value $g=\sqrt{3.6}$ is equal to
the value $\lambda=3.6$ in that paper). This implies that the weak-repulsion
case corresponds to a charge-ordered insulator, whereas the strong-repulsion
case corresponds to an AF insulator.

In the high-temperature case ($\beta t_{h}=1$), the peak position for
intermediate Coulomb repulsion in Fig. 1(b) is blue-shifted, indicating that
the insulator-like behavior is enhanced. By contrast, the peak positions in
the small- and large-$U$ cases in Figs. 1(a) and 1(c), respectively, are red-
shifted. Although the present system is not yet a conductor, the blue-shifted
peak in the intermediate-repulsion case does suggest the presence of the
conductor–insulator phase transition that is observed in a bulk HH system.

####  IV.2.2 Effects of different forms of system–bath coupling and of system
size

Above, we considered the situation in which each LO phonon mode is coupled to
its own bath. Moreover, we assumed that the system part of the system–bath
interaction was a linear function of phonon coordinate, as
$\hat{V}_{i}=\hat{b}_{i}^{\dagger}+\hat{b}_{i}\propto\hat{q}_{i}$, which leads
to population relaxation and excitation.Tanimura (2006) This form of
interaction has been used in studies of vibrational spectroscopy and electron
transfer on the basis of a Brownian oscillator model.Tanimura (2006, 2020);
Garg, Onuchic, and Ambegaokar (1985); Ueno and Tanimura (2021); Tanimura
(1990); Tanaka and Tanimura (2009, 2010); Tanimura (2012); Cainelli and
Tanimura (2021) When the vibrational excitation energy becomes close to the
thermal excitation energy, this process becomes important. If the vibrational
excitation energy is much higher than the thermal excitation energy, however,
vibrational dephasing described by $\hat{V}_{i}\propto\hat{q}_{i}^{2}$ plays a
significant role. Okumura and Tanimura (1997); Ishizaki and Tanimura (2006)
Under such conditions, this interaction can be approximated in diagonal
modulation form as
$\hat{V}_{i}\approx\hat{b}^{\dagger}_{i}\hat{b}_{i}$.Tanimura (2006, 2020) In
this subsection, we compare the different forms of interactions for
vibrational dephasing ($\hat{V}_{i}=\hat{b}^{\dagger}_{i}\hat{b}_{i}$) and
vibrational relaxation ($\hat{V}_{i}=\hat{b}^{\dagger}_{i}+\hat{b}_{i}$).
Moreover, we examine the difference between the local and nonlocal
($\hat{V}_{0}=\sum_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}$) system–bath
interactions, as well as the effects of system size.

In the following, we consider only a high-temperature case ($\beta t_{h}=0.1$)
and set $K=0$ and $c_{i,0}=\eta_{i}\gamma_{i}/\beta$ to reduce the numerical
costs of integrating the HEOM. In the present case, owing to the strong
relaxation effects, it is not easy to apply Eq. (8) to calculate the optical
conductivity. Therefore, we evaluate this from Eq. (9).

Figure 2:  Optical conductivity $\operatorname{Re}\\{\sigma(\omega]\\}$ in the
local vibrational-relaxation case,
$\hat{V}_{i}=\hat{b}^{\dagger}_{i}+\hat{b}_{i}$, in which the LO mode of each
site is modulated independently (blue curve), and that in the nonlocal bath
case, $\hat{V}_{0}=\sum_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}$, in which both of
the LO modes are coupled to a single bath (red curve). The Coulomb repulsion
and inverse temperature are set as $U/t_{h}=1$ and $\beta t_{h}=0.1$, and the
other parameter values are the same as in Fig. 1, except those for the ADOs.

In Fig. 2, we present the calculated results of
$\operatorname{Re}\\{\sigma[\omega]\\}$ for the independent-local (blue curve)
and nonlocal (red curve) bath coupling cases. The peaks in the nonlocal bath
case are narrower and not shifted in comparison with the peaks in the local
bath case, because the interaction
$\hat{V}_{0}=\sum_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}$ contributes only
indirectly to energy relaxation, and the nonlocal heat bath enhances the
intersite quantum coherence among the phonon modes.

Figure 3:  Optical conductivity $\operatorname{Re}\\{\sigma[\omega]\\}$ for
diagonal system–bath coupling with $\eta=0.05$ for the two-site ($N=2$) case
(solid black curve) and for four-site ($N=4$) cases with open (blue curve) and
periodic (red curve) boundary conditions. For the purpose of illustration, we
also depict the two-site case with stronger coupling $\eta=1$ (dashed black
curve). The Coulomb repulsion and inverse temperature are set as $U/t_{h}=1$
and $\beta t_{h}=0.1$. The other parameter values are the same as in Fig. 1,
except those for the ADOs. Note that owing to numerical error (marked in the
gray area), the four-site result with periodic boundary condition does not
approach zero at zero frequency.

Next, we discuss the effects of system size. To reduce the computational cost,
we consider the nonlocal bath case only. In Fig. 3, the blue and red curves
represent the results with $\eta=0.05$ in four-site cases with open and
periodic boundary conditions, respectively. We set $N_{\max}=5$ and
$M_{LO}=7$, with the other parameter values being the same as in Fig. 2. Note
that we set slightly smaller values of $M_{LO}$ to reduce the computational
cost for the investigation of size effects; we have verified that the line
shapes for $M_{LO}=9$ and $7$ are qualitatively the same in the two-site case
(results are not shown). For reference purposes, we also re-plot here as the
solid black curve the two-site result depicted in Fig. 2. As these results
indicate, the number of peaks is suppressed for a larger system. This is
because we consider the system part of the system–bath coupling in the form
$\hat{V}_{0}=\sum_{i}^{N}\hat{b}^{\dagger}_{i}\hat{b}_{i}$, and the number of
system–bath couplings, which determines the effective coupling strength for
the system, increases as the number of sites increases.

To confirm this, we depict the two-site result for the case of large
system–bath coupling ($\eta=1$) as the dashed black curve in Fig. 3; as $\eta$
increases, the peaks in the high-frequency region are suppressed, while those
in the low-frequency region are merged and enhanced. As the number of sites
increases, the energy gap decreases because of the increased number of
eigenstates, which is reminiscent of band theory. Thus, the lowest peak in the
four-site case is red-shifted in comparison with that in the two-site case.
This tendency is more pronounced for the periodic boundary condition than for
the open boundary condition, because the effective size of the system is
larger in the former case.

Figure 4: Time evolution of the expectation value of $\hat{H}_{S}$ under
external electric fields for (a) electric excitation and (b) LO phonon
excitation.

####  IV.2.3 Continuous wave response

Finally, we discuss the steady-state response under periodic external electric
fields. Here, we consider the case of off-diagonal system–bath coupling
($\hat{V}_{i}=\hat{b}_{i}^{\dagger}+\hat{b}_{i}$) for (a) the electronic
excitation

$\displaystyle\begin{aligned}
A(t)&=A_{\mathrm{ex}}\cos(\Omega_{\mathrm{ex}}t),\\\\[3.0pt]
E_{\mu}(t)&=0,\end{aligned}$ (13)

and (b) the LO phonon excitation,

$\displaystyle\begin{aligned} A(t)&=0,\\\\[3.0pt] E_{\mu}(t)&=\mu
E_{\mathrm{ex}}\cos(\Omega_{\mathrm{ex}}t).\end{aligned}$ (14)

Under an external time-dependent perturbation, we cannot evaluate the optical
conductivity from Eq. (9), because there is no thermal equilibrium state in
this situation. We therefore simulate Eq. (8) explicitly. For an external
field $A^{\prime}(t)=A_{\mathrm{probe}}\theta(t-t_{\mathrm{on}})$, Eq. (8) is
expressed as (see Appendix B):

$\displaystyle\sigma[\omega]=-\frac{1}{A_{\mathrm{probe}}}\int_{0}^{\infty}dt\,e^{i\omega(t-t_{\mathrm{on}})}[\braket{\hat{j}(t)}-\braket{\hat{j}(t)}_{0}].$
(15)

For sufficiently small $A_{\mathrm{probe}}$ (from $5\times 10^{-3}$ to
$1\times 10^{-2}$), the optical conductivity defined by Eq. (8) agrees with
that defined by Eq. (9) in the equilibrium case.Shao  _et al._ (2016) Because
of the prefactor $e^{-i\omega t_{\mathrm{on}}}$, the profile of the optical
conductivity changes depending on $t_{\mathrm{on}}$, but that dependence is
minor anyway, and so we ignore it in the following discussion.

Here, we consider the high-temperature case ($\beta t_{h}=1$) in order to find
signs of the photoinduced phase transition from insulator to conductor that is
observed in bulk systems. When the HEOM are integrated over a sufficiently
long time, the system driven by the periodic external fields in Eqs. (13) and
(14) reaches a time-dependent steady state. The transient behaviors as
characterized by the system energy are depicted in Fig. 4. In the present
theory of open quantum dynamics, the energy supplied by the external fields
dissipates to the heat bath, whose time-irreversible behavior is characterized
by the fluctuation–dissipation theorem, whereas a wavefunction-based
Schrödinger approach cannot properly treat thermal effects characterized by a
temperature.

Figure 5: Optical conductivity $\operatorname{Re}\\{\sigma[\omega]\\}$ under
the electronic excitation presented in Eq. (13) for (a) weak ($U/t_{h}=1$),
(b) intermediate ($U/t_{h}=3$), and (c) strong ($U/t_{h}=5$) Coulomb
repulsion. The gray curves represent the optical conductivity of the
equilibrium states at high temperature (the same results as shown by the red
curve in Fig. 1). The amplitude $A_{\mathrm{ex}}$ is set to $0.5$.  Figure 6:
Optical conductivity $\operatorname{Re}\\{\sigma[\omega]\\}$ under the LO
phonon excitation presented in Eq. (14) for various values of the excitation
frequency $\Omega_{\mathrm{ex}}$ for (a) weak ($U/t_{h}=1$), (b) intermediate
($U/t_{h}=3$), and (c) strong ($U/t_{h}=5$) Coulomb repulsion. The inverse
temperature is set as $\beta t_{h}=1$. The black lines represent the optical
conductivity of the equilibrium states at high temperature (the same results
as shown by the red curve in Fig. 1). None of the results are normalized. The
amplitude $\mu E_{\mathrm{ex}}$ is set to $0.5$.

In Fig. 5, we present the calculated optical conductivity under the electronic
excitation given by Eq. (13). For a bulk Hubbard system, it has been reported
that by tuning the frequency of the external electronic field to the magnitude
of Coulomb repulsion $U$, the insulating Hubbard model becomes
superconducting.Kaneko  _et al._ (2019) Here, we examine this behavior using
the HH-plus-bath model.

For this purpose, we have to estimate the effective strength of the Coulomb
repulsion in the HH-plus-bath model. In contrast to the Hubbard model, this
estimation is difficult for both the HH and HH-plus-bath cases, because the
transition energies are altered from the Hubbard case owing to the
contribution from electron–phonon coupling. In the strong-repulsion case,
however, the electron–phonon contribution is minor, and we find that the
effective Coulomb repulsion can be evaluated from the peak position of the
optical conductivity, in the same way as in the two-site Hubbard model (see
Appendix C).

Therefore, in Fig. 5(c) in the strong-repulsion case ($U/t_{h}=5$), we apply
an external electric field with frequency $\Omega_{\mathrm{ex}}/t_{h}=0.71$
and find that the peak shifts to the red by approximately $0.11$, which
indicates that the driven system is in the conducting state. This peak shift
can be attributed to $\eta$ pairing, as in the Hubbard case.Kaneko  _et al._
(2019)

In the intermediate-repulsion ($U/t_{h}=3$) and weak-repulsion ($U/t_{h}=1$)
cases, the frequency of the main peak in
$\operatorname{Re}\\{\sigma[\omega]\\}$ does not agree with that given by the
effective Coulomb repulsion, because of the contribution from the LO phonon
interaction. Thus, the peak does not shift, even when the frequency of the
external field is tuned to the dominant peak position, with
$\Omega_{\mathrm{ex}}/t_{h}=0.75$ (weak-repulsion case) and
$\Omega_{\mathrm{ex}}/t_{h}=0.63$ (intermediate-repulsion case). We then vary
the resonant frequency $\Omega_{\mathrm{ex}}$ to find out whether $\eta$
pairing occurs in the weak- and intermediate-repulsion cases, and we set
$\Omega_{\mathrm{ex}}/t_{h}=0.3$. We display
$\operatorname{Re}\\{\sigma[\omega]\\}$ for the weak- and intermediate-
repulsion cases in Figs. 5(a) and 5(b), respectively. In the intermediate-
repulsion case, the peak is red-shifted from $0.13$ to $0.05$, while the peak
profile is similar to that in the strong-repulsion case: We find that $\eta$
pairing occurs in this situation, although estimation of the effective Coulomb
repulsion remains difficult. In the weak-repulsion case, although a peak shift
can be observed, the profiles of the two peaks at $0.05$ and $0.18$ are
different from those in Figs. 5(b) and 5(c). This indicates the suppression of
$\eta$ pairing.

Although our system is small, we may conclude here that $\eta$ pairing, which
plays a key role in the superconductivity of the Hubbard system, is suppressed
when the electron–phonon coupling becomes stronger than the Coulomb repulsion,
while superconductivity may appear when the Coulomb repulsion is stronger or
comparable to the electron–phonon coupling, even when the HH system is coupled
to a heat bath.

In Fig. 6, we present the calculated optical conductivity under the LO phonon
excitation given by Eq. (14) for various values of the laser frequency
$\Omega_{\mathrm{ex}}$ and Coulomb repulsion $U$. In all three cases with
different strengths of repulsion $U$, the position of the lowest peak near
$\omega/t_{h}=0.2$ does not change very much as a function of
$\Omega_{\mathrm{ex}}/t_{h}$, while the peak intensities are suppressed for
$\Omega_{\mathrm{ex}}/t_{h}=2.5$, which is approximately the same as the
effective resonant frequency of the LO phonon mode under off-diagonal
system–bath coupling. In this case, the LO phonon mode is strongly excited by
the external electric field. As a result, a polaron is created, while electron
transfer is suppressed. This indicates that the LO phonon excitation does not
lead to photoinduced superconductivity.

To observe superconductivity with photoexcited phonons, which has been
reported in previous studies, it may be necessary to include Coulomb repulsion
between nearest-neighbor sites Matsueda  _et al._ (2011); Hashimoto and
Ishihara (2017); Seo, Tanaka, and Ishihara (2018) and to change the distance
between molecules.Yonemitsu (2012)

##  V CONCLUDING REMARKS

In this paper, we introduced the HH-plus-bath model to investigate the effects
of a short-range repulsive Coulomb interaction and a local optical vibrational
mode for the investigation of electron dynamics under nonequilibrium
conditions. We employed the HEOM formalism to simulate the dynamics of
electrons described by this model under various physical conditions for local
and nonlocal LO phonon–bath interactions in a numerically rigorous manner. We
calculated the optical conductivity at finite temperature by integrating the
HEOM and found that the present model exhibits conductor-like behavior when
the Coulomb repulsion and the coupling strength between the electrons and the
LO phonons are comparable. This finding agrees with previous results obtained
for an HH bulk system.

We then calculated the optical properties of the system under continuous-wave
laser excitation. We found that the excitation of the electrons promotes
photoinduced conductivity when the Coulomb repulsion is comparable to or
larger than the electron–phonon interaction, while the electrical properties
do not change in the region in which electron–phonon coupling plays a role. We
found that excitation of LO phonons always suppresses the conductivity.
Although we restricted our analysis here to the case of a steady-state
response, there is no inherent restriction on the dynamical simulation with
use of the HEOM. Thus, it is straightforward to compute ultrafast nonlinear
spectra, including two-dimensional electronic–vibrational spectra, on the
basis of the present formulation.

In the present work, we limited our analysis to a small system with a specific
parameter set focusing on the role of local electron–phonon interactions.
Because electron transfer is a long-range effect, an extension of the present
investigation to a larger system would provide deeper insight, in particular,
with regard to studies of superconductivity. Thus, to make the present
approach more useful, further computational efforts need to be made to treat
larger systems that consist of many sites, for example, by employing the
hierarchical Schrödinger equations of motion Nakamura and Tanimura (2018) and
the tensor-train method Borrelli (2019); Borrelli and Dolgov (2021); Shi  _et
al._ (2018); Yan  _et al._ (2021). Such investigations are left for future
work. Nevertheless, we believe that the present results elucidate the key
features of the HH model under thermal conditions with regard to the
fundamental nature of strongly correlated electron systems.

###### Acknowledgements.

The calculations reported here are supported in part by the Research Center
for Computational Science, National Institute of Natural Sciences. Y.T. is
supported by JSPS KAKENHI Grant No. B 21H01884.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A DERIVATION OF THE HEOM FOR THE HH-PLUS-BATH SYSTEM

In this appendix, we construct the HEOM for the HH-plus-bath system. The
reduced density operator of the system,
$\hat{\rho}(t)=\operatorname{tr}_{B}\\{e^{-i\hat{H}_{\mathrm{tot}}t}\hat{\rho}_{\mathrm{tot}}(0)e^{i\hat{H}_{\mathrm{tot}}t}\\}$,
is expressed in the coherent state representation of the path integral as

$\braket{\bm{\xi},\bm{\alpha}}{\hat{\rho}(t)}{\bm{\xi}^{\prime},\bm{\alpha}^{\prime}}=\\\
\int\mathcal{D}^{2}\bm{\xi}\,\mathcal{D}^{2}\bm{\xi}^{\prime}\int\mathcal{D}^{2}\bm{\alpha}\,\mathcal{D}^{2}\bm{\alpha}^{\prime}\int
d^{2}\bm{\xi}_{i}d^{2}\bm{\xi}^{\prime}_{i}\int\frac{d^{2}\bm{\alpha}_{i}}{\pi^{N}}\frac{d^{2}\bm{\alpha}^{\prime}_{i}}{\pi^{N}}e^{iS_{S}(\bm{\xi},\bm{\alpha};t)}\braket{\bm{\xi}_{i},\bm{\alpha}_{i}}{\hat{\rho}(0)}{\bm{\xi}^{\prime}_{i},\bm{\alpha}^{\prime}_{i}}e^{-iS_{S}(\bm{\xi}^{\prime},\bm{\alpha}^{\prime};t)}\mathcal{F}(\bm{\alpha},\bm{\alpha}^{\prime};t).$
(16)

Here, the electronic states are described as a set of Grassmann numbers
$\bm{\xi}=(\xi_{1,\uparrow},\xi_{1,\downarrow},\xi_{2,\uparrow},\ldots,\xi_{N,\downarrow})$,
and $\ket{\bm{\xi}}$ is the eigenstate of fermions (electrons) in the coherent
state representation. The LO phonons are described as a set of complex numbers
$\bm{\alpha}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{N})$, and
$\ket{\bm{\alpha}}$ are the eigenstates of bosons (LO phonons) in the coherent
state representation. The functional $S_{S}(\bm{\xi},\bm{\alpha};t)$ is the
action of the system and is given by

$\displaystyle
S_{S}(\bm{\xi},\bm{\alpha};t)=\int_{0}^{t}dt^{\prime}\left[i\big{\\{}\bm{\xi}^{*}(t^{\prime})\cdot\dot{\bm{\xi}}(t^{\prime})+\tfrac{1}{2}\big{[}\bm{\alpha}^{*}(t^{\prime})\cdot\dot{\bm{\alpha}}(t^{\prime})-\dot{\bm{\alpha}}^{*}(t^{\prime})\cdot\bm{\alpha}(t^{\prime})\big{]}\big{\\}}-\left(H_{S}(\bm{\xi},\bm{\alpha};t^{\prime})+\sum_{i}\lambda_{i}V_{i}^{2}(\alpha_{i};t^{\prime})\right)\right],$
(17)

and $H_{S}(\bm{\xi},\bm{\alpha};t^{\prime})$ and $V_{i}(\alpha_{i};t)$ are the
path-integral representations of $\hat{H}_{S}$ and $\hat{V}_{i}$. The
influence functional $\mathcal{F}(\bm{\alpha},\bm{\alpha}^{\prime};t)$ is
evaluated asTanimura (2014, 2015, 2006)

$\displaystyle\mathcal{F}(\bm{\alpha},\bm{\alpha}^{\prime};t)=\prod_{i}\exp\\!\left(-\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}\,V_{i}^{\times}(\alpha_{i},\alpha_{i}^{\prime};t^{\prime})\\!\left[C^{\prime}_{i}(t^{\prime}-t^{\prime\prime})V_{i}^{\times}(\alpha_{i},\alpha_{i}^{\prime};t^{\prime\prime})-iC^{\prime\prime}_{i}(t^{\prime}-t^{\prime\prime})V_{i}^{\circ}(\alpha_{i},\alpha_{i}^{\prime};t^{\prime\prime})\right]\right).$
(18)

The hyper-operators are defined as
$V_{i}^{\times}(\alpha_{i},\alpha_{i}^{\prime};t)\equiv
V_{i}(\alpha_{i};t)-V_{i}(\alpha_{i}^{\prime};t)$ and
$V_{i}^{\circ}(\alpha_{i},\alpha_{i}^{\prime};t)\equiv
V_{i}(\alpha_{i};t)+V_{i}(\alpha_{i}^{\prime};t)$. The real and imaginary
parts of the bath correlation functions of the $i$th heat bath,
$C^{\prime}_{i}(t)$ and $C^{\prime\prime}_{i}(t)$, are expressed using the SDF
as

$\displaystyle C^{\prime}_{i}(t)$
$\displaystyle=\int_{0}^{\infty}d\omega\,J_{i}(\omega)\coth\\!\left(\frac{\beta\omega}{2}\right)\cos\omega
t,$ (19) $\displaystyle C^{\prime\prime}_{i}(t)$
$\displaystyle=\int_{0}^{\infty}d\omega\,J_{i}(\omega)\sin\omega t,$ (20)

where $\beta$ is the inverse temperature divided by the Boltzmann constant,
$\beta=1/k_{B}T$. By using the Drude–Lorentz SDF in Eq. (6), the
reorganization energy and the two-time correlation function of the heat bath
are characterized by the parameters $\eta_{i}$ and $\gamma_{i}$, as
$\lambda_{i}=\eta_{i}\gamma_{i}/2$,
$C^{\prime}_{i}(t)=c_{i,0}e^{-\gamma_{i}|t|}+\sum_{k}c_{i,k}e^{-\nu_{i,k}|t|}$,
and
$C^{\prime\prime}_{i}(t)=\mathrm{sgn}(t)\eta_{i}\gamma_{i}^{2}e^{-\gamma_{i}|t|}/2$.
The parameters $\\{c_{i,k}\\}$ and $\\{\nu_{i,k}\\}$ are determined by a
Matsubara spectral decomposition scheme Ishizaki and Tanimura (2005) and a
Padé spectral decomposition scheme.Hu  _et al._ (2011) In this paper, we use
the $[N-1/N]$ Padé spectral decomposition scheme.

To obtain differential equations in time, we consider the time derivative
defined as

$\displaystyle\frac{\partial}{\partial t}\hat{\rho}(t)\equiv\int
d^{2}\bm{\xi}\,d^{2}\bm{\xi}^{\prime}\int\frac{d^{2}\bm{\alpha}}{\pi^{N}}\frac{d^{2}\bm{\alpha}^{\prime}}{\pi^{N}}\ket{\bm{\xi},\bm{\alpha}}\lim_{\Delta
t\to 0}\frac{\braket{\bm{\xi},\bm{\alpha}}{\hat{\rho}(t+\Delta
t)}{\bm{\xi}^{\prime},\bm{\alpha}^{\prime}}-\braket{\bm{\xi},\bm{\alpha}}{\hat{\rho}(t)}{\bm{\xi}^{\prime},\bm{\alpha}^{\prime}}}{\Delta
t}\bra{\bm{\xi}^{\prime},\bm{\alpha}^{\prime}}.$ (21)

The details of the evaluation of the right-hand side of Eq. (21) are described
in Refs. Tanimura and Kubo, 1989; Tanimura, 2006, 2020, 2014, 2015. The
auxiliary density operators (ADOs) presented in Eq. (III) are defined as

$\displaystyle\braket{\bm{\xi},\bm{\alpha}}{\hat{\rho}_{\vec{n}_{1},\ldots,\vec{n}_{N}}(t)}{\bm{\xi}^{\prime},\bm{\alpha}^{\prime}}={}$
$\displaystyle\int\mathcal{D}^{2}\bm{\xi}\,\mathcal{D}^{2}\bm{\xi}^{\prime}\int\mathcal{D}^{2}\bm{\alpha}\,\mathcal{D}^{2}\bm{\alpha}^{\prime}\int
d^{2}\bm{\xi}_{i}\,d^{2}\bm{\xi}^{\prime}_{i}\int\frac{d^{2}\bm{\alpha}_{i}}{\pi^{N}}\frac{d^{2}\bm{\alpha}^{\prime}_{i}}{\pi^{N}}$
$\displaystyle\times\prod_{i=1}^{N}\left\\{\left[\int_{0}^{t}dt^{\prime\prime}e^{-\nu_{i,0}(t-t^{\prime\prime})}\left(-ic_{i,0}V_{i}^{\times}(\alpha_{i},\alpha_{i}^{\prime};t^{\prime\prime})-\frac{\eta_{i}\nu_{i,0}^{2}}{2}V^{\circ}(\alpha_{i},\alpha_{i}^{\prime};t^{\prime\prime})\right)\right]^{n_{i,0}}\right.$
$\displaystyle\left.\times\prod_{k=0}^{K}\left(\int_{0}^{t}dt^{\prime\prime}\,e^{-\nu_{i,k}(t-t^{\prime\prime})}\big{[}-ic_{i,k}V_{i}^{\times}(\alpha_{i},\alpha_{i};t^{\prime\prime})\big{]}\right)^{n_{i,k}}\right\\}$
$\displaystyle\times
e^{iS_{S}(\bm{\xi},\bm{\alpha};t)}\braket{\bm{\xi}_{i},\bm{\alpha}_{i}}{\hat{\rho}(0)}{\bm{\xi}^{\prime}_{i},\bm{\alpha}^{\prime}_{i}}e^{-iS_{S}(\bm{\xi}^{\prime},\bm{\alpha}^{\prime};t)}\mathcal{F}(\bm{\alpha},\bm{\alpha}^{\prime};t),$
(22)

in which the zeroth-order term $\hat{\rho}_{\vec{0},\ldots,\vec{0}}(t)$
corresponds to Eq. (16).

Figure 7: Same results as in Fig. 1, but with all of the peaks depicted.

##  Appendix B DERIVATION OF OPTICAL CONDUCTIVITY FOR EQUILIBRIUM STATE

Here, we illustrate the derivation of the optical conductivity for equilibrium
systems, where external fields are set to zero: $A(t)=E_{\mu}(t)=0$. We omit
the superscript $i$ from $t_{h}^{i}$, $U^{i}$, $g^{i}$, and $\Omega^{i}$ for
the sake of simplicity. The current operator is defined as the functional
derivative of the Hamiltonian with respect to the vector potential:

$\displaystyle\hat{j}(t)=-\frac{\delta\hat{H}_{S}(t)}{\delta
A^{\prime}(t)}=i\sum_{i,\sigma}(t_{h}e^{iA^{\prime}(t)}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}-\mathrm{H.c.}).$
(23)

For sufficiently small $A^{\prime}(t)$, we have
$\hat{j}(t)\simeq\hat{j}^{p}-\hat{j}^{d}A^{\prime}(t)$. In the same way, the
system Hamiltonian is expanded as

$\displaystyle\hat{H}_{S}(t)\simeq{}$
$\displaystyle-\sum_{i,\sigma}(t_{h}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i+1,\sigma}+\mathrm{H.c.})+\sum_{i}U\hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}$
$\displaystyle+\sum_{i}g(\hat{n}_{i,\uparrow}+\hat{n}_{i,\downarrow})(\hat{b}_{i}^{\dagger}+\hat{b}_{i})+\sum_{i}\Omega\hat{b}_{i}^{\dagger}\hat{b}_{i}$
$\displaystyle-\hat{j}^{p}A^{\prime}(t)$ $\displaystyle={}$
$\displaystyle\hat{H}^{0}_{S}-\hat{j}^{p}A^{\prime}(t).$ (24)

The total density operator for $\hat{H}_{\mathrm{tot}}$ is then expressed in
the framework of first-order perturbation theory for $A^{\prime}(t)$ as

$\displaystyle\hat{\rho}_{\mathrm{tot}}(t)=\hat{\rho}^{eq}_{\mathrm{tot}}+i\int_{-\infty}^{t}dt^{\prime}\,G_{0}(t-t^{\prime})[\hat{j}^{p},\hat{\rho}^{eq}_{\mathrm{tot}}]A^{\prime}(t^{\prime}).$
(25)

Here,
$G_{0}(t)\hat{O}=e^{-i\hat{H}_{\mathrm{tot}}^{0}t}\hat{O}e^{i\hat{H}_{\mathrm{tot}}^{0}t}$
is the time evolution operator without the perturbation defined as
$\hat{H}_{\mathrm{tot}}^{0}\equiv\hat{H}_{S}^{0}+\hat{H}_{I}+\hat{H}_{B}$. The
expectation value of the current is then expressed in terms of the linear
response function as Mahan (2000); Lenarčič  _et al._ (2014)

$\displaystyle\braket{\hat{j}(t)}$
$\displaystyle=\operatorname{tr}\\{\hat{j}(t)\hat{\rho}^{eq}_{\mathrm{tot}}(t)\\}$
$\displaystyle=\braket{\hat{j}^{p}}_{\mathrm{eq}}-\braket{\hat{j}^{d}}_{\mathrm{eq}}A^{\prime}(t)+\int_{-\infty}^{t}dt^{\prime}R^{(1)}(t-t^{\prime})A^{\prime}(t),$
(26)

where
$R^{(1)}(t)=i\theta(t)\braket{[\hat{j}^{p}(t),\hat{j}^{p}(0)]}_{\mathrm{eq}}$
is the response function expressed in terms of the time-dependent current
operator, defined as
$\hat{j}^{p}(t)=e^{i\hat{H}_{\mathrm{tot}}^{0}t}\hat{j}^{p}e^{-i\hat{H}_{\mathrm{tot}}^{0}t}$,
and
$\braket{\hat{O}}_{\mathrm{eq}}=\operatorname{tr}\\{\hat{O}\hat{\rho}^{eq}_{\mathrm{tot}}\\}$
is the expectation value in the equilibrium state. Without external fields,
the current does not flow and we have $\braket{j^{p}}_{\mathrm{eq}}=0$.

Fourier expressions of $\braket{\hat{j}(t)}$, $A^{\prime}(t)$, and
$R^{(1)}(t)$ are expressed as $\braket{\hat{j}[\omega]}$,
$A^{\prime}[\omega]$, and $R^{(1)}[\omega]$, respectively. Equation (26) is
then expressed as

$\displaystyle\braket{\hat{j}[\omega]}=-\big{(}\\!\braket{\hat{j}^{d}}-R^{(1)}[\omega]\big{)}A^{\prime}[\omega].$

Because the optical conductivity is defined as
$\sigma[\omega]=\braket{\hat{j}[\omega]}/E^{\prime}[\omega]$, and because the
Fourier transform of the electric field $E^{\prime}(t)=-\partial
A^{\prime}(t)/\partial t$ is expressed as $E^{\prime}[\omega]=i\omega
A^{\prime}[\omega]$, we obtain Eq. (9).

The optical conductivity is a response under external perturbation. Thus, we
exclude the equilibrium contribution $\braket{\hat{j}^{p}}_{\mathrm{eq}}$ as
in Eq. (9), and the steady-state contribution $\braket{\hat{j}(t)}_{0}$ as in
Eq. (8). For the direct evaluation of Eq. (8), we introduce an artificial
damping factor $\epsilon$ in the Fourier transform of the vector potential, as
$A^{\prime}[\omega]=\int_{-\infty}^{\infty}dt\,e^{(i\omega+\epsilon)t}A^{\prime}(t)$.
Equation (8) is then modified and takes the following form:

$\displaystyle\sigma[\omega]=\frac{1}{(i\omega+\epsilon)A^{\prime}[\omega]}\int_{0}^{\infty}dt\,e^{i\omega
t}\big{[}\braket{\hat{j}(t)}-\braket{\hat{j}(t)}_{0}\big{]}.$ (27)

When we set $A^{\prime}(t)=A_{\mathrm{probe}}\theta(t-t_{\mathrm{on}})$, we
obtain
$A^{\prime}[\omega]=-A_{\mathrm{probe}}e^{(i\omega+\epsilon)t_{\mathrm{on}}}/(i\omega+\epsilon)$.
By substituting this into Eq. (27), we obtain Eq. (15) for $\epsilon\to 0$.

##  Appendix C ENTIRE PROFILES OF OPTICAL CONDUCTIVITY PRESENTED IN FIG. 1

Figure 7 presents all of the peak profiles of optical conductivity calculated
from Eq. (9) for (a) weak ($U/t_{h}=1$), (b) intermediate ($U/t_{h}=3$), and
(c) strong ($U/t_{h}=5$) Coulomb repulsion in the low-temperature case $\beta
t_{h}=4$ (blue curves) and the high-temperature case $\beta t_{h}=1$ (red
curves). The bath parameters are given by $\gamma_{i}/t_{h}=3$ and
$\eta_{i}=0.05$ $(i=1,2)$. To conduct numerical calculations, we choose
$M_{LO}=9$, $N_{\max}=2$, and $K=2$ for $\beta t_{h}=4$, and $K=1$ for $\beta
t_{h}=1$. The low-frequency part, which is important for characterizing the
conductivity, is presented in Fig. 1.

We find that while the peak positions in the high-frequency regions do not
change very much, regardless of the temperature, those in the low-frequency
region are temperature-dependent, as discussed in Sec. IV.2.1. This is because
these low-frequency peaks arise from the energy gap between the highest
occupied and lowest unoccupied molecular orbitals (HOMO and LUMO). To
illustrate this point, we consider the half-filled two-site Hubbard model
whose eigenenergies are given by $E_{0}=U[1-\sqrt{1+(4t_{h}/U)^{2}}]/2$,
$E_{1}=0$, $E_{2}=U$, and $E_{3}=U[1+\sqrt{1+(4t_{h}/U)^{2}}]/2$. Here,
$E_{1}$ corresponds to the triplet state, whereas the others correspond to the
singlet states. The eigenenergies are depicted in Fig. 8 as functions of $U$.
Because the total spin number
$\hat{S}^{2}=[\sum_{i}(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})/2]^{2}$ is
conserved in our calculations, transitions between the singlet and triplet
states do not occur.

Figure 8: Eigenenergies of half-filled two-site Hubbard model.

In the low-temperature case ($\beta t_{h}=4$), almost all populations are in
the ground state. Thus, the major peak is due to the $E_{0}\to E_{2}$
transition. In the large-$U$ region, the energy gap between these states is
approximately $U$, while the effects from the LO phonons are minor: thus, in
this region, the effective Coulomb repulsion can be estimated from the
position of the maximum peak. Because the contribution of the LO phonons
become larger when the Coulomb repulsion becomes smaller, the peak position
and profile are not easy to predict, as explained by Fig. 5. The side peak
around $0.5$ in the strong-repulsion case also results from the contribution
of the LO phonons. The side peak around $1$ originates from the oscillations
of the LO phonons.

## REFERENCES

## References

  * Saito and Yoshida (2007) G. Saito and Y. Yoshida, Bull. Chem. Soc. Jpn. 80, 1 (2007).
  * Gunnarsson (1997) O. Gunnarsson, Rev. Mod. Phys. 69, 575 (1997).
  * Holstein (1959a) T. Holstein, “Studies of polaron motion: Part i. the molecular-crystal model,” Ann. Phys. 8, 325–342 (1959a).
  * Holstein (1959b) T. Holstein, “Studies of polaron motion: Part ii. the “small” polaron,” Ann. Phys 8, 343–389 (1959b).
  * Beni, Pincus, and Kanamori (1974) G. Beni, P. Pincus, and J. Kanamori, “Low-temperature properties of the one-dimensional polaron band. i. extreme-band-narrowing regime,” Phys. Rev. B 10, 1896 (1974).
  * Aubry (1995) S. Aubry, “Exact results and conjectures on the adiabatic holstein-hubbard model at large electron-phonon coupling,” in  _The Hubbard Model_ (Springer, 1995) pp. 125–143.
  * Takada and Chatterjee (2003) Y. Takada and A. Chatterjee, “Possibility of a metallic phase in the charge-density-wave–spin-density-wave crossover region in the one-dimensional hubbard-holstein model at half filling,” Phys. Rev. B 67, 081102 (2003).
  * Sykora, Hübsch, and Becker (2008) S. Sykora, A. Hübsch, and K. W. Becker, “Coexistence of superconductivity and charge-density waves in a two-dimensional holstein model at half-filling,” (2008), arXiv:0808.4037 .
  * Miyao (2018) T. Miyao, “Ground state properties of the holstein–hubbard model,” Ann. Henri Poincaré 19, 2543–2555 (2018).
  * Berciu (2007) M. Berciu, “Exact green’s functions for the two-site hubbard-holstein hamiltonian,” Phys. Rev. B 75, 081101 (2007).
  * Yamada, Ishizuka, and Ōno (2014) T. Yamada, J. Ishizuka, and Y. Ōno, “Metal–insulator transition and superconductivity in the two-orbital hubbard–holstein model for iron-based superconductors,” J. Phys. Soc. Jpn 83, 044711 (2014).
  * Hebard  _et al._ (1991) A. F. Hebard, M. J. Rosseinsky, R. C. Haddon, D. W. Murphy, S. H. Glarum, T. T. M. Palstra, A. P. Ramirez, and A. R. Kortan, “Superconductivity at 18 k in potassium-doped c60,” Nature 350, 600–601 (1991).
  * Rosseinsky  _et al._ (1991) M. J. Rosseinsky, A. Ramirez, S. Glarum, D. Murphy, R. Haddon, A. Hebard, T. Palstra, A. Kortan, S. Zahurak, and A. Makhija, “Superconductivity at 28 k in rb x c 60,” Phys. Rev. Lett. 66, 2830 (1991).
  * Palstra  _et al._ (1995) T. Palstra, O. Zhou, Y. Iwasa, P. Sulewski, R. Fleming, and B. Zegarski, “Superconductivity at 40k in cesium doped c60,” Solid State Commun. 93, 327–330 (1995).
  * Lanzara  _et al._ (2001) A. Lanzara, P. V. Bogdanov, X. J. Zhou, S. A. Kellar, D. L. Feng, E. D. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, _et al._ , “Evidence for ubiquitous strong electron–phonon coupling in high-temperature superconductors,” Nature 412, 510–514 (2001).
  * Huang  _et al._ (2003) Z. B. Huang, W. Hanke, E. Arrigoni, and D. J. Scalapino, “Electron-phonon vertex in the two-dimensional one-band hubbard model,” Phys. Rev. B 68, 220507 (2003).
  * Johnston  _et al._ (2012) S. Johnston, I. M. Vishik, W. S. Lee, F. Schmitt, S. Uchida, K. Fujita, S. Ishida, N. Nagaosa, Z. X. Shen, and T. P. Devereaux, “Evidence for the importance of extended coulomb interactions and forward scattering in cuprate superconductors,” Phys. Rev. Lett. 108, 166404 (2012).
  * Ishihara and Nagaosa (2004) S. Ishihara and N. Nagaosa, “Interplay of electron-phonon interaction and electron correlation in high-temperature superconductivity,” Phys. Rev. B 69, 144520 (2004).
  * Mishchenko  _et al._ (2008) A. S. Mishchenko, N. Nagaosa, Z.-X. Shen, G. De Filippis, V. Cataudella, T. P. Devereaux, C. Bernhard, K. W. Kim, and J. Zaanen, “Charge dynamics of doped holes in high t c cuprate superconductors: a clue from optical conductivity,” Phys. Rev. Lett. 100, 166401 (2008).
  * De Filippis  _et al._ (2009) G. De Filippis, V. Cataudella, A. S. Mishchenko, C. A. Perroni, and N. Nagaosa, “Optical conductivity of a doped mott insulator: The interplay between correlation and electron-phonon interaction,” Phys. Rev. B 80, 195104 (2009).
  * Takada (1996) Y. Takada, “Superconductivity in the half-filled hubbard-holstein model in the antiadiabatic region,” J. Phys. Soc. Jpn. 65, 1544–1547 (1996).
  * Chatterjee (2010) A. Chatterjee, “Existence of an intermediate metallic phase at the sdw-cdw crossover region in the one-dimensional holstein-hubbard model at half-filling,” Adv. Cond. Matter Phys. 2010 (2010).
  * Tezuka, Arita, and Aoki (2007) M. Tezuka, R. Arita, and H. Aoki, “Phase diagram for the one-dimensional hubbard-holstein model: a density-matrix renormalization group study,” Phys. Rev. B 76 (2007).
  * Von der Linden, Berger, and Valášek (1995) W. Von der Linden, E. Berger, and P. Valášek, “The hubbard-holstein model,” J. Low. Temp. Phys. 99, 517–525 (1995).
  * Macridin, Sawatzky, and Jarrell (2004) A. Macridin, G. A. Sawatzky, and M. Jarrell, “Two-dimensional hubbard-holstein bipolaron,” Phys. Rev. B 69, 245111 (2004).
  * Clay and Hardikar (2005) R. T. Clay and R. P. Hardikar, “Intermediate phase of the one dimensional half-filled hubbard-holstein model,” Phys. Rev. Lett. 95, 096401 (2005).
  * Nowadnick  _et al._ (2012) E. A. Nowadnick, S. Johnston, B. Moritz, R. T. Scalettar, and T. P. Devereaux, “Competition between antiferromagnetic and charge-density-wave order in the half-filled hubbard-holstein model,” Phys. Rev. Lett. 109, 246404 (2012).
  * Johnston  _et al._ (2013) S. Johnston, E. A. Nowadnick, Y. F. Kung, B. Moritz, R. T. Scalettar, and T. P. Devereaux, “Determinant quantum monte carlo study of the two-dimensional single-band hubbard-holstein model,” Phys. Rev. B 87, 235133 (2013).
  * Ohgoe and Imada (2017) T. Ohgoe and M. Imada, “Competition among superconducting, antiferromagnetic, and charge orders with intervention by phase separation in the 2d holstein-hubbard model,” Phys. Rev. Lett. 119, 197001 (2017).
  * Weber and Hohenadler (2018) M. Weber and M. Hohenadler, “Two-dimensional holstein-hubbard model: Critical temperature, ising universality, and bipolaron liquid,” Phys. Rev. B 98, 085405 (2018).
  * Tsuji, Oka, and Aoki (2009) N. Tsuji, T. Oka, and H. Aoki, “Nonequilibrium steady state of photoexcited correlated electrons in the presence of dissipation,” Phys. Rev. Lett. 103, 047403 (2009).
  * Mitrano  _et al._ (2016) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, _et al._ , “Possible light-induced superconductivity in k 3 c 60 at high temperature,” nature 530, 461–464 (2016).
  * Mazza and Georges (2017) G. Mazza and A. Georges, “Nonequilibrium superconductivity in driven alkali-doped fullerides,” Phys. Rev. B 96, 064515 (2017).
  * Sentef  _et al._ (2016) M. A. Sentef, A. Kemper, A. Georges, and C. Kollath, “Theory of light-enhanced phonon-mediated superconductivity,” Phys. Rev. B 93, 144506 (2016).
  * Kennes  _et al._ (2017) D. M. Kennes, E. Y. Wilner, D. R. Reichman, and A. J. Millis, “Transient superconductivity from electronic squeezing of optically pumped phonons,” Nat. Phys. 13, 479 (2017).
  * Kaneko  _et al._ (2019) T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, “Photoinduced $\eta$ pairing in the hubbard model,” Phys. Rev. Lett. 122, 077002 (2019).
  * Ido, Ohgoe, and Imada (2017) K. Ido, T. Ohgoe, and M. Imada, “Correlation-induced superconductivity dynamically stabilized and enhanced by laser irradiation,” Sci. Adv. 3, e1700718 (2017).
  * Kim  _et al._ (2016) M. Kim, Y. Nomura, M. Ferrero, P. Seth, O. Parcollet, and A. Georges, “Enhancing superconductivity in a 3 c 60 fullerides,” Phys. Rev. B 94, 155152 (2016).
  * Werner and Eckstein (2015) P. Werner and M. Eckstein, “Field-induced polaron formation in the holstein-hubbard model,” EPL (Europhysics Letters) 109, 37002 (2015).
  * Golež  _et al._ (2012) D. Golež, J. Bonča, L. Vidmar, and S. A. Trugman, “Relaxation dynamics of the holstein polaron,” Phys. Rev. Lett. 109, 236402 (2012).
  * Cantaluppi  _et al._ (2018) A. Cantaluppi, M. Buzzi, G. Jotzu, D. Nicoletti, M. Mitrano, D. Pontiroli, M. Riccò, A. Perucchi, P. Di Pietro, and A. Cavalleri, “Pressure tuning of light-induced superconductivity in k 3 c 60,” Nat. Phys. 14, 837–841 (2018).
  * Nicoletti  _et al._ (2017) D. Nicoletti, E. Casandruc, D. Fu, P. Giraldo-Gallo, I. R. Fisher, and A. Cavalleri, “Anomalous relaxation kinetics and charge-density-wave correlations in underdoped bapb1- xbixo3,” Proc. Natl. Acad. Sci. U.S.A. 114, 9020–9025 (2017).
  * Mankowsky  _et al._ (2017) R. Mankowsky, M. Fechner, M. Först, A. von Hoegen, J. Porras, T. Loew, G. L. Dakovski, M. Seaberg, S. Möller, G. Coslovich, _et al._ , “Optically induced lattice deformations, electronic structure changes, and enhanced superconductivity in yba2cu3o6. 48,” Struct. Dyn. 4, 044007 (2017).
  * Weiss (2012) U. Weiss,  _Quantum Dissipative Systems_ (World Scientific, Singapore, 2012).
  * Garg, Onuchic, and Ambegaokar (1985) A. Garg, J. N. Onuchic, and V. Ambegaokar, “Effect of friction on electron transfer in biomolecules,” J. Chem. Phys. 83, 4491–4503 (1985).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, “Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature,” Proc. Natl. Acad. Sci. U.S.A. 106, 17255 (2009).
  * Hein  _et al._ (2012) B. Hein, C. Kreisbeck, T. Kramer, and M. Rodríguez, “Modelling of oscillations in two-dimensional echo-spectra of the fenna–matthews–olson complex,” New Journal of Physics 14, 023018 (2012).
  * Tanimura and Kubo (1989) Y. Tanimura and R. Kubo, “Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath,” Journal of the Physical Society of Japan 58, 101–114 (1989), https://doi.org/10.1143/JPSJ.58.101 .
  * Tanimura (2006) Y. Tanimura, “Stochastic liouville, langevin, fokker–planck, and master equation approaches to quantum dissipative systems,” Journal of the Physical Society of Japan 75, 082001 (2006), https://doi.org/10.1143/JPSJ.75.082001 .
  * Tanimura (2020) Y. Tanimura, “Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM),” J. Chem. Phys. 153, 020901 (2020).
  * Chen, Zhao, and Tanimura (2015) L. Chen, Y. Zhao, and Y. Tanimura, “Dynamics of a one-dimensional holstein polaron with the hierarchical equations of motion approach,” J. Phys. Chem. Lett. 6, 3110–3115 (2015).
  * Zhou  _et al._ (2016) N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, “Fast, accurate simulation of polaron dynamics and multidimensional spectroscopy by multiple davydov trial states,” The Journal of Physical Chemistry A 120, 1562–1576 (2016), pMID: 26871592, https://doi.org/10.1021/acs.jpca.5b12483 .
  * Dunn, Tempelaar, and Reichman (2019) I. S. Dunn, R. Tempelaar, and D. R. Reichman, “Removing instabilities in the hierarchical equations of motion: Exact and approximate projection approaches,” J. Chem. Phys. 150, 184109 (2019).
  * Cainelli and Tanimura (2021) M. Cainelli and Y. Tanimura, “Exciton transfer in organic photovoltaic cells: A role of local and nonlocal electron-phonon interactions in a donor domain,” J. Chem. Phys. 154, 034107 (2021).
  * Ueno and Tanimura (2021) S. Ueno and Y. Tanimura, “Modeling and simulating the excited-state dynamics of a system with condensed phases: A machine learning approach,” J. Chem. Theory Comput. 17, 3618 (2021).
  * Tanimura (1990) Y. Tanimura, “Nonperturbative expansion method for a quantum system coupled to a harmonic-oscillator bath,” Phys. Rev. A 41, 6676–6687 (1990).
  * Tanaka and Tanimura (2009) M. Tanaka and Y. Tanimura, “Quantum dissipative dynamics of electron transfer reaction system: Nonperturbative hierarchy equations approach,” J. Phys. Soc. Jpn. 78, 073802 (2009).
  * Tanaka and Tanimura (2010) M. Tanaka and Y. Tanimura, “Multistate electron transfer dynamics in the condensed phase: Exact calculations from the reduced hierarchy equations of motion approach,” J. Chem. Phys. 132, 214502 (2010).
  * Tanimura (2012) Y. Tanimura, “Reduced hierarchy equations of motion approach with drude plus brownian spectral distribution: Probing electron transfer processes by means of two-dimensional correlation spectroscopy,” J. Chem. Phys. 137, 22A550 (2012).
  * Tanimura (2014) Y. Tanimura, “Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities,” The Journal of Chemical Physics 141, 044114 (2014), https://doi.org/10.1063/1.4890441 .
  * Tanimura (2015) Y. Tanimura, “Real-time and imaginary-time quantum hierarchal fokker-planck equations,” The Journal of Chemical Physics 142, 144110 (2015), https://doi.org/10.1063/1.4916647 .
  * Zhang  _et al._ (2017) H.-D. Zhang, Q. Qiao, R.-X. Xu, X. Zheng, and Y. Yan, “Efficient steady-state solver for hierarchical quantum master equations,” The J. Chem. Phys. 147, 044105 (2017).
  * Mahan (2000) G. D. Mahan,  _Many-particle physics_ (Kluwer Academic/Plenum Publishers, New York, 2000).
  * Lenarčič  _et al._ (2014) Z. Lenarčič, D. Golež, J. Bonča, and P. Prelovšek, “Optical response of highly excited particles in a strongly correlated system,” Phys. Rev. B 89, 125123 (2014).
  * Bergeron  _et al._ (2011) D. Bergeron, V. Hankevych, B. Kyung, and A.-M. S. Tremblay, “Optical and dc conductivity of the two-dimensional hubbard model in the pseudogap regime and across the antiferromagnetic quantum critical point including vertex corrections,” Phys. Rev. B 84, 085128 (2011).
  * Castella, Zotos, and Prelovšek (1995) H. Castella, X. Zotos, and P. Prelovšek, “Integrability and ideal conductance at finite temperatures,” Phys. Rev. Lett. 74, 972 (1995).
  * Zotos and Prelovšek (1996) X. Zotos and P. Prelovšek, “Evidence for ideal insulating or conducting state in a one-dimensional integrable system,” Phys. Rev. B 53, 983 (1996).
  * Okumura and Tanimura (1997) K. Okumura and Y. Tanimura, “Two-time correlation functions of a harmonic system nonbilinearly coupled to a heat bath: Spontaneous raman spectroscopy,” Phys. Rev. E 56, 2747–2750 (1997).
  * Ishizaki and Tanimura (2006) A. Ishizaki and Y. Tanimura, “Modeling vibrational dephasing and energy relaxation of intramolecular anharmonic modes for multidimensional infrared spectroscopies,” The Journal of Chemical Physics 125, 084501 (2006), https://doi.org/10.1063/1.2244558 .
  * Shao  _et al._ (2016) C. Shao, T. Tohyama, H.-G. Luo, and H. Lu, “Numerical method to compute optical conductivity based on pump-probe simulations,” Phys. Rev. B 93, 195144 (2016).
  * Matsueda  _et al._ (2011) H. Matsueda, S. Sota, T. Tohyama, and S. Maekawa, “Relaxation dynamics of photocarriers in one-dimensional mott insulators coupled to phonons,” J. Phys. Soc. Jpn. 81, 013701 (2011).
  * Hashimoto and Ishihara (2017) H. Hashimoto and S. Ishihara, “Photoinduced charge-order melting dynamics in a one-dimensional interacting holstein model,” Phys. Rev. B 96, 035154 (2017).
  * Seo, Tanaka, and Ishihara (2018) H. Seo, Y. Tanaka, and S. Ishihara, “Photoinduced collective mode, inhomogeneity, and melting in a charge-order system,” Phys. Rev. B 98, 235150 (2018).
  * Yonemitsu (2012) K. Yonemitsu, “Theory of photoinduced phase transitions in molecular conductors: Interplay between correlated electrons, lattice phonons and molecular vibrations,” Crystals 2, 56–77 (2012).
  * Nakamura and Tanimura (2018) K. Nakamura and Y. Tanimura, “Hierarchical schrödinger equations of motion for open quantum dynamics,” Phys. Rev. A 98, 012109 (2018).
  * Borrelli (2019) R. Borrelli, “Density matrix dynamics in twin-formulation: An efficient methodology based on tensor-train representation of reduced equations of motion,” The Journal of Chemical Physics 150, 234102 (2019), https://doi.org/10.1063/1.5099416 .
  * Borrelli and Dolgov (2021) R. Borrelli and S. Dolgov, “Expanding the range of hierarchical equations of motion by tensor-train implementation,” The Journal of Physical Chemistry B 125, 5397–5407 (2021), pMID: 33980011, https://doi.org/10.1021/acs.jpcb.1c02724 .
  * Shi  _et al._ (2018) Q. Shi, Y. Xu, Y. Yan, and M. Xu, “Efficient propagation of the hierarchical equations of motion using the matrix product state method,” The Journal of Chemical Physics 148, 174102 (2018), https://doi.org/10.1063/1.5026753 .
  * Yan  _et al._ (2021) Y. Yan, M. Xu, T. Li, and Q. Shi, “Efficient propagation of the hierarchical equations of motion using the tucker and hierarchical tucker tensors,” The Journal of Chemical Physics 154, 194104 (2021), https://doi.org/10.1063/5.0050720 .
  * Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, “Quantum dynamics of system strongly coupled to low-temperature colored noise bath: Reduced hierarchy equations approach,” J. Phys. Soc. Jpn. 74, 3131–3134 (2005).
  * Hu  _et al._ (2011) J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y.-J. Yan, J. Chem. Phys. 134, 244106 (2011).