# Linear and non-linear infrared response of one-dimensional vibrational
Holstein polarons in the anti-adiabatic limit: optical and acoustical phonon
models

Cyril Falvo<EMAIL_ADDRESS>Institut des Sciences Moléculaires d’Orsay
(ISMO), CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91405 Orsay, France
Univ. Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France

###### Abstract

The theory of linear and non-linear infrared response of vibrational Holstein
polarons in one-dimensional lattices is presented in order to identify the
spectral signatures of self-trapping phenomena. Using a canonical
transformation the optical response is computed from the small polaron point
of view which is valid in the anti-adiabatic limit. Two types of phonon baths
are considered: optical phonons and acoustical phonons, and simple expressions
are derived for the infrared response. It is shown that for the case of
optical phonons, the linear response can directly probe the polaron density of
states. The model is used to interpret the experimental spectrum of crystaline
actetanilide in the C$=$O range. For the case of acoustical phonons, it is
shown that two bound states can be observed in the two-dimensional infrared
spectrum at low temperature. At high temperature, analysis of the time-
dependence of the two-dimensional infrared spectrum indicates that bath
mediated correlations slow down spectral diffusion. The model is used to
interpret the experimental linear-spectroscopy of model $\alpha$-helix and
$\beta$-sheet polypeptides. This work shows that the Davydov Hamiltonian
cannot explain the observations in the NH stretching range.

##  I Introduction

The dynamics of electronic or vibrational excitons in quasi one-dimensional
lattices has been an open topic of research for the past 60 years.Mahan
(1981); Holstein (1959a, b) The interplay between exciton delocalization and
coupling with the lattice vibrations results in the self-trapping phenomena,
i.e. the formation of a polaron. A polaron usually refers to a quasi-particle
that comprises the exciton and the lattice deformation created by the exciton
which modifies its dynamics.Mahan (1981); Holstein (1959a, b) The concept of
self-trapping in one-dimensional lattices has a large number of applications
for example in molecular aggregates,Spano (2010); Huynh  _et al._ (2013); Lu
and Mukamel (1991); Sun  _et al._ (2015); Chorošajev  _et al._ (2014);
Chorosajev, Rancova, and Abramavicius (2016) conjugated polymers,Yamagata and
Spano (2014); Barford and Marcus (2014) halogen-bridged metal
complexes,Okamoto  _et al._ (1992) molecular crystalsFillaux (1981); Barthes
_et al._ (1998); Herrebout, Clou, and Desseyn (2001); Careri  _et al._ (1983,
1984); Eilbeck, Lomdahl, and Scott (1984); Alexander and Krumhansl (1986);
Edler, Hamm, and Scott (2002); Edler and Hamm (2002, 2003); Hamm and Edler
(2006) and biological macromolecules.Davydov (1973, 1985); Scott (1982, 1992);
Edler  _et al._ (2004, 2005); Hamm (2009); Brown and Ivić (1989); Ivić  _et
al._ (1997); Pouthier (2003); Pouthier and Falvo (2004); Falvo and Pouthier
(2005a, b, c); Tsivlin, Meyer, and May (2006); Tsivlin and May (2006, 2007);
Bodis  _et al._ (2009); Cruzeiro (2009); Goj and Bittner (2011) Vibrational
excitons emerge in molecular crystal or within biological macromolecules by
the delocalization of high-frequency vibrations through dipole-dipole
interactions. This is the case for example, in $\alpha$-helix polypeptides
where the amide-I band corresponds to the delocalization of the C$=$O
vibrations of each peptide group along the backbone of the peptide that forms
a quasi one-dimensional lattice.Miyazawa (1960) These excitons are strongly
coupled to the CO$\cdots$NH hydrogen-bonds that stabilize the helix. This
coupling was first introduced by A. S. Davydov which speculated the formation
of a soliton able to transfer energy from one side of the $\alpha$-helix to
the other.Davydov (1973, 1985); Scott (1982, 1992) It appeared later that this
coupling results into the formation of a vibrational polaron rather than a
soliton.Brown and Ivić (1989); Ivić  _et al._ (1997); Pouthier (2003);
Pouthier and Falvo (2004); Falvo and Pouthier (2005a, b, c)

Just a few years after the work of Davydov, the infrared (IR) spectroscopy of
the molecular crystals acetanilide (ACN)Careri  _et al._ (1983, 1984) and
N-methylacetamide (NMA)Fillaux (1981); Barthes  _et al._ (1998); Herrebout,
Clou, and Desseyn (2001) showed some anomalous temperature dependance that was
interpreted as a signature of self-trapping. These molecular crystals which
consist of quasi-one-dimensional chains of hydrogen-bonded peptide groups
resembling the hydrogen-bond network of $\alpha$-helices were then considered
as model systems for polypeptides. In ACN, at ambiant temperature, the amide-I
band is characterized by a single band located at 1666 $\textrm{cm}^{-1}$,
while at lower temperature a second band appears at 1650 $\textrm{cm}^{-1}$.
The amide-A band corresponding to the N$-$H stretching vibration is
characterized by a main band located at 3295 $\textrm{cm}^{-1}$ with a series
of 9 satelite peaks towards low frequency. These observations show that a
strong coupling occurs between the C$=$O and N$-$H vibrations with some low
frequency optical phonons. From a theoretical point of view the dynamics of
vibrational excitons in molecular crystals and in $\alpha$-helices can be
described by the same Holstein Hamiltonian, which was first introduced to
describe the dynamics of electrons in molecular crystals. The main difference
is that in molecular crystals the vibrational excitons are coupled to optical
phonons while in the Davydov Hamiltonian the vibrational excitons are coupled
to acoustical phonons.Eilbeck, Lomdahl, and Scott (1984); Alexander and
Krumhansl (1986); Scott (1982, 1992)

Two decades after these observations, Hamm and Edler shed new light on the
dynamics of vibrational polarons by performing non-linear IR spectroscopy of
ACN and NMA crystals as well as a model $\alpha$-helix.Edler, Hamm, and Scott
(2002); Edler and Hamm (2002, 2003); Hamm and Edler (2006); Edler  _et al._
(2004, 2005) This work was reviewed in Ref. 30. Developed over the past two-
decades, time-resolved nonlinear IR spectroscopy, in particular two-
dimensional (2D) IR spectroscopy have allowed researchers to study the
vibrational dynamics of condensed phase systems including peptides, proteins,
and water.Fayer (2013); Hamm and Zanni (2011); Khalil, Demirdöven, and
Tokmakoff (2003); Loparo, Roberts, and Tokmakoff (2006); Wong  _et al._
(2013); Bloem  _et al._ (2012); Middleton  _et al._ (2012); Bandaria  _et al._
(2010); Ghosh  _et al._ (2014); Kim and Hochstrasser (2009) 2D-IR spectroscopy
can probe vibrational anharmonic couplings, vibrational relaxation, population
transport, chemical-exchange dynamics and spectral diffusion therefore
providing much more information than absorption spectroscopy.Jansen and
Knoester (2009); Zheng  _et al._ (2005); Falvo  _et al._ (2008); Cho (2008);
Kim and Hochstrasser (2009) In ACN and NMA, Edler and Hamm used 2D-IR
spectroscopy to show that vibrational self-trapping and the formation of
vibrational polarons occured in these molecular crystals.Edler, Hamm, and
Scott (2002); Edler and Hamm (2002, 2003); Hamm and Edler (2006); Hamm (2009)
They also performed pump-probe spectroscopy on a model $\alpha$-helix in the
N$-$H spectral range.Edler  _et al._ (2004, 2005) They show the appearance in
the two-exciton spectrum of two bound states. These two bound states were
interpreted as the signature of a strong coupling between the vibrational
exciton and a set of accoustical phonons in accordance with the original model
of Davydov. A similar observation was made a few years later in the spectrum
of a model $\beta$-sheet peptide.Bodis  _et al._ (2009)

A large number of theoretical studies have been devoted to the Holstein
Hamiltonian (electronic or vibrational). Holstein polarons are usually
described within two limiting cases that depend on the size of the polaron
wavefunction: small and large polarons.Holstein (1959a, b) For the former the
discreteness of the lattice plays a key role while for the later a continuum
approximation can be used. Large polarons are often described within the
adiabatic limit, i.e. the case when the lattice deformation remains static
when the exciton moves along the lattice. In this case, variational
approachesLu and Mukamel (1991); Huynh  _et al._ (2013); Sun  _et al._ (2015)
or mixed quantum-classical simulationsScott (1992); Cruzeiro (2009) gives in
general good results.Brown and Ivić (1989); Ivić  _et al._ (1997) In contrast,
small polarons are usually described within the anti-adiabatic limit which
corresponds to a weak exciton coupling.Tempelaar  _et al._ (2013) In this
limit the lattice deformation follows the exciton modifying its effective
mass. In this case, a canonical transformation allows to switch to the small
polaron point of view.Lang and Firsov (1963); Brown and Ivić (1989); Ivić  _et
al._ (1997); Pouthier (2003); Pouthier and Falvo (2004); Falvo and Pouthier
(2005a, b, c); Yalouz, Falvo, and Pouthier (2017) It has been shown that for
the case of vibrational excitons, because the hopping constant between
nearest-neighbor lattice sites is in general small compared to the phonon
frequency, the anti-adiabatic limit gives good results.Brown and Ivić (1989);
Ivić  _et al._ (1997); Pouthier (2003) The Holstein Hamiltonian has been also
solved by a variety of numerical methods that include the Density Matrix
Renormalization GroupJeckelmann and White (1998), the Multi-Configuration
Time-dependant Hartree method,Tsivlin and May (2007) the Hierarchy Equation Of
Motion (HEOM)Chen, Zhao, and Tanimura (2015), using the two-particle
approximationPhilpott (1971); Spano (2002) or using a direct exact
diagonalization.Hamm and Edler (2006); Yalouz, Falvo, and Pouthier (2017);
Yalouz, Pouthier, and Falvo (2017) Most theoretical studies focused on the
energy transport properties of polarons and on linear spectroscopy, very few
are dedicated to predict non-linear spectroscopy. This is particularly true
for the case of vibrational polarons where to my knowledge the few studies
were conducted on pump-probe spectroscopy,Edler  _et al._ (2004); Tsivlin and
May (2006); Tsivlin, Meyer, and May (2006); Woutersen (2007) and none were
conducted on 2D-IR spectroscopy. Note that two recent studies focused on the
2D spectroscopy of electronic excitons in molecular aggregates using a
variational approach.Huynh  _et al._ (2013); Sun  _et al._ (2015) However, as
mentioned earlier this approach is mostly adapted for the case of large
polarons and is not adapted for vibrational excitons. Therefore, for
vibrational polarons there is a clear lack of theoretical work to predict the
non-linear IR response.

In this article, the theory of linear and non-linear spectroscopy of
vibrational polarons in one-dimensional lattices is presented in order to
establish a physical framework to identify the spectral signatures of self-
trapping phenomena. The case of both optical and acoustical phonons are
considered allowing to cover both Davydov and molecular crystals models. This
theoretical work relies on the anti-adiabatic approximation which assumes that
the vibrational excitons are slow compared to the phonon bath. Note that in
this article, the simple case of a one-dimensional (1D) lattice is
investigated keeping the Holstein Hamiltonian as simple as possible in order
to set up the framework for the nonlinear response of vibrational polarons and
present analytical results. In section II, simple expressions are derived for
the linear and non-linear optical response of vibrational polarons. These
expressions are used in the section III for a variety of parameters values in
the Holstein model. In section IV where further theoretical derivations are
performed, the model results are discussed within the context of experimental
observations. Finally, future experiments to probe self-trapping phenomena are
suggested in addition to theoretical developments needed in the future as well
as conclusions are presented in section V.

##  II Theoretical model

In this section, the theoretical framework describing vibrational excitons
coupled to optical and acoustical phonons is presented within the context of
the anti-adiabatic limit. Using this approximation the linear and third-order
response is given.

###  II.1 Vibrational Holstein Hamiltonian

A one-dimensional chain of $N$ identical high frequency vibrations coupled to
a phonon bath is considered. The system Hamiltonian $\hat{H}$ is written as

$\hat{H}=\hat{H}_{v}+\hat{H}_{b}+\hat{H}_{vb},$ (1)

where $\hat{H}_{v}$ is the vibrational Hamiltonian, $\hat{H}_{b}$ the bath
Hamiltonian and $\hat{H}_{vb}$ the coupling between the vibrations and the
bath. The vibrational Hamiltonian $\hat{H}_{v}$ is described by an excitonic
Hamiltonian written as

$\hat{H}_{v}=\sum_{n}\omega_{0}b_{n}^{\dagger}b_{n}-Ab_{n}^{\dagger
2}b_{n}^{2}+J\left(b_{n+1}^{\dagger}+b_{n-1}^{\dagger}\right)b_{n},$ (2)

where $\omega_{0}$ is the fundamental frequency, $A$ is the anharmonicity, $J$
is the hopping constant and where $b_{n}^{\dagger}$ and $b_{n}$ are the vibron
creation and annihilation operators. In Eq. (2) and in the remaining of this
paper, the convention $\hbar=1$ is used. The bath is described by a set of $N$
phonons of frequencies $\Omega_{q}$ and wavevector $q=2\pi p/N$ with
$p=-(N-1)/2,\dots,(N-1)/2$. Using the phonon creation and annihilation
operators $a^{\dagger}_{q}$ and $a_{q}$ the bath Hamiltonian is written

$\hat{H}_{b}=\sum_{q}\Omega_{q}(a_{q}^{\dagger}a_{q}+1/2).$ (3)

To describe the coupling between the high frequency vibrations and the bath
modes it is assumed that each bath mode induces fluctuations of the
fundamental frequencies. The coupling hamiltonian is then written as

$\hat{H}_{vp}=\frac{1}{\sqrt{N}}\sum_{n}\sum_{q}\left(\Delta_{q}e^{-\textrm{i}qn}a_{q}^{\dagger}+\Delta_{q}^{*}e^{\textrm{i}qn}a_{q}\right)b_{n}^{\dagger}b_{n}.$
(4)

Note that here each bath modes are coupled to all the vibrations therefore
introducing strong bath mediated correlations between different vibrations. In
this article, two types of phonon models are considered, a model of optical
phonons and a model of acoustical phonons. Derivation of the optical and
acoustical models are detailed in appendix A. For the optical phonon model the
phonon frequency and coupling are written as

$\displaystyle\Omega^{\text{opt}}_{q}=\Omega_{\text{opt}},$ (5)
$\displaystyle\Delta^{\text{opt}}_{q}=\Delta_{\text{opt}},$ (6)

where $\Omega_{\text{opt}}$ is the frequency of the phonon and
$\Delta_{\text{opt}}$ is the coupling strength. The acoustical phonon model is
derived from the Davydov Hamiltonian and is given by the following parameters

$\displaystyle\Omega_{q}^{\text{ac}}=\Omega_{\text{ac}}\left|\sin q/2\right|,$
(7)
$\displaystyle\Delta^{\text{ac}}_{q}=-2\textrm{i}\Delta_{{\text{ac}}}\sqrt{|\sin
q/2|}\cos q/2,$ (8)

where $\Omega_{\text{ac}}$ is the cutoff frequency and where
$\Delta_{\text{ac}}$ is the coupling strength.

###  II.2 Effective Hamiltonian in the anti adiabatic limit

To partially remove the vibron-bath coupling Hamiltonian, a Lang-Firsov
transformation is applied.Lang and Firsov (1963) A “full dressing” is
considered and the following unitary transformation is introduced

$\hat{U}=\exp\left(\sum_{n}\hat{X}_{n}b_{n}^{\dagger}b_{n}\right),$ (9)

where the operator $\hat{X}_{n}$ is defined by

$\hat{X}_{n}=\frac{1}{\sqrt{N}}\sum_{q}\left(\frac{\Delta_{q}e^{-\textrm{i}qn}}{\Omega_{q}}a_{q}^{\dagger}-\frac{\Delta_{q}^{*}e^{\textrm{i}qn}}{\Omega_{q}}a_{q}\right).$
(10)

By using Eq. (9), the transformed Hamiltonian $\tilde{H}=U\hat{H}U^{\dagger}$
is written as

$\tilde{H}=\sum_{n}\left(\omega_{0}-\epsilon_{0}\right)b_{n}^{\dagger}b_{n}-\left(A+\epsilon_{0}\right)b_{n}^{\dagger
2}b_{n}^{2}-2\sum_{n<m}\epsilon_{|n-m|}b_{n}^{\dagger}b_{n}b_{m}^{\dagger}b_{m}\\\
+\sum_{n}J\left(\Theta^{\dagger}_{n+1}b_{n+1}^{\dagger}+\Theta^{\dagger}_{n-1}b_{n-1}^{\dagger}\right)\Theta_{n}b_{n}+\hat{H}_{b},$
(11)

where the dressing operators $\Theta^{\dagger}_{n}$ are defined by the
transformation of the vibron creation operators $b^{\dagger}_{n}$

$\hat{U}b^{\dagger}_{n}\hat{U}^{\dagger}=b^{\dagger}_{n}\Theta^{\dagger}_{n},$
(12)

and are written as

$\Theta^{\dagger}_{n}=\exp\left(\hat{X}_{n}\right).$ (13)

The parameters $\epsilon_{n}$ characterize the reorganizational energies of
the bath, they are written as

$\epsilon_{n}=\frac{1}{N}\sum_{q}\frac{\left|\Delta_{q}\right|^{2}}{\Omega_{q}}\cos(nq).$
(14)

In the small polaron point of view, the vibrational exciton are dressed by the
bath. The remaining coupling between the polaron and the bath now operates
through the hopping term which is modulated by bath coherent states. The main
advantage of this procedure is that the exciton-phonon coupling has been
strongly reduced and a mean field approach can then be used.Ivić  _et al._
(1997) The final Hamiltonian is written as a sum of three contribution

$\tilde{H}=\hat{H}_{0}+\hat{H}_{b}+\Delta\hat{H},$ (15)

where $\hat{H}_{0}=\langle\tilde{H}-\hat{H}_{b}\rangle_{b}$ is the effective
Hamiltonian of the dressed excitons and
$\Delta\hat{H}=\tilde{H}-\hat{H}_{b}-\hat{H}_{0}$ is the remaining part of the
exciton-bath interaction. The symbol $\langle\dots\rangle_{b}$ stands for the
thermal average over the bath degrees of freedom which are assumed to be in
equilibrium at temperature $T$. After straightforward calculation, the
effective polaron hamiltonian is finally written as

$\hat{H}_{0}=\sum_{n}\left(\omega_{0}-\epsilon_{0}\right)b_{n}^{\dagger}b_{n}-\left(A+\epsilon_{0}\right)b_{n}^{\dagger
2}b_{n}^{2}\\\
-2\sum_{n<m}\epsilon_{|n-m|}b_{n}^{\dagger}b_{n}b_{m}^{\dagger}b_{m}+\sum_{n}Je^{-S(\beta)}\left(b_{n+1}^{\dagger}+b_{n-1}^{\dagger}\right)b_{n},$
(16)

where the temperature dependent coupling constant $S(\beta)$ is the nearest-
neighbor dressing factor given by

$S(\beta)=\frac{1}{N}\sum_{q}\left|\frac{\Delta_{q}}{\Omega_{q}}\right|^{2}\coth\left(\frac{\beta\Omega_{q}}{2}\right)\left(1-\cos(q)\right),$
(17)

where $\beta=1/k_{\text{B}}T$. In the following, the effect of the remaining
coupling $\Delta\hat{H}$ is disregarded and the linear and nonlinear optical
responses of polarons will be computed under the effective Hamiltonian given
by Eq. (16). This approximation is relevant in the anti-adiabatic limit where
the hopping constant $J$ is small. This approach can be improved by treating
the remaining coupling using perturbation theoryPouthier and Falvo (2004);
Pouthier (2013); Yalouz and Pouthier (2016) which can give very reliable
results on a large range of parameters provided that no accidental resonances
occur.Pouthier (2013) However, as a first step this work will only consider
the effective Hamiltonian and the effect of the remaining coupling will be the
subject of a separate study.

Since $\hat{H}_{0}$ commute with the number operator
$\hat{N}=\sum_{n}b_{n}^{\dagger}b_{n}$, $\hat{H}_{0}$ is block diagonal in the
eigenvalues of the operator $\hat{N}$, $v=0,1,2,\dots$. This article focus on
the third-order nonlinear optical response and therefore only the blocks
$v=0,1$ and $v=2$ need to be considered. The one-exciton states block $v=1$ is
trivially diagonalized and the eigenstates are expressed by plane-waves

$\left|k\right\rangle=\frac{1}{\sqrt{N}}\sum_{n}e^{\textrm{i}kn}b^{\dagger}_{n}\left|\varnothing\right\rangle,$
(18)

and the eigenvalues are given by

$\omega_{k}=\tilde{\omega}_{0}+2\tilde{J}(\beta)\cos k,$ (19)

where $\tilde{\omega}_{0}=\omega_{0}-\epsilon_{0}$ is the shifted frequency
and $\tilde{J}(\beta)=Je^{-S(\beta)}$ is the effective hopping constant. The
two-excitons states block $v=2$ can be simplified by using the periodicity of
the lattice. Introducing the following center of mass plane-wave basisPouthier
(2003)

$\left|k\
m\right\rangle=\frac{1}{\sqrt{N}}\sum_{n}e^{\textrm{i}k\left(n+m/2\right)}\xi_{m}b_{n}^{\dagger}b_{n+m}^{\dagger}\left|\varnothing\right\rangle,$
(20)

where $m=0,\dots,(N-1)/2$ is the distance between the two-exciton and where
$\xi_{m}$ is defined by

$\xi_{m}=\begin{cases}0&\mbox{if }m<0,\\\ 1/\sqrt{2}&\mbox{if }m=0,\\\
1&\mbox{if }m>0.\end{cases}$ (21)

Using this basis-set, one can show that the Hamiltonian is block diagonal in
the wave vector $k$. The $k$-th block can be deduced from the
equationsPouthier (2003)

$\displaystyle\hat{H}_{0}\left|k\ 0\right\rangle$
$\displaystyle=\left(2\tilde{\omega}_{0}-2A-2\epsilon_{0}\right)\left|k\
0\right\rangle+\sqrt{2}\tilde{J}_{k}\left|k\ 1\right\rangle,$ (22)
$\displaystyle\hat{H}_{0}\left|k\ 1\right\rangle$
$\displaystyle=\left(2\tilde{\omega}_{0}-2\epsilon_{1}\right)\left|k\
1\right\rangle+\sqrt{2}\tilde{J}_{k}\left|k\
0\right\rangle+\tilde{J}_{k}\left|k\ 2\right\rangle,$ (23)
$\displaystyle\hat{H}_{0}\left|k\ m\right\rangle$
$\displaystyle=\left(2\tilde{\omega}_{0}-2\epsilon_{m}\right)\left|k\
m\right\rangle+\tilde{J}_{k}\left|k\ m-1\right\rangle+\tilde{J}_{k}\left|k\
m+1\right\rangle,$ $\displaystyle\mbox{if }m>1,$ (24)

with $\tilde{J}_{k}=2\tilde{J}\cos(k/2)$. Each block $k$ of $\hat{H}_{0}$ can
then be easily diagonalized numerically giving a set of eigenvalues
$\omega_{k\sigma}$ and eigenvectors $\psi_{k\sigma}(m)$ where
$\sigma=0,\dots,(N-1)/2$ labels the different eigenvalues.

###  II.3 Linear optical response

The coupling of the vibrations to the optical field $E(\textbf{r},t)$ is given
by

$\hat{H}_{\text{int}}=-E(\textbf{r},t)\hat{V},$ (25)

where $\hat{V}$ is the dipole operator expressed for a set of identical
molecules as a function of the projection of the transition dipole moments
$\mu$ on the electric field as

$\hat{V}=\sum_{n}\mu\left(b_{n}+b_{n}^{\dagger}\right).$ (26)

The linear optical response is given by the response function written as

$R^{(1)}(t)=\textrm{i}\Theta(t)\left\langle\left[\hat{V}(t),\hat{V}\right]\right\rangle=\textrm{i}\Theta(t)\left(J(t)-J^{*}(t)\right),$
(27)

where $\Theta(t)$ is the Heaviside function,
$V(t)=e^{\textrm{i}\hat{H}t}\hat{V}e^{-\textrm{i}\hat{H}t}$ is the time
evolution of the dipole operator in the Heisenberg picture and where
$\langle\dots\rangle$ is the thermal average over all degrees of freedom. The
function $J(t)$ can be expressed as a function of the total density matrix at
equilibrium $\hat{\rho}=\exp(-\beta\hat{H})/Z(\beta)$ as

$J(t)=\left\langle\hat{V}(t)\hat{V}\right\rangle=\textrm{Tr}\left[\hat{\rho}e^{\textrm{i}\hat{H}t}\hat{V}e^{-\textrm{i}\hat{H}t}\hat{V}\right].$
(28)

By introducing the Lang-Firsov unitary transformation $\hat{U}$ in the
correlation function, neglecting the remaining coupling $\Delta\hat{H}$,
assuming that the harmonic frequencies of the vibrations are much higher than
the temperature and using the rotating wave approximation, the function $J(t)$
is now written

$J(t)=\sum_{n,m}\mu^{2}\left\langle\varnothing\right|b_{n}e^{-\textrm{i}\hat{H}_{0}t}b^{\dagger}_{m}\left|\varnothing\right\rangle
C_{|n-m|}(t)$ (29)

where $C_{|n-m|}(t)$ is the bath correlation function given by

$C_{|n-m|}(t)=\left\langle\theta_{n}(t)\theta_{m}^{\dagger}\right\rangle_{b}=\exp\left(-g_{|n-m|}(t)\right),$
(30)

where the linebroadening function $g_{n}(t)$ is given by

$g_{n}(t)=\frac{1}{N}\sum_{q}\left|\frac{\Delta_{q}}{\Omega_{q}}\right|^{2}\left\\{\coth\left(\frac{\beta\Omega_{q}}{2}\right)\left(1-\cos\left(\Omega_{q}t-qn\right)\right)+\textrm{i}\sin\left(\Omega_{q}t-qn\right)\right\\}.$
(31)

Note that the linebroadening function for $n=1$ at $t=0$ is simply the
nearest-neighbor dressing factor $g_{1}(0)=S(\beta)$. After straightforward
calculation, the optical response is written

$J(t)=\mu^{2}\sum_{k}e^{-\textrm{i}\omega_{k}t}C_{k}(t),$ (32)

where $C_{k}(t)$ is the spatial Fourier transform of the bath correlation
function

$C_{k}(t)=\sum_{n}e^{\textrm{i}kn}\exp\left(-g_{n}(t)\right).$ (33)

Finally, the absorption spectrum $\alpha(\omega)$ is then directly
proportional to the Fourier transform of the response fonction $R^{(1)}(t)$
given by Eq. (27). For a periodic and isolated system, assuming that the laser
wavelength is larger than the system’s typical size, only the excitation
energy corresponding to a vanishing wavevector $k\rightarrow 0$ should
contribute to the linear optical response. Here, because of the coupling to
the phonon bath, all modes $k$ contribute to the optical response with
different weights corresponding to the spatial Fourier transform of the bath
correlation function. Note that the expression for the linebroadening function
$g_{n}(t)$ is very close to the expression for the linebroadening function of
a single isolated transition coupled to a harmonic bath.Duke and Mahan (1965)
In Eq. (31), the dephasing of the vibrations takes into account the
delocalized nature of the phonons and the correlation induced by the bath. In
addition, the Stokes shift which is usually included in the definition of the
linebroadening functionMukamel (1995) is not present in Eq. (31). This Stokes
shift is in fact included directly in the definition of the polaron
Hamiltonian via the reorganizational energies $\epsilon_{n}$ defined in Eq.
(14).

###  II.4 Third-order nonlinear optical response

The third-order response function is given by

$R^{(3)}(t_{1},t_{2},t_{3})=\textrm{i}^{3}\Theta(t_{1})\Theta(t_{2})\Theta(t_{3})\left\langle\left[\left[\left[\hat{V}(t_{1}+t_{2}+t_{3}),\hat{V}(t_{1}+t_{2})\right],\hat{V}(t_{1})\right],\hat{V}(0)\right]\right\rangle.$
(34)

The three nested commumators yield eight Liouville space pathways.Mukamel
(1995); Abramavicius  _et al._ (2009) Each nonlinear technique is based on a
specific phase matching condition which selects a subgroup of pathways. For
simplification, only the expressions for the signal corresponding to the
direction
$\textbf{k}_{\text{I}}=-\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3}$ are
presented. The expression for the direction
$\textbf{k}_{\text{II}}=\textbf{k}_{1}-\textbf{k}_{2}+\textbf{k}_{3}$ is given
in appendix B. The total response function for $\textbf{k}_{\text{I}}$ is
given as a sum of three contributions

$R_{\textbf{k}_{\text{I}}}(t_{1},t_{2},t_{3})=R_{1}+R_{2}-R_{3},$ (35)

where $R_{1}$ is the ground state bleaching (GSB), $R_{2}$ is the excited
state emission (ESE) and $R_{3}$ is the excited state absorption (ESA). Using
the same approach as for the linear response function each contributions to
the $\textbf{k}_{\text{I}}$ signal can be written as

$\displaystyle R_{1}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N}\sum_{k_{1}k_{2}}e^{\textrm{i}\omega_{k_{1}}t_{1}-\textrm{i}\omega_{k_{2}}t_{3}}C^{(1)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3}),$
(36) $\displaystyle R_{2}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N}\sum_{k_{1}k_{2}}e^{\textrm{i}\omega_{k_{1}}(t_{1}+t_{2})-\textrm{i}\omega_{k_{2}}(t_{3}+t_{2})}C^{(2)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3}),$
(37) $\displaystyle R_{3}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N^{2}}\sum_{k_{1}k_{2}k_{3}\sigma}e^{\textrm{i}\omega_{k_{1}}(t_{1}+t_{2}+t_{3})-\textrm{i}\omega_{k_{2}}t_{2}-\textrm{i}\omega_{k_{3}\sigma}t_{3}}C^{(3)}_{k_{1}k_{2}k_{3}}(t_{1},t_{2},t_{3})A_{k_{1}k_{3}\sigma}A_{k_{2}k_{3}\sigma},$
(38)

with the functions $C^{(i)}(t_{1},t_{2},t_{3})$ defined by

$\displaystyle C^{(1)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{-\textrm{i}k_{1}m_{1}+\textrm{i}k_{2}m_{3}}e^{-g^{(1)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})},$
(39) $\displaystyle C^{(2)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{-\textrm{i}k_{1}(m_{1}+m_{2})+\textrm{i}k_{2}(m_{2}+m_{3})}e^{-g^{(2)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})},$
(40) $\displaystyle C^{(3)}_{k_{1}k_{2}k_{3}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{-\textrm{i}k_{1}(m_{1}+m_{2}+m_{3})+\textrm{i}k_{2}m_{2}+\textrm{i}k_{3}m_{3}}e^{-g^{(3)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})},$
(41)

and where the linebroadening functions are given by

$\displaystyle
g^{(1)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g^{*}_{m_{1}}(t_{1})-g^{*}_{m_{2}}(t_{2})+g_{m_{3}}(t_{3})$
$\displaystyle+g^{*}_{m_{1}+m_{2}}(t_{1}+t_{2})+g^{*}_{m_{2}+m_{3}}(t_{2}+t_{3})-g^{*}_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3}),$
(42) $\displaystyle
g^{(2)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g^{*}_{m_{1}}(t_{1})-g_{m_{2}}(t_{2})+g^{*}_{m_{3}}(t_{3})$
$\displaystyle+g^{*}_{m_{1}+m_{2}}(t_{1}+t_{2})+g_{m_{2}+m_{3}}(t_{2}+t_{3})-g^{*}_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3}),$
(43) $\displaystyle
g^{(3)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g^{*}_{m_{1}}(t_{1})-g_{m_{2}}(t_{2})+g_{m_{3}}(t_{3})$
$\displaystyle+g^{*}_{m_{1}+m_{2}}(t_{1}+t_{2})+g_{m_{2}+m_{3}}(t_{2}+t_{3})-g^{*}_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3}).$
(44)

The tensor $A_{kk^{\prime}\sigma}$ is expressed as a function of the two-
excitons wave function

$A_{kk^{\prime}\sigma}=2\sum_{m}\Psi_{k^{\prime}\sigma}(m)\xi_{m}\cos\left(\left(k^{\prime}/2-k\right)m\right).$
(45)

Note that the discrete spatial Fourier transform in Eqs. (33), (39), (40) and
(41) can be easily computed numerically using the 1D, 2D and 3D Fast Fourier
Transform (FFT) algorithm.Frigo and Johnson (2005)

##  III Results

In this section, the previous formalism is applied to compute the linear and
nonlinear spectroscopy of vibrational polarons in a 1D lattice. The parameters
range used here corresponds to the amide-I vibration in $\alpha$-helix
polypeptides and molecular crystals such as ACN or NMA often modeled as quasi-
one-dimensional chains. The intramolecular anharmonicity value is fixed to
$A=8~{}\textrm{cm}^{-1}$.Hamm, Lim, and Hochstrasser (1998) The hopping
constant ranges between -10 $\textrm{cm}^{-1}$ to 10 $\textrm{cm}^{-1}$. Hamm
and Edler (2006) For the optical phonon model the optical phonon frequency is
fixed to $\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$ corresponding to the
crystalline acetanilide (ACN) optical frequency. Hamm and Edler (2006) For the
acoustical phonon model, the cutoff frequency is fixed to
$\Omega_{\text{ac}}=100~{}\textrm{cm}^{-1}$ corresponding to the
$\alpha$-helices cutoff frequency.Scott (1982) The coupling between the
vibration and the phonons will take a typical value of
$\Delta_{\text{opt}}=\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$. Hamm and Edler
(2006) A phenomenological life-time of $T_{1}=1.5$ ps was added to the
calculation of the linear and non-linear spectra. This value is chosen close
to the relaxation time of 1.2 ps measured for amide-I vibration in
peptides.Hamm, Lim, and Hochstrasser (1998) All numerical calculations were
performed using a number of sites of $N=51$. This number was found to be large
enough to obtain results close to an infinite system. For the case of the
acoustical phonon model, to avoid spurious effects due to the finite size used
in the numerical calculations, the sum over the phonon wavevector $q$ in Eq.
(31) is performed using a larger number of phonon modes $N_{\text{ph}}=5001$.
This number is chosen so that no recursion is observed in the behavior of the
linebroadening function $g_{n}(t)$ given by Eq. (31). Finally the harmonic
frequency is set to the value $\omega_{0}=0$ without any loses of generality.
In the following subsections, the influence of the structure of the bath on
the linear and non-linear vibrational responses is investigated by using
optical and acoustical phonon models.

###  III.1 Optical phonon model

Figure 1: One-polaron (left panel) and two-polarons (right panel) energy
spectra for the optical phonon model as a function of the wave vector $k$ for
$\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$,
$\Delta_{\text{opt}}=25~{}\textrm{cm}^{-1}$, $T=0$ K and for
$J=-10~{}\textrm{cm}^{-1}$ (black solid lines) and $J=0~{}\textrm{cm}^{-1}$
(blue dashed lines).

First, the case of the optical phonon model is considered. The one-polaron
energy spectrum which controls the behavior of the linear absorption spectrum
is depicted on the left panel of Fig. 1. As seen in Eq. (19), the one-polaron
eigenfrequencies are centered around the shifted frequency
$\tilde{\omega}_{0}$ with a width of $4\tilde{J}(\beta)$. For
$\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$,
$\Delta_{\text{opt}}=25~{}\textrm{cm}^{-1}$, $T=0$ K and
$J=-10~{}\textrm{cm}^{-1}$, the shifted frequency takes the value
$\tilde{\omega}_{0}=-12.5~{}\textrm{cm}^{-1}$ and the total dispersion is
$4\tilde{J}=31.2~{}\textrm{cm}^{-1}$. The two-polarons energy spectrum which
controls the behavior of the non-linear optical response is reported on the
right panel of Fig. 1 for the same set of parameters. The energy spectrum
shows one continuum band characterizing the two-polaron free states with a
total bandwidth of $8\tilde{J}=62.4~{}\textrm{cm}^{-1}$ and an isolated band
corresponding to the two-polaron bound states which is determined by the
anharmonicity.Kimball, Fong, and Shen (1981)

Figure 2: Upper panel: Linear absorption spectrum of the optical phonon model
for $\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$,
$\Delta_{\text{opt}}=25~{}\textrm{cm}^{-1}$ and $T=0$ K as a function of the
value of the hopping constant $J$. Lower panel: Linear absorption spectrum for
$\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$,
$\Delta_{\text{opt}}=25~{}\textrm{cm}^{-1}$ and $J=-10~{}\textrm{cm}^{-1}$ as
a function of the temperature $T$.

The upper panel of Fig. 2 shows the dependence of the linear absorption
spectrum for the optical phonon model as a function of the hopping constant
$J$. For $J=-10$ $\textrm{cm}^{-1}$, the spectrum exhibits a sharp and strong
zero-phonon line (ZPL) located at $\omega=-28.1~{}\textrm{cm}^{-1}$
corresponding the the position of the lowest polaron energy
$\omega_{k=0}=\tilde{\omega}_{0}-2\tilde{J}$. A broad second band is also
present in the absorption spectrum. It corresponds to the one-phonon
excitation band (0-1 transition in the Franck-Condon (FC) picture) and is
shifted by $\Omega_{\text{opt}}$ with respect to the ZPL. This broad band
exhibits a double peak shape and has a bandwidth of 31 $\textrm{cm}^{-1}$
corresponding to the total polaron dispersion $4\tilde{J}$. As decreasing the
value of $|J|$ to 0 the ZPL shifts towards $\omega=-12.5~{}\textrm{cm}^{-1}$
while the second band does not shift but its bandwidth reduces significantly.
Upon increasing the value of $J$ to 10 $\textrm{cm}^{-1}$, the ZPL continues
to shift to $\omega=0~{}\textrm{cm}^{-1}$ and the one-phonon band bandwidth
increases again to recover the bandwidth for $J=-10$ $\textrm{cm}^{-1}$. This
behavior of the linear absorption spectrum is almost identical to the
numerical calculation performed by Hamm and Edler based on a direct
diagonalization of the Holstein Hamiltonian. Hamm and Edler (2006) The main
difference between this result and the result of Ref. 23 is seen in the one-
phonon band which is symmetric in this calculation but is stronger on the
lower energy side of the band in Ref. 23. This difference can be explained by
the remaining coupling term $\Delta\hat{H}$ which was neglected here and which
introduce a residual coupling between the ZPL and the one-phonon band.
However, this result captures the main features observed in Ref. 23. The lower
panel of Fig. 2 shows the temperature dependance of the linear absorption
spectrum for the optical phonon model. Upon increasing the temperature, the
linear absorption spectrum exhibits new bands on the blue side of the spectrum
corresponding to $n$th phonon bands as well as hot bands on the red side of
the spectrum. Also by increasing the temperture, the ZPL decreases sharply.
All $n$th phonon bands and hot-bands exhibit the same broad double peak shape
with a bandwidth which corresponds to the polaron dispersion $4\tilde{J}$ and
decreases with temperature. The ZPL shape is also modified as the temperature
is increased.

Figure 3: 2D spectrum for the optical phonon model for
$\Omega_{\text{opt}}=~{}50~{}\textrm{cm}^{-1}$,
$\Delta_{\text{opt}}=25~{}\textrm{cm}^{-1}$ and $T=0$ K and for $J=0$ and
$J=-10$ $\textrm{cm}^{-1}$. Vertical red lines and horizontal red lines mark
the position of the ZPL and the $n$th-phonon bands. The blue vertical lines
correspond to the ZPL and $n$th-phonon lines shifted by the anharmonicity

Fig. 3 shows the absorptive part of the 2D spectrum (sum of the rephasing
$\textbf{k}_{\text{I}}$ and non-rephasing $\textbf{k}_{\text{II}}$ signalsHamm
and Zanni (2011)) for $J=0~{}\textrm{cm}^{-1}$ and $J=-10~{}\textrm{cm}^{-1}$.
For $J=0~{}\textrm{cm}^{-1}$ the 2D spectrum shows multiple negative and
positive peaks. Horizontal red lines have been added to mark the position of
the ZPL and the one-phonon band in the absorption spectrum. The vertical red
lines mark the position of the ZPL, the one-phonon band as well as the one-
phonon hot band (1-0 transition in the FC picture). The one-phonon hot band
originates from the excited state emission. The strongest negative peak
corresponds to the anharmonically shifted ZPL originating from the excited
state absorption. The position of this transition as well as the shifted one-
phonon and two-phonons bands and the shifted hot-band are marked by vertical
blue lines. Note that the shifted hot-bands appear as positive peaks and not
negative peaks even though they originate from the ESA contribution. Analysis
of the response functions shows that the vibrational overlaps corresponding to
these peaks are negative therefore changing the sign of the peaks.

For $J=-10~{}\textrm{cm}^{-1}$ the 2D spectrum is strongly modified as most
bands disappear. Only the ZPL and the one-phonon band remain. As in the linear
absorption spectrum, the one-phonon band is broad with a bandwidth
corresponding to the polaron dispersion. However, this can only being seen
along the $\omega_{1}$ axis of the 2D spectrum. An additional peak appears on
the blue side of the ZPL along the $\omega_{3}$ axis, this peak originates
from the exciton-exciton scattering.Abramavicius (2013) A similar peak is also
visible on the side of the one-phonon band.

###  III.2 Acoustical phonons

Figure 4: One-polaron (left panel) and two-polarons (right panel) energy
spectra for the acoustical phonon model as a function of the wavevector $k$
for $\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ , $T=0$ K and for
$J=-10~{}\textrm{cm}^{-1}$ (black solid lines) and $J=0~{}\textrm{cm}^{-1}$
(blue dashed lines). Figure 5: Upper panel: Linear absorption spectrum of the
acoustical phonon model for $\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ and $T=0$ K as a function of the
value of the phonon coupling $\Delta_{\text{ac}}$ and for $J=0$ (dashed line)
and $J=-10~{}\textrm{cm}^{-1}$ (full line). Lower panel: Linear absorption
spectrum for $\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ and $J=-10~{}\textrm{cm}^{-1}$ as a
function of the temperature $T$. The spectra are normalized to the maximum of
the peak.

Next, the case of the acoustical phonon model is considered. The one-polaron
and two-polarons energy spectra are depicted in Fig. 4 for the parameters
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ , $T=0$ K and for
$J=-10~{}\textrm{cm}^{-1}$ and $J=0~{}\textrm{cm}^{-1}$. The main difference
in the two-polarons energy spectrum with respect to the optical phonon model
can be seen by the appearance of a second bound state. This bound state is
present for all wavevector $k$ for $J=0$ but is only present near $k=\pi/2$
for $J=-10~{}\textrm{cm}^{-1}$. The presence of these two bound states are
related to the occurence of two type of anharmonicities of local and non-local
nature due to the coupling with the phonon bath. Detailed studies on the
nature of these bound states have been performed.Pouthier (2003); Falvo,
Pouthier, and Eilbeck (2006) For example, a phase diagram for the appearance
of the bound states as a function of the anharmonicity, coupling to the bath
and hopping constant, have been drawn using decimation methods.Pouthier (2003)
The presence of two bound states have been suggested to occur for the N$-$H
vibrations in $\alpha$-helix and $\beta$-sheet peptides trough the appearance
of two excited state absorption peaks in their pump-probe spectrum.Edler  _et
al._ (2004, 2005); Bodis  _et al._ (2009)

The upper panel of Fig. 5 shows the dependence of the linear absorption
spectrum for the acoustical phonon model as a function of the coupling
constant $\Delta_{\text{ac}}$ at $T=0$ K and for a hopping constant
$J=-10~{}\textrm{cm}^{-1}$ and $J=0~{}\textrm{cm}^{-1}$. In Fig. 5, the
spectra were normalized with respect to the maximum of the band to highlight
the increase in bandwidth as a function of the coupling. For small values of
the coupling the linear absorption is characterized by a single asymmetric
band. Upon increasing the coupling constant the absorption band bandwidth
increases keeping an asymmetric shape. Only at very large value of the
coupling the shape of the band tends to be more symmetric. The full width at
half maximum (FWHM) of the spectrum increases from $4~{}\textrm{cm}^{-1}$ for
$\Delta_{\text{ac}}=20~{}\textrm{cm}^{-1}$ to $215~{}\textrm{cm}^{-1}$ for
$\Delta_{\text{ac}}=100~{}\textrm{cm}^{-1}$. For small values of
$\Delta_{\text{ac}}$, the effect of the hopping constant $J$ is a simple shift
of the band by $2\tilde{J}$. Upon increasing the value of the coupling this
shift decreases. The lower panel of Fig. 5 shows the temperature dependance of
the linear absorption spectrum for the acoustical phonon model. Upon
increasing the temperature, the bandwidth of the absorption band increases
rapidly and its shape become more symmetric with a typical Lorentzian shape.
The bandwidth (full width at half maximum) increases from
$5~{}\textrm{cm}^{-1}$ for $T=0$ K to $120~{}\textrm{cm}^{-1}$ for $T=200$ K.

Figure 6: 2D spectrum for the acoustical phonon model for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ and $T=0$ K and for $J=0$ and
$J=-10$ $\textrm{cm}^{-1}$.

Next, the nonlinear optical response in the low temperature regime is
discussed. Fig. 6 shows the absorptive part of the 2D spectrum for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$, $T=0$ K and for $J=0$ and $J=-10$
$\textrm{cm}^{-1}$. For $J=0$ the 2D spectrum is characterized by a pair of
negative-positive peaks located on the diagonal of the spectrum originating
from the interference of the ESA and the GSB and ESE pathways and two negative
peaks red shifted along the $\omega_{3}$ axis. These two peaks are the
signature of the two bound states visible in the two-polarons spectrum. For
$J=-10~{}\textrm{cm}^{-1}$, the pair of negative-positive peak is red shifted
by $2\tilde{J}$ and only one negative peak is present. This results from the
fact that only one bound state is present for the wavevector $k=0$.

Figure 7: 2D spectra for the acoustical phonon model for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ and $T=300$ K, $J=-10$
$\textrm{cm}^{-1}$ as a function of the delay time $t_{2}$, for a 1D chain
(left panel) and for a single isolated site $N=1$.

At high temperature, the 2D spectrum is strongly modified. The left panel of
Fig. 7 shows the 2D spectrum for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$, $T=300$ K and for $J=0$ and
$J=-10$ $\textrm{cm}^{-1}$ and for different time delay $t_{2}$. For $t_{2}=0$
the 2D spectrum is characterized by a pair of negative-positive peaks
elongated along the diagonal. The width along the diagonal is much larger than
the width of the low temperature 2D spectrum and corresponds essentially to
the linear absorption bandwidth. The shape of the 2D spectrum is tightly
connected to the separation of homogeneous and inhomogeneous broadening. As
the waiting time $t_{2}$ increases, the shape of the two peaks is strongly
modified toward a more circular shape. This is a well known phenomenon which
has been observed in many molecular systems and is a signature of the spectral
diffusion due to fluctuations of the frequencies.Ishikawa  _et al._ (2007);
Fecko  _et al._ (2005); Falvo  _et al._ (2015) In particular the shape of the
peaks can be directly related to the frequency-frequency correlation function
(FFCF) through metrics computed from 2D lineshape such as the center line
slope (CLS).Kwak  _et al._ (2007); Kwak, Rosenfeld, and Fayer (2008); Falvo
(2016) To understand the effect of the bath-induced correlations, the 2D
spectrum for the same time delays $t_{2}$ and for a single site $N=1$, but
still assuming a infinite phonon model, is computed. This would correspond,
for example, to the spectroscopy of an impurity.Duke and Mahan (1965) For a
single site, coupling to the bath induces a similar elongated shape. But as
$t_{2}$ increases the shape of the peaks is more rapidly circular. This is a
signature of a faster decay of the FFCF.

##  IV Interpretation and discussion

In the previous section, the numerical results have shown that the nature of
the bath can strongly modify the linear and non-linear optical lineshape.
Specific approximations and analytical expressions can be developed to fully
understand the relation between the spectroscopic signature and the processes
involved. In this section I will also discuss the physical meaning and the
experimental implications of these results.

###  IV.1 Optical phonons

For the optical phonon model, the bath coupling and optical frequency are
independent of the wavevector $q$. In this case, very simple expressions can
be derived for the linear absorption. Introducing the Huang-Rhys coupling
constant $S_{\text{opt}}=\Delta_{\text{opt}}^{2}/\Omega_{\text{opt}}^{2}$, the
Stokes shift and the dressing factor are then expressed as

$\displaystyle\epsilon_{n}=\frac{\Delta_{\text{opt}}^{2}}{\Omega_{\text{opt}}}\delta_{n,0}=\Omega_{\text{opt}}S_{\text{opt}}\delta_{n,0},$
(46) $\displaystyle
S(\beta)=S_{\text{opt}}\coth\left(\frac{\beta\Omega_{\text{opt}}}{2}\right).$
(47)

The linebroadening functions are written as

$g_{n}(t)=S_{\text{opt}}\left\\{\coth\left(\frac{\beta\Omega_{\text{opt}}}{2}\right)\left(1-\delta_{n,0}\cos\Omega_{\text{opt}}t\right)+\textrm{i}\delta_{n,0}\sin\Omega_{\text{opt}}t\right\\}.$
(48)

After straightforward calculations, the linear response function can be
written

$J(t)=N\mu^{2}e^{-S(\beta)}\left[e^{-\textrm{i}\omega_{k=0}t}+\tilde{\rho}(t)\sum_{n=-\infty}^{\infty}M_{n}(\beta)e^{-\textrm{i}n\Omega_{\text{opt}}t}\right],$
(49)

where $\tilde{\rho}(t)=\frac{1}{N}\sum_{k}e^{-\textrm{i}\omega_{k}t}$ is the
Fourier transform of the polaron density of states and where the constants
$M_{n}(\beta)$ are defined by

$M_{n}(\beta)=I_{n}\left(S_{0}/\sinh(\beta\Omega_{\text{opt}}/2)\right)e^{n\Omega_{\text{opt}}\beta/2}-\delta_{n,0},$
(50)

where $I_{n}(x)$ are the modified Bessel functions of the first
kind.Abramowitz and Stegun (1972) The linear absorption $\alpha(\omega)$ is
directly proportional to the Fourier transform of the function $J(t)$ and is
written

$\alpha(\omega)=N\mu^{2}e^{-S(\beta)}\left(\delta\left(\omega-\omega_{k=0}\right)+\sum_{n=-\infty}^{\infty}M_{n}(\beta)\rho\left(\omega-n\Omega_{\text{opt}}\right)\right).$
(51)

Therefore the absorption spectrum is given by the sum of a delta-like ZPL and
a series of peaks corresponding to the Franck-Condon vibrational progression.
For $T=0$ K, the ZPL is not broadened by the bath while the shape of the other
bands is given by the polaron density of states
$\rho(\omega)=\frac{1}{N}\sum_{k}\delta(\omega-\omega_{k})$. From Eq. (19), it
is easy to deduce the polaron density of states in the limit
$N\rightarrow\infty$, it is written

$\rho(\omega)=\frac{1}{\pi\sqrt{4\tilde{J}^{2}-\left(\omega-\tilde{\omega}_{0}\right)^{2}}}\quad\text{if}\quad|\omega-\tilde{\omega}_{0}|<2\tilde{J},$
(52)

which diverge when $\omega=\tilde{\omega}_{0}\pm 2\tilde{J}$ resulting in a
double peak shape as observed in Fig. 2. The expressions for the linear
response obtained here are very similar to the expressions given in Ref. 85
which uses a similar derivation. The main difference appears for the ZPL which
in Ref. 85 is broadened by the phonons even at $T=0$ K while in our case it is
not since $M_{0}(\infty)=0$. Note that in Ref. 85, the regime explored
corresponds to a strong Huang-Rhys coupling constant for which the polaron
bandwidth vanishes. Eqs. (51) and (52) fully capture the behavior of the
linear absorption in Fig. 2 and in Ref. 23. Therefore, this work gives a
theoretical basis to understand the nature of the peaks observed in the linear
absorption.

Next, the previous results are used to understand the physical meaning of the
experimental measurements of ACN. In ACN, the experimental linear absorption
spectrum shows at low temperature one peak located at
$1666~{}\textrm{cm}^{-1}$ and a second band at $1650~{}\textrm{cm}^{-1}$. The
band at $1666~{}\textrm{cm}^{-1}$ shows some marked features with three
subbands which disappear at high temperature and the band at
$1650~{}\textrm{cm}^{-1}$ decreases rapidly as a function of temperature (see
for example Refs.17, 18, and 30). This behavior is very similar to the result
presented in the present study using the vibrational Holstein model. As one
can expects from the Franck-Condon picture the intensity of the low frequency
band which corresponds to the ZPL decreases sharply as a function of
temperature while the intensity of the one-phonon band increases. In addition,
the shape of the one-phonon band which corresponds to the polaron density of
states decreases as a function of temperature due to an increase of the
dressing factor (Eq. 52). This is consistent with the disappearance of the
substructure observed in the $1666~{}\textrm{cm}^{-1}$ band of ACN. Note that
in the experiment, the frequency difference between the two peaks is much
smaller than the frequency difference in the 1D Hostein model. This difference
could originate from the fact that we did not include the 3D structure of the
crystal and additional dipole-dipole couplings need to be included.Hamm and
Edler (2006); Hamm (2009) Similarly, the experimental 2D-IR spectrum of ACN
have some marked similarities with the 2D-IR septrum predicted by the
vibrational Holstein model (see for example Refs. 22 and 30). It shows clearly
a pair of negative/positive peaks located at $1650~{}\textrm{cm}^{-1}$ which
would correspond to the ZPL with an aditionnal side peak which we have
interpreted as a result of exciton-exciton scattering. In addition there is no
peaks on the diagonal corresponding to the $1666~{}\textrm{cm}^{-1}$ band
(one-phonon band in the Holstein model).

###  IV.2 Acoustical phonons

For the acoustical phonon model, introducing the coupling constant
$S_{\text{ac}}=2\Delta_{\text{ac}}^{2}/\Omega_{c}^{2}$, the dressing factor
and the Stokes shifts are given by

$\displaystyle S(\beta)=\frac{8S_{\text{ac}}}{\pi}\int_{0}^{1}\textrm{d}x\
x\sqrt{1-x^{2}}\coth\left(\beta\Omega_{c}x/2\right),$ (53)
$\displaystyle\epsilon_{n}=\Omega_{c}S_{\text{ac}}\left(\delta_{n,0}+\frac{1}{2}\delta_{n,1}+\frac{1}{2}\delta_{n,-1}\right),$
(54)

giving rise to two types of couplings: a local anharmonic coupling and a
nearest-neighbor anharmonic coupling. The effect of these two types of
anharmonic couplings on the two-exciton states and the transport properties
has been studies in details.Pouthier (2003); Pouthier and Falvo (2004); Falvo,
Pouthier, and Eilbeck (2006) It is these couplings that are responsible for
the appearance of two bound states in the polaron energy spectrum as seen in
Fig. 4.

####  IV.2.1 Linear absorption

Using Eqs. (32) and (33) one can easily show that in the limit
$N\rightarrow\infty$, the linear optical response is written

$J(t)=N\mu^{2}\sum_{n}e^{-\textrm{i}\tilde{\omega}_{0}t}e^{-g_{n}(t)}(-\textrm{i})^{n}J_{n}\left(2\tilde{J}(\beta)t\right),$
(55)

where $J_{n}(x)$ are the Bessel functions of the first kind. For the case of a
vanishing hopping constant $J=0$, i.e. the anti-adiabatic limit, the optical
response is written

$J(t)=N\mu^{2}\exp\left(-\textrm{i}\tilde{\omega}_{0}t-g_{0}(t)\right),$ (56)

where the linebroadening function $g_{0}(t)$ is written

$g_{0}(t)=\frac{4S_{\text{ac}}}{\pi}\int_{0}^{1}\frac{\textrm{d}x}{x}\sqrt{1-x^{2}}\left\\{\coth\left(\beta\Omega_{c}x/2\right)(1-\cos\Omega_{c}tx)+\textrm{i}\sin\Omega_{c}tx\right\\}.$
(57)

Next, two situations are considered: the case of the strong coupling limit
where $S_{\text{ac}}\gg 1$ and the case of weak coupling limit
$S_{\text{ac}}\ll 1$.

In the strong coupling limit $S_{\text{ac}}\gg 1$, the dynamics is controlled
by the behavior of the correlation function at short times. Therefore, a
second-order Taylor expansion of the linebroadening function can be performed
in the time variable $\tau=\Omega_{c}t$, giving

$g_{0}(\tau)=\textrm{i}S_{\text{ac}}\tau+S(\beta)\tau^{2}/4+\mathcal{O}(\tau^{3}).$
(58)

The absorption spectrum is then Gaussian

$\alpha(\omega)=N\mu^{2}\sqrt{\frac{4\pi}{S(\beta)\Omega_{c}^{2}}}\exp\left(-\frac{(\omega-\omega_{0})^{2}}{S(\beta)\Omega_{c}^{2}}\right).$
(59)

This expression is valid for all temperature range.

In the weak coupling limit $S_{\text{ac}}\ll 1$, the dynamics is controlled by
the behavior of the correlation function at long times which strongly depend
on the temperature range.Duke and Mahan (1965) For the high temperature case
$\beta\Omega_{c}\ll 1$ the linebroadening function is written

$g_{0}(\tau)=\frac{8S_{\text{ac}}}{\pi\beta\Omega_{c}}\int_{0}^{1}\frac{\textrm{d}x}{x^{2}}\sqrt{1-x^{2}}(1-\cos\tau
x),$ (60)

which gives an analytical but cumbersome expression as a function of the
Bessel functions and the Struve functions.Abramowitz and Stegun (1972) At long
timescales one can use the continuum approximation which only consider the
effect of the low-frequency phonons $x=\Omega_{q}/\Omega_{c}\sim
q/2\rightarrow 0$. With this approximation the linebroadening function is
written

$g_{0}(\tau)=\frac{4S_{\text{ac}}}{\beta\Omega_{c}}|\tau|,$ (61)

The linebroadening function therefore increases linearly with time when
$\tau\rightarrow\infty$. The absorption spectrum has a Lorentzian shape with a
width proportional to the temperature and shifted by the bath reorganizational
energy

$\alpha(\omega)=N\mu^{2}\frac{8S_{\text{ac}}/\beta}{\left(\omega-\omega_{0}-\epsilon_{0}\right)^{2}+\left(4S_{\text{ac}}/\beta\right)^{2}}.$
(62)

In the low temperature regime $\beta\Omega_{c}\gg 1$, the linebroadening
function is written as

$g_{0}(\tau)=\frac{4S_{\text{ac}}}{\pi}\int_{0}^{1}\frac{\textrm{d}x}{x}\sqrt{1-x^{2}}\left(1-e^{-\textrm{i}\tau
x}\right).$ (63)

An asymptotic expansion of this expression gives

$g_{0}(\tau)=\frac{4S_{\text{ac}}}{\pi}\left(\gamma-1+\log
2+\log\tau\right)+2\textrm{i}S_{\text{ac}}\operatorname{sgn}\tau,$ (64)

where $\operatorname{sgn}\tau$ is the sign function and where $\gamma\approx
0.577216$ is Euler’s constant. The linebroadening function increases
logarithmically with time. Introducing the coupling constant
$\delta=4S_{\text{ac}}/\pi$, for a weak coupling $\delta<1$ the Fourier
transform of the function $J(t)$ can be calculated and the absorption spectrum
is then given by a power law for
$\omega>\tilde{\omega}_{0}=\omega_{0}-\epsilon_{0}$

$\alpha(\omega)=\frac{2\sin(\delta\pi)\Gamma(1-\delta)e^{-\delta(\gamma-1)}}{(2\Omega_{c})^{\delta}(\omega-\tilde{\omega}_{0})^{1-\delta}}\Theta(\omega-\tilde{\omega}_{0})$
(65)

where $\Gamma(x)$ is the gamma function. Fig. 8 shows a direct comparison of
the absorption spectrum at $T=0$ K computed directly from the linear response
function with the approximation of Eq. (65). The spectrum was centered around
its maximum located at $\tilde{\omega}_{0}=\omega_{0}-\epsilon_{0}$. To
highlight the effect of the phonon broadening, the lifetime $T_{1}$ was
increased to 50 ps. For a coupling constant
$\Delta_{\text{ac}}=50~{}\textrm{cm}^{-1}$, the long time approximation
perfectly match the full calculation for frequencies lower then 50
$\textrm{cm}^{-1}$. As decreasing the coupling constant $\Delta_{\text{ac}}$,
discrepancies occur around the maximum of the absorption band corresponding to
long time behavior from the response function. This originate from the ad-hoc
relaxation $T_{1}$ introduced in the numerical calculation which is not
included in Eq. (65). Note that in the low temperature regime, the effect of
the remaining coupling neglected in this work might control the spectral
bandwidth.

For $J\neq 0$, effects from the polaron dispersion should arise through the
sum over the index $n$ in Eq. (55). However as $n$ increases the
linebroadening function increases quickly and its contribution decreases in
the expression of the response function $J(t)$. For example, in the weak-
coupling limit, at high-temperature and using the continuum approximation the
linebroadening function is written

$g_{n}(\tau)=\frac{2S_{\text{ac}}}{\beta\Omega_{c}}\left(|\tau-2n|+|\tau+2n|\right),$
(66)

which increases linearly with $n$. In this limit, the contribution of the
polaron bandwidth is then negligible. In the low-temperature limit, $g_{n}(0)$
increases with $\log n$ and only small effects of $J$ on the absorption
spectrum should be expected. Fig. 9 shows the linebroadening function
$g_{n}(t)$ as a function of time $t$ and distance $n$ in the low temperature
regime ($T=0$ K) and in the high temperature regime ($T=300$ K). In the low
temperature regime, the real part of the linebroadening function shows clearly
a logarithmic behavior for all $n$ at long timescales. The imaginary part has
a step function behavior. In the high-temperature limit, the real part of the
linebroadening function is constant for $\Omega_{c}t<2m$ and then increases
linearly with time.

Figure 8: Linear absorption spectrum as a function of the phonon coupling
$\Delta_{\text{ac}}$ constant for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$, $J=0~{}\textrm{cm}^{-1}$,
$T=0$ K. To show the specific effect of the phonon broadening the relaxation
time was increased to $T_{1}=50$ ps. The dashed line corresponds to the
absorption spectrum computed from Eq. (65). Figure 9: Linebroadening function
$g_{n}(t)$ computed from the acoustical phonon model for
$\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$ and
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ as a function of the time $t$ and
distance $n$ for $T=0$ K and $T=300$ K.

####  IV.2.2 Nonlinear spectra

To simplify, the case $J=0$ is considered. Since acoustical phonons are
considered, even if the sites are not directly coupled in this limit through
the hopping constant, coupling with the phonon bath can introduce long-range
phonon-mediated correlations. For $J=0$, the sum of GSB and ESE response
functions is given by

$R_{1}+R_{2}=N\mu^{4}e^{\textrm{i}\tilde{\omega}_{0}(t_{1}-t_{3})}\sum_{m}\left(e^{-g^{(1)}_{0,m,0}}+e^{-g^{(2)}_{m,-m,m}}\right).$
(67)

If there is no coupling to the bath, this contribution to the nonlinear signal
scales as $N^{2}$. For a weak coupling with the bath, the scaling of this
contribution will strongly depend on the temperature. In the low temperature
regime, the linebroadening function scales as $g_{m}\approx\delta\log m$. In
this case, the sum of the ESE and GSB contributions scales as $N^{2-\delta}$.
At high-temperature, the linebroadening function scales as $g_{m}\sim m$ and
the sum of the ESE and GSB contributions scales as $N$. The ESA response
function is written as

$\displaystyle R_{3}$
$\displaystyle=2N\mu^{4}e^{\textrm{i}\tilde{\omega}_{0}(t_{1}-t_{3})}\left(e^{-2\textrm{i}(\epsilon_{0}+A)t_{3}}-1\right)e^{-g^{(3)}_{0,0,0}}$
$\displaystyle+2N\mu^{4}e^{\textrm{i}\tilde{\omega}_{0}(t_{1}-t_{3})}\left(e^{-2\textrm{i}\epsilon_{1}t_{3}}-1\right)\left(e^{-g^{(3)}_{1,-1,1}}+e^{-g^{(3)}_{0,1,0}}\right)$
$\displaystyle+N\mu^{4}e^{\textrm{i}\tilde{\omega}_{0}(t_{1}-t_{3})}\sum_{m}\left(e^{-g^{(3)}_{0,m,0}}+e^{-g^{(3)}_{m,-m,m}}\right).$
(68)

The first two terms scale as $N$ while the last term is almost identical to
the sum of the ESE and GSB contributions. In fact this term completely
compensates Eq. (67) if $g^{(1)}_{0,m,0}=g^{(3)}_{0,m,0}$ and
$g^{(2)}_{m,-m,m}=g^{(3)}_{m,-m,m}$. In this case, the spectra is given by one
positive peak located on the diagonal and two negative peaks shifted due to
the two types of anharmonicity. However if $g^{(1)}_{0,m,0}\neq
g^{(3)}_{0,m,0}$ or $g^{(2)}_{m,-m,m}\neq g^{(3)}_{m,-m,m}$ this is not true
anymore. In fact for the ESE contribution the difference can be expressed as

$\displaystyle e^{-g^{(2)}_{m,-m,m}}-e^{-g^{(3)}_{m,-m,m}}$
$\displaystyle\propto e^{-g_{m}^{*}(t_{3})}-e^{-g_{m}(t_{3})}$
$\displaystyle\propto\textrm{i}\sin(g^{\prime\prime}_{m}(t_{3})),$ (69)

where $g^{\prime\prime}_{m}(t)$ is the imaginary part of the linebroadening
function. The imaginary part of the linebroadening function in the long time
approximation is given by

$g^{\prime\prime}_{m}(t_{3})\approx
2S_{\text{ac}}\operatorname{sgn}(\Omega_{c}t_{3}-2m),$ (70)

and do not vanish if $\Omega_{c}t_{3}$ is larger than $2m$. Note that the
difference appears with a $\pi/2$ dephasing which originates from the complex
factor $\textrm{i}=\sqrt{-1}$ in Eq. (69). Consequently, this contribution to
the spectrum is therefore dispersive and results in a pair of positive and
negative peaks located on the diagonal and along the $\omega_{3}$ axis as
noticed in Fig. 6. This effect however will decrease quickly as temperature
increases and the lineshape function decay faster. It can only be seen at very
low temperature as observed in Fig. 6.

On Fig. 7, strong differences have been observed between the case of an
isolated site $N=1$ and the case of a 1D chain coupled to an acoustical phonon
bath. These differences originate mostly from the bath mediated correlations.
To quantify the time evolution of the 2D spectrum, the CLS has been computed.
The CLS is commonly used as a metric to quantify the fluctuation timescales
and extract the FFCF from 2D spectra.Kwak  _et al._ (2007); Kwak, Rosenfeld,
and Fayer (2008); Roy, Pshenichnikov, and Jansen (2011); Falvo (2016) The CLS
for the isolated site and the 1D chain are represented on the left panel of
Fig. 10. For $N=1$ the CLS decays quickly over the first 200 fs with an
oscillation corresponding to the frequency $\Omega_{c}$. This behavior is
characteristic of an underdamped Brownian oscillator. For the 1D chain, the
CLS still exhibits a fast decay over the first 200 fs but then it decays on a
much slower timescale. A convenient way to interpret these results is to
introduce the FFCF for a delocalized system

$D_{n}(t)=\langle\hat{v}_{n}(t)\hat{v}_{0}(0)\rangle,$ (71)

where
$\hat{v}_{n}=\frac{1}{\sqrt{N}}\sum_{q}(\Delta_{q}e^{-\textrm{i}qn}a_{q}^{\dagger}+\Delta_{q}^{*}e^{\textrm{i}qn}a_{q})$
is the frequency fluctuation operator of site $n$. This function measures the
correlation of the fluctuations between two sites separated by the distance
$n$ and after the time $t$. It can be written as

$D_{n}(t)=\frac{1}{N}\sum_{q}|\Delta_{q}|^{2}\left(\coth\left(\frac{\beta\Omega_{q}}{2}\right)\cos(\Omega_{q}t-qn)-\textrm{i}\sin(\Omega_{q}t-qn)\right).$
(72)

For the acoustical phonon model, in the high temperature limit and for
$N\rightarrow\infty$, the bath correlation function is written as a function
of the variable $\tau=\Omega_{c}t$ as

$D_{n}(\tau)=\frac{16\Delta_{\text{ac}}^{2}}{\pi\beta\Omega_{c}}\int_{0}^{1}\textrm{d}x\sqrt{1-x^{2}}\cos(\tau
x)T_{2n}\left(\sqrt{1-x^{2}}\right),$ (73)

where $T_{m}(x)$ are the Chebyshev polynomials of the first kind.Abramowitz
and Stegun (1972) For example for $n=0,1,2$, the first three functions
$D_{n}(\tau)$ are written

$\displaystyle
D_{0}(\tau)=\frac{8\Delta_{\text{ac}}^{2}}{\beta\Omega_{c}}\frac{J_{1}(\tau)}{\tau},$
(74) $\displaystyle
D_{1}(\tau)=\frac{8\Delta_{\text{ac}}^{2}}{\beta\Omega_{c}}\left(\frac{6J_{2}(\tau)}{\tau^{2}}-\frac{J_{1}(\tau)}{\tau}\right),$
(75) $\displaystyle
D_{2}(\tau)=\frac{8\Delta_{\text{ac}}^{2}}{\beta\Omega_{c}}\left(\frac{J_{1}(\tau)(\tau^{2}-120)}{\tau^{3}}-\frac{24J_{2}(\tau)(\tau^{2}-20)}{\tau^{4}}\right).$
(76)

For larger distances $n$ the correlation function can be in principle
calculated analytically but the corresponding expressions are cumbersome and
involve polynomials in $\tau$ of order $n$. The FFCFs are represented on the
right panel of Fig. 10. For $n=0$, the correlation function decays in a
similar manner as the CLS for $N=1$ with the same underdamped oscillation. The
CLS appears just shifted compared to the FFCF. This is not a surprise because
it has been shown that by including fast fluctuation processes the CLS
measures directly the scaled and shifted FFCF.Falvo (2016) As $n$ increases
the maximum of the FFCF is then located at $\Omega_{c}t=n$ and then decreases
quickly with underdamped oscillations. Looking at the maximum value of the
FFCF as a function of $n$ one can see that this maximum decreases slowly in a
similar fashion as the CLS for the 1D chain. This shows that the bath-mediated
long-range correlations measured by the FFCFs is directly responsible for the
slow decay of the CLS in the 2D-IR spectrum.

Figure 10: Left panel: Center line slope of the 2D-IR spectra for the
acoustical model for $\Omega_{\text{ac}}=~{}100~{}\textrm{cm}^{-1}$,
$\Delta_{\text{ac}}=25~{}\textrm{cm}^{-1}$ and $T=300$ K, $J=-10$
$\textrm{cm}^{-1}$ as a function of the delay time $t_{2}$, for a 1D chain
(dashed line) and for a single isolated site $N=1$ (full line). Right panel:
bath correlation functions $D_{n}(t)$ as a function of time.

####  IV.2.3 Experimental implications

In Ref. 28, Edler and coworkers measured the pump-probe spectrum of a model
$\alpha$-helix in the NH spectral range. They observed two bound states that
were interpreted as the signature of the two anharmonic couplings from Eq.
(54): one anharmonic coupling that create a pair of excitons located on the
same site and one anharmonic coupling that creates a pair of excitons located
on two nearest-neighbor sites. A similar observation was then made by Bodis
and coworkers on a model $\beta$-sheet peptide.Bodis  _et al._ (2009) The
model used to interpret the experiment of Ref. 28, relied on a slightly
different Hamiltonian than the one presented in this paper as it included also
the fluctuations of the anharmonicity due to the phonons. To explain the
experimental observations it was assumed a very strong coupling between the NH
vibrations and the accoustical phonon bath, which would correspond with the
present model to a coupling constant of $S_{\text{ac}}=1.2$. However, the
intrepretation used in Ref. 28 did not take into account the spectral
broadening induced by the bath on the linear and the non-linear spectra. In
the high temperature limit, the dressing factor is given by
$S(\beta)=4S_{\text{ac}}k_{\text{B}}T/\Omega_{\text{ac}}$. Using Eq. (59) we
can compute the FWHM of the linear absorption spectrum in the strong coupling
limit. It is given by

$\Delta\omega=4\sqrt{S_{\text{ac}}\Omega_{\text{ac}}k_{\text{B}}T\ln 2}$ (77)

Therefore the FWHM is proportional to the square root of the temperature. For
a temperature of 300 K and a coupling constant $S_{\text{ac}}=1.2$ the FWHM is
$525~{}\textrm{cm}^{-1}$. This value, in complete disagreement with the
measurements, demonstrates that the model is unable to explain all the
experimental observations. In addition, as was observed in Ref. 28, the
spectral bandwidth in the NH range do not change significantly with
temperature. This shows that the Davydov Hamiltonian cannot explain the
emergence of two bound states in the NH pump-probe spectra of $\alpha$-helices
and $\beta$-sheets. Therefore additional theoretical work is needed to explain
these bound states.

##  V Conclusions

In this article, a new methodology to calculate the non-linear response of
vibrational systems based on the small polaron approach is presented. This
approach relies on a unitary transformation which dresses the vibrational
excitation by a phonon cloud and is valid in the anti-adiabatic limit where
the hopping constant is small with respect to the phonon frequency. This
method allows to calculate the optical response of large systems and can
describe explicitly bath-mediated correlation. This method was used to
calculate the linear and non-linear spectroscopy of 1D model chains
considering both optical and acoustical phonon bath.

For the case of an optical phonon bath, a simple expression was derived for
the linear absorption. It shows that the absorption spectrum is characterized
by a main zero-phonon line and a series of $n$th-phonon lines. The $n$th-
phonon lines lineshape is given by the polaron density of states. Here, the
case of a 1D lattice has been considered and the density of states is
characterized by a double peak. This result can be transposed to the more
complicate case of an $n$-dimensional (nD) lattice. The density of states will
strongly depend on the dimensionality of the problem and therefore impact the
shape of the absorption spectrum. This confirms the result of Hamm and
EdlerHamm and Edler (2006) which states that 3D effects might play an
important role to explain linear absorption spectrum of ACN. The approach
presented in this article represents a first step towards a full understanding
of the linear and non-linear spectra of crystaline acetanilide. It will be
extended to include explicitly the 3D structure of ACN allowing therefore a
direct comparison with experiment.

For the case of the acoustical phonon bath, this article shows that in the
C$=$O spectral range, two bound states can be clearly visible in the 2D-IR
spectrum at low temperature. However only one is visible at ambiant
temperature. Up to now, no experimental measurements of 2D-IR spectroscopy of
$\alpha$-helix polypeptides have been made at very low temperature in the
C$=$O spectral range. New measurements in this temperature range could
therefore bring new information on the vibrational dynamics in this system.
This article shows that for the case of acoustical phonons at ambiant
temperature, the spectral diffusion measured from the 2D-IR spectrum appears
much slower for the case of a lattice than for a single site. This can be
explained by bath mediated correlations between distant sites. Experimental
observation of such process could be performed on a model $\alpha$-helix for
example by measuring the 2D-IR spectrum of the $\alpha$-helix and compare it
to an isotope labeled $\alpha$-helix for which the spectrum of a single
isolated vibration coupled to the full phonon bath could be obtained. Finally,
the present model puts in question the validity of the model used to explain
the pump-probe spectra measured in the N$-$H spectral range for a model
$\alpha$-helix and a model $\beta$-sheep peptides. The model predicts a
spectral bandwidth that increases with the square root of the temperature, a
behavior in disagreement with the experimental measurements.

The main limitation of this work resides in the fact that the residual
coupling between the polaron and the bath has been neglected. This coupling
can modify the absorption spectrum and the 2D-IR spectrum. Further theoretical
developments are necessary to fully account for this coupling in the linear
and non-linear response. In particular, inclusion of the remaining coupling by
perturbation theory seems a very promising approach.Pouthier and Falvo (2004);
Pouthier (2013); Yalouz and Pouthier (2016)

###### Acknowledgements.

The author gratefully acknowledges financial support by the Agence Nationale
de la Recherche (ANR) grants ANR-11-BS04-0027 and ANR-16-CE29-0025 as well as
the use of the computing center MésoLUM of the LUMAT research federation (FR
LUMAT 2764).

##  Appendix A Optical and Acoustical phonon models

In the optical phonon model, the bath is characterized by a set of vibrational
coordinates $u_{n}$ and corresponding momentum $p_{n}$ and of harmonic
frequency $\Omega_{\text{opt}}$. The bath Hamiltonian is then written as

$\hat{H}_{b}=\sum_{n}\frac{p_{n}^{2}}{2}+\frac{\Omega_{\text{opt}}^{2}}{2}u_{n}^{2}.$
(78)

I will consider that the exciton frequency of site $n$ is linearly coupled to
the $n$th bath mode. Therefore the system-bath coupling Hamiltonian is written

$\hat{H}_{vb}=\sum_{n}\chi u_{n}b_{n}^{\dagger}b_{n},$ (79)

where $\chi$ is a coupling constant. Using a plane wave basis and the phonon
creation and annihiliation operators the expression for the bath coordinates
$u_{n}$ is given by

$\displaystyle
u_{n}=\frac{1}{\sqrt{2N\Omega_{\text{opt}}}}\sum_{q}\left(e^{\textrm{i}qn}a_{q}+e^{-\textrm{i}qn}a_{q}^{\dagger}\right).$
(80)

Similarly the corresponding momentum $p_{n}$ is written

$p_{n}=\textrm{i}\sqrt{\frac{\Omega_{\text{opt}}}{2N}}\sum_{q}\left(e^{-\textrm{i}qn}a^{\dagger}_{q}-e^{\textrm{i}qn}a_{q}\right).$
(81)

Using these expressions, one can immediately obtain the expression for the
bath and coupling Hamiltonians Eqs. (3) and (4) with
$\Omega_{q}=\Omega_{\text{opt}}$ and
$\Delta_{q}=\chi/\sqrt{2\Omega_{\text{opt}}}$.

For the model of acoustical phonons, the bath Hamiltonian is written as a
function of the bath mode $u_{n}$ as

$\hat{H}_{b}=\sum_{n}\frac{p_{n}^{2}}{2}+\frac{W}{2}\left(u_{n+1}-u_{n}\right)^{2},$
(82)

where $W$ is a coupling constant. Following Davydov,Davydov (1985) the system-
bath coupling Hamiltonian is written as

$\hat{H}_{vb}=\sum_{n}\chi\left(u_{n+1}-u_{n-1}\right)b_{n}^{\dagger}b_{n}.$
(83)

Using the phonon creation and annihilation operators, the expression for the
bath coordinates and momentum are given by

$\displaystyle
u_{n}=\sum_{q}\frac{1}{\sqrt{2N\Omega_{q}}}\left(e^{\textrm{i}qn}a_{q}+e^{-\textrm{i}qn}a_{q}^{\dagger}\right),$
(84) $\displaystyle
p_{n}=\sum_{q}\textrm{i}\sqrt{\frac{\Omega_{q}}{2N}}\left(e^{-\textrm{i}qn}a^{\dagger}_{q}-e^{\textrm{i}qn}a_{q}\right),$
(85)

where $\Omega_{q}=\Omega_{\text{ac}}|\sin q/2|$ and where the cutoff frequency
is given by $\Omega_{\text{ac}}=\sqrt{4W}$. One can then obtain immediately
the expression for the coupling Hamiltonian Eq. (4) with

$\Delta_{q}=-2\textrm{i}\chi W^{-1/4}\sqrt{|\sin q/2|}\cos q/2.$ (86)

##  Appendix B $\textbf{k}_{\text{II}}$ response functions

The signal in the direction $\textbf{k}_{\text{II}}$ can be written as the sum
of three contributions

$R_{\textbf{k}_{\text{II}}}(t_{1},t_{2},t_{3})=R_{4}+R_{5}-R_{6},$ (87)

where $R_{4}$ and $R_{5}$ are respectively the GSB and ESE while $R_{6}$ is
the ESA. These response functions are written

$\displaystyle R_{4}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N}\sum_{k_{1}k_{2}}e^{-\textrm{i}\omega_{k_{1}}t_{1}-\textrm{i}\omega_{k_{2}}t_{3}}C^{(4)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
(88) $\displaystyle R_{5}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N}\sum_{k_{1}k_{2}}e^{-\textrm{i}\omega_{k_{1}}(t_{1}+t_{2}+t_{3})+\textrm{i}\omega_{k_{2}}t_{2}}C^{(5)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
(89) $\displaystyle R_{6}(t_{1},t_{2},t_{3})$
$\displaystyle=\frac{\mu^{4}}{N^{2}}\sum_{k_{1}k_{2}k_{3}\sigma}e^{-\textrm{i}\omega_{k_{1}}(t_{1}+t_{2})+\textrm{i}\omega_{k_{2}}(t_{2}+t_{3})-\textrm{i}\omega_{k_{3}\sigma}t_{3}}C^{(6)}_{k_{1}k_{2}k_{3}}(t_{1},t_{2},t_{3})A_{k_{1}k_{3}\sigma}A_{k_{2}k_{3}\sigma}$
(90)

with the functions $C^{(i)}(t_{1},t_{2},t_{3})$ defined by

$\displaystyle C^{(4)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{\textrm{i}k_{1}m_{1}+\textrm{i}k_{2}m_{3}}e^{-g^{(4)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})}$
(91) $\displaystyle C^{(5)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{\textrm{i}k_{1}(m_{1}+m_{2}+m_{3})-\textrm{i}k_{2}m_{2}}e^{-g^{(5)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})}$
(92) $\displaystyle C^{(6)}_{k_{1}k_{2}}(t_{1},t_{2},t_{3})$
$\displaystyle=\sum_{m_{1}m_{2}m_{3}}e^{\textrm{i}k_{1}(m_{1}+m_{2})-\textrm{i}k_{2}(m_{2}+m_{3})+\textrm{i}k_{3}m_{3}}e^{-g^{(6)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})}$
(93)

and where the linebroadening functions are given by

$\displaystyle
g^{(4)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g_{m_{1}}(t_{1})+g_{m_{2}}(t_{2})+g_{m_{3}}(t_{3})$
$\displaystyle-
g_{m_{1}+m_{2}}(t_{1}+t_{2})-g_{m_{2}+m_{3}}(t_{2}+t_{3})+g_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3})$
(94) $\displaystyle
g^{(5)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g_{m_{1}}(t_{1})+g^{*}_{m_{2}}(t_{2})+g^{*}_{m_{3}}(t_{3})$
$\displaystyle-
g_{m_{1}+m_{2}}(t_{1}+t_{2})-g^{*}_{m_{2}+m_{3}}(t_{2}+t_{3})+g_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3})$
(95) $\displaystyle
g^{(6)}_{m_{1}m_{2}m_{3}}(t_{1},t_{2},t_{3})=g_{m_{1}}(t_{1})+g^{*}_{m_{2}}(t_{2})+g_{m_{3}}(t_{3})$
$\displaystyle-
g_{m_{1}+m_{2}}(t_{1}+t_{2})-g^{*}_{m_{2}+m_{3}}(t_{2}+t_{3})+g_{m_{1}+m_{2}+m_{3}}(t_{1}+t_{2}+t_{3})$
(96)

## References

  * Mahan (1981) G. D. Mahan,  _Many-particle physics_ (Kluwer Academic / Plenum Publishers, New York, 1981).
  * Holstein (1959a) T. Holstein, Ann. Phys. 8, 325 (1959a).
  * Holstein (1959b) T. Holstein, Ann. Phys. 8, 343 (1959b).
  * Spano (2010) F. C. Spano, Acc. Chem. Res. 43, 429 (2010).
  * Huynh  _et al._ (2013) T. D. Huynh, K.-W. Sun, M. Gelin, and Y. Zhao, J. Chem. Phys. 139, 104103 (2013).
  * Lu and Mukamel (1991) N. Lu and S. Mukamel, J. Chem. Phys. 95, 1588 (1991).
  * Sun  _et al._ (2015) K.-W. Sun, M. F. Gelin, V. Y. Chernyak, and Y. Zhao, J. Chem. Phys. 142, 212448 (2015).
  * Chorošajev  _et al._ (2014) V. Chorošajev, A. Gelzinis, L. Valkunas, and D. Abramavicius, J. Chem. Phys. 140, 244108 (2014).
  * Chorosajev, Rancova, and Abramavicius (2016) V. Chorosajev, O. Rancova, and D. Abramavicius, Phys. Chem. Chem. Phys. 18, 7966 (2016).
  * Yamagata and Spano (2014) H. Yamagata and F. C. Spano, J. Phys. Chem. Lett. 5, 622 (2014).
  * Barford and Marcus (2014) W. Barford and M. Marcus, J. Chem. Phys. 141, 164101 (2014).
  * Okamoto  _et al._ (1992) H. Okamoto, T. Mitani, K. Toriumi, and M. Yamashita, Phys. Rev. Lett. 69, 2248 (1992).
  * Fillaux (1981) F. Fillaux, Chem. Phys. 62, 287 (1981).
  * Barthes  _et al._ (1998) M. Barthes, H. N. Bordallo, J. Eckert, O. Maurus, G. de Nunzio, and J. Léon, J. Phys. Chem. B 102, 6177 (1998).
  * Herrebout, Clou, and Desseyn (2001) W. A. Herrebout, K. Clou, and H. O. Desseyn, J. Phys. Chem. A 105, 4865 (2001).
  * Careri  _et al._ (1983) G. Careri, U. Buontempo, F. Carta, E. Gratton, and A. C. Scott, Phys. Rev. Lett. 51, 304 (1983).
  * Careri  _et al._ (1984) G. Careri, U. Buontempo, F. Galluzzi, A. C. Scott, E. Gratton, and E. Shyamsunder, Phys. Rev. B 30, 4689 (1984).
  * Eilbeck, Lomdahl, and Scott (1984) J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott, Phys. Rev. B 30, 4703 (1984).
  * Alexander and Krumhansl (1986) D. M. Alexander and J. A. Krumhansl, Phys. Rev. B 33, 7172 (1986).
  * Edler, Hamm, and Scott (2002) J. Edler, P. Hamm, and A. C. Scott, Phys. Rev. Lett. 88, 067403 (2002).
  * Edler and Hamm (2002) J. Edler and P. Hamm, J. Chem. Phys. 117, 2415 (2002).
  * Edler and Hamm (2003) J. Edler and P. Hamm, J. Chem. Phys. 119, 2709 (2003).
  * Hamm and Edler (2006) P. Hamm and J. Edler, Phys. Rev. B 73, 094302 (2006).
  * Davydov (1973) A. S. Davydov, J. Theor. Biol. 38, 559 (1973).
  * Davydov (1985) A. S. Davydov,  _Solitons in Molecular Systems_ (Springer Science, Dordrecht, 1985).
  * Scott (1982) A. C. Scott, Phys. Rev. A 26, 578 (1982).
  * Scott (1992) A. C. Scott, Phys. Rep. 217, 1 (1992).
  * Edler  _et al._ (2004) J. Edler, R. Pfister, V. Pouthier, C. Falvo, and P. Hamm, Phys. Rev. Lett. 93, 106405 (2004).
  * Edler  _et al._ (2005) J. Edler, V. Pouthier, C. Falvo, R. Pfister, and P. Hamm, in  _Ultrafast Phenomena XIV_ , Springer Series in Chemical Physics, Vol. 79, edited by T. Kobayashi, T. Okada, T. Kobayashi, K. A. Nelson, and S. Silvestri (Springer, 2005) pp. 401–403.
  * Hamm (2009) P. Hamm, J. Biol. Phys. 35, 17 (2009).
  * Brown and Ivić (1989) D. W. Brown and Z. Ivić, Phys. Rev. B 40, 9876 (1989).
  * Ivić  _et al._ (1997) Z. Ivić, D. Kostić, Željko Pržulj, and D. Kapor, J. Phys. Cond. Mat. 9, 413 (1997).
  * Pouthier (2003) V. Pouthier, Phys. Rev. E 68, 021909 (2003).
  * Pouthier and Falvo (2004) V. Pouthier and C. Falvo, Phys. Rev. E 69, 041906 (2004).
  * Falvo and Pouthier (2005a) C. Falvo and V. Pouthier, J. Chem. Phys. 123, 184709 (2005a).
  * Falvo and Pouthier (2005b) C. Falvo and V. Pouthier, J. Chem. Phys. 123, 184710 (2005b).
  * Falvo and Pouthier (2005c) C. Falvo and V. Pouthier, J. Chem. Phys. 122, 014701 (2005c).
  * Tsivlin, Meyer, and May (2006) D. V. Tsivlin, H.-D. Meyer, and V. May, J. Chem. Phys. 124, 134907 (2006).
  * Tsivlin and May (2006) D. V. Tsivlin and V. May, J. Chem. Phys. 125, 224902 (2006).
  * Tsivlin and May (2007) D. V. Tsivlin and V. May, Chem. Phys. 338, 150 (2007).
  * Bodis  _et al._ (2009) P. Bodis, E. Schwartz, M. Koepf, J. J. L. M. Cornelissen, A. E. Rowan, R. J. M. Nolte, and S. Woutersen, J. Chem. Phys. 131, 124503 (2009).
  * Cruzeiro (2009) L. Cruzeiro, J. Biol. Phys. 35, 43 (2009).
  * Goj and Bittner (2011) A. Goj and E. R. Bittner, J. Chem. Phys. 134, 205103 (2011).
  * Miyazawa (1960) T. Miyazawa, J. Chem. Phys. 32, 1647 (1960).
  * Fayer (2013) M. D. Fayer, ed.,  _Ultrafast Infrared Vibrational Spectroscopy_ (CRC Press, 2013).
  * Hamm and Zanni (2011) P. Hamm and M. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, Cambridge, 2011).
  * Khalil, Demirdöven, and Tokmakoff (2003) M. Khalil, N. Demirdöven, and A. Tokmakoff, J. Phys. Chem. A 107, 5258 (2003).
  * Loparo, Roberts, and Tokmakoff (2006) J. J. Loparo, S. T. Roberts, and A. Tokmakoff, J. Chem. Phys. 125, 194522 (2006).
  * Wong  _et al._ (2013) D. B. Wong, C. H. Giammanco, E. E. Fenn, and M. D. Fayer, J. Phys. Chem. B 117, 623 (2013).
  * Bloem  _et al._ (2012) R. Bloem, K. Koziol, S. A. Waldauer, B. Buchli, R. Walser, B. Samatanga, I. Jelesarov, and P. Hamm, J. Phys. Chem. B 116, 13705 (2012).
  * Middleton  _et al._ (2012) C. T. Middleton, P. Marek, P. Cao, C.-C. Chiu, S. Singh, A. M. Woys, J. J. de Pablo, D. P. Raleigh, and M. T. Zanni, Nat. Chem. 4, 355 (2012).
  * Bandaria  _et al._ (2010) J. N. Bandaria, S. Dutta, M. W. Nydegger, W. Rock, A. Kohen, and C. M. Cheatum, Proc. Natl. Acad. Sci. USA 107, 17974 (2010).
  * Ghosh  _et al._ (2014) A. Ghosh, J. Wang, Y. S. Moroz, I. V. Korendovych, M. Zanni, W. F. DeGrado, F. Gai, and R. M. Hochstrasser, J. Chem. Phys. 140, 235105 (2014).
  * Kim and Hochstrasser (2009) Y. S. Kim and R. M. Hochstrasser, J. Phys. Chem. B 113, 8231 (2009).
  * Jansen and Knoester (2009) T. L. C. Jansen and J. Knoester, Acc. Chem. Res. 42, 1405 (2009).
  * Zheng  _et al._ (2005) J. Zheng, K. Kwak, J. Asbury, X. Chen, I. R. Piletic, and M. D. Fayer, Science 309, 1338 (2005).
  * Falvo  _et al._ (2008) C. Falvo, T. Hayashi, W. Zhuang, and S. Mukamel, J. Phys. Chem. B 112, 12479 (2008).
  * Cho (2008) M. Cho, Chem. Rev. 108, 1331 (2008).
  * Tempelaar  _et al._ (2013) R. Tempelaar, A. Stradomska, J. Knoester, and F. C. Spano, J. Phys. Chem. B 117, 457 (2013).
  * Lang and Firsov (1963) I. G. Lang and Y. A. Firsov, Sov. Phys. JETP 16, 1301 (1963).
  * Yalouz, Falvo, and Pouthier (2017) S. Yalouz, C. Falvo, and V. Pouthier, Quantum Inf. Process. 16, 143 (2017).
  * Jeckelmann and White (1998) E. Jeckelmann and S. R. White, Phys. Rev. B 57, 6376 (1998).
  * Chen, Zhao, and Tanimura (2015) L. Chen, Y. Zhao, and Y. Tanimura, J. Phys. Chem. Lett. 6, 3110 (2015).
  * Philpott (1971) M. R. Philpott, J. Chem. Phys. 55, 2039 (1971).
  * Spano (2002) F. C. Spano, J. Chem. Phys. 116, 5877 (2002).
  * Yalouz, Pouthier, and Falvo (2017) S. Yalouz, V. Pouthier, and C. Falvo, Phys. Rev. E 96, 022304 (2017).
  * Woutersen (2007) S. Woutersen, J. Chem. Phys. 126, 226101 (2007).
  * Pouthier (2013) V. Pouthier, J. Chem. Phys. 138, 044108 (2013).
  * Yalouz and Pouthier (2016) S. Yalouz and V. Pouthier, Phys. Rev. E 93, 052306 (2016).
  * Duke and Mahan (1965) C. Duke and G. Mahan, Phys. Rev. 139, A1965 (1965).
  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ (Oxford University Press, 1995).
  * Abramavicius  _et al._ (2009) D. Abramavicius, B. Palmieri, D. V. Voronine, F. Šanda, and S. Mukamel, Chem. Rev. 109, 2350 (2009).
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, Proc. IEEE 93, 216 (2005).
  * Hamm, Lim, and Hochstrasser (1998) P. Hamm, M. Lim, and R. M. Hochstrasser, J. Phys. Chem. B 102, 6123 (1998).
  * Kimball, Fong, and Shen (1981) J. C. Kimball, C. Y. Fong, and Y. R. Shen, Phys. Rev. B 23 (1981).
  * Abramavicius (2013) D. Abramavicius, Europhys. Lett. 101, 57007 (2013).
  * Falvo, Pouthier, and Eilbeck (2006) C. Falvo, V. Pouthier, and J. C. Eilbeck, Physica D 221, 58 (2006).
  * Ishikawa  _et al._ (2007) H. Ishikawa, I. J. Finkelstein, S. Kim, K. Kwak, J. K. Chung, K. Wakasugi, A. M. Massari, and M. D. Fayer, Proc. Natl. Acad. Sci. USA 104, 16116 (2007).
  * Fecko  _et al._ (2005) C. J. Fecko, J. J. Loparo, S. T. Roberts, and A. Tokmakoff, J. Chem. Phys. 122, 054506 (2005).
  * Falvo  _et al._ (2015) C. Falvo, L. Daniault, T. Vieille, V. Kemlin, J.-C. Lambry, C. Meier, M. H. Vos, A. Bonvalet, and M. Joffre, J. Phys. Chem. Lett. 6, 2216 (2015).
  * Kwak  _et al._ (2007) K. Kwak, S. Park, I. J. Finkelstein, and M. D. Fayer, J. Chem. Phys. 127, 124503 (2007).
  * Kwak, Rosenfeld, and Fayer (2008) K. Kwak, D. E. Rosenfeld, and M. D. Fayer, J. Chem. Phys. 128, 204505 (2008).
  * Falvo (2016) C. Falvo, J. Chem. Phys. 144, 234103 (2016).
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun, eds.,  _Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables_ (Dover, New York, 1972).
  * Čevizović, Zeković, and Ivić (2009) D. Čevizović, S. Zeković, and Z. Ivić, Chem. Phys. Lett. 480, 75 (2009).
  * Roy, Pshenichnikov, and Jansen (2011) S. Roy, M. S. Pshenichnikov, and T. L. C. Jansen, J. Phys. Chem. B 115, 5431 (2011).