AA YYYY

# The Optical Signatures of Stochastic Processes in Many-Body Exciton
Scattering

Hao Li  1 S. A. Shah  1 Ajay Ram Srimath Kandada2 Carlos Silva  3,4,5 Andrei
Piryatinski  6 Eric R. Bittner1 1 Department of Chemistry, University of
Houston, Houston, Texas 77204, United States  2 Department of Physics and
Center for Functional Materials, Wake Forest University, 1834 Wake Forest
Road, Winston-Salem, North Carolina 27109, United States  3 School of
Chemistry and Biochemistry, Georgia Institute of Technology, 901 Atlantic
Drive, Atlanta, GA 30332, United States  4 School of Physics, Georgia
Institute of Technology, 837 State Street, Atlanta, GA 30332, United States  5
School of Materials Science and Engineering, Georgia Institute of Technology,
North Avenue, Atlanta, GA 30332, United States  6 Theoretical Division, Los
Alamos National Laboratory, Los Alamos, NM, 87545 United States

###### Abstract

We review our recent quantum stochastic model for spectroscopic lineshapes in
the presence of a co-evolving and non-stationary background population of
excitations. Starting from a field theory description for interacting bosonic
excitons, we derive a reduced model whereby optical excitons are coupled to an
incoherent background via scattering as mediated by their screened Coulomb
coupling. The Heisenberg equations of motion for the optical excitons are then
driven by an auxiliary stochastic population variable, which we take to be the
solution of an Ornstein–Uhlenbeck process. Here we discuss an overview of the
theoretical techniques we have developed as applied to predicting coherent
non-linear spectroscopic signals. We show how direct (Coulomb) and exchange
coupling to the bath give rise to distinct spectral signatures and discuss
mathematical limits on inverting spectral signatures to extract the background
density of states.

###### doi:

10.1146/((please add article doi))

###### keywords:

excitation-induced dephasing, many-body effects in quantum dynamics, coherent
non-linear spectroscopy, quantum stochastic calculus

††journal: Xxxx. Xxx. Xxx. Xxx.

###### Contents

  1. 1 INTRODUCTION
  2. 2 NONLINEAR SPECTROSCOPIC SIGNATURES OF EID
  3. 3 NONLINEAR COHERENT SPECTROSCOPY OF NONSTATIONARY SYSTEMS
    1. 3.1 Optical Bloch Equations
    2. 3.2 Stochastic Many-body Processes
    3. 3.3 Exciton scattering contributions in the mean-field (Hartree) limit and quantum Langevin equations.
    4. 3.4 Predictions from the stochastic model
      1. 3.4.1 Linear response
      2. 3.4.2 Two-dimensional coherent spectroscopy
    5. 3.5 Exciton/polaron formation dynamics due to exchange interactions
      1. 3.5.1 Effect on 2D spectroscopy
      2. 3.5.2 Comparison to the Anderson-Kubo model and our previous excitation-induced dephasing (EID) theory
  4. 4 Perspective

##  1 INTRODUCTION

It is well recognized that many-body phenomena have a profound effect on the
linear and non-linear optical lineshapes of semiconductors with reduced
dimensionality, in which Coulomb correlations can be particularly strong due
to decreased screening and quantum confinement effects. One such effect is
biexciton formation, in which Coulomb binding of two electron-hole pairs
results in new two-electron, two-hole quasiparticles [1, 2, 3, 4, 5, 6, 7, 8,
9, 10, 11]. Another important process that is highly relevant in exciton
quantum dynamics is excitation induced dephasing (EID) [12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], primarily investigated in
two-dimensional (2D) systems such as III-V quantum wells [13, 18, 19, 20, 22,
23, 24], single-layer transition-metal dichalchogenides [25, 26], quantum dot
photocells [28], and two-dimensional metal-halide perovskite derivatives [27].
This can be described as the incoherent Coulomb elastic scattering between
multiple excitons or between excitons and an electron-hole plasma generated
with the excitation optical field. The scattering process gives rise to faster
dephasing dynamics compared to the low-density pure-dephasing limit, and may
be the dominant dephasing pathway at sufficiently high densities. In many
systems, especially those with strong exciton-phonon coupling, the background
excitations are transient and co-evolve with optical modes of the system and
consequently a strictly incoherent kinetic description such as this mesoscopic
approach or a kinetic Markovian Boltzmann-like scattering theory [15] cannot
describe coherence dynamics. [] EIDExcitation
induced dephasing EISExcitation
induced shift EID can be effectively rationalized from a mesoscopic
perspective by means of the optical Bloch equations, which capture the effect
of many-body exciton scattering on both population and coherence dynamics
derived from coherent spectroscopy of semiconductors [20, 21].

Recent advances towards a more microscopic perspective has been presented by
Katsch et al., in which excitonic Heisenberg equations of motion are used to
describe linear excitation line broadening in two-dimensional transition-metal
dichalchogenides [29]. Their results indicate exciton-exciton scattering from
a dark background as a dominant mechanism in the power-dependent broadening
EID and sideband formation. Similar theoretical modelling on this class of
materials and their van der Waals bilayers have yielded insight into the role
of effective mass asymmetry on EID processes [30]. These modelling works
highlight the need for microscopic approaches to understand nonlinear quantum
dynamics of complex 2D semiconductors, but the computational expense could
become considerable if other many-body details such as polaronic effects are
to be included [31]. As an alternative general approach, we recently developed
an analytical theory of dephasing in the same vein as Anderson-Kubo lineshape
theory [32, 33], but that includes transient EID and Coulomb screening
effects, would be valuable to extract microscopic detail on screened exciton-
exciton scattering from time-dependent nonlinear coherent ultrafast
spectroscopy, via direct and unambiguous measurement of the homogeneous
excitation linewidth [34, 35].

Here we present an overview of our work that employs a quantum stochastic
approach, derived from a first-principles many-body theory of interacting
excitons, to develop a mostly analytical model that describes linear and
nonlinear spectral lineshapes that result from exciton-exciton scattering
processes, and, importantly, their dependence on population time due to the
evolution of a non-stationary/non-equilibrium excitation background (see Fig.
1(a)). Our approach is similar in spirit to the celebrated Anderson-Kubo
theory [32, 33] and reduces to that in the limit of a stationary background
population at sufficiently long times [36, 36]. The model captures a
microscopic picture of EID by integrating over the interactions of excitons
produced via a well-defined coherent pathway (Fig. 10 below).The background
excitons that do not have a well-defined phase relationship induced by the
optical field and can be treated as a non-stationary source of quantum noise.
In doing so, we can directly insert the spectral density of the bath into non-
linear spectral response functions and obtain fully analytical expressions for
the coherent exciton lineshapes.

We implement the model to investigate the evolution of the two-dimensional
coherent excitation lineshape in a polycrystalline thin film of a prototypical
two-dimensional single-layer metal-halide perovskite derivative,
phenylethylammonium lead iodide [(PEA)2PbI4] (see Fig. 1(c) for the crystal
structure). We have selected this material as a model system because of its
well-defined exciton lineshape that we have modeled quantitatively within a
Wannier-Mott framework [37] and because it displays strong many-body phenomena
— strongly bound biexcitons at room temperature [38], and robust EID effects
[27]. Furthermore, we have concluded that the primary excitations are exciton
polarons [39, 31] — quasiparticles with Coulomb correlations that are
renormalized by lattice dynamics via polaronic effects; both electron-hole and
photocarrier-lattice correlations are ingredients of the system Hamiltonian
such that the lattice dressing constitutes an integral component of its
eigenstates and eigenvalues. This renders the system in a highly dynamically
disordered state such that lattice screening effects play an important role in
shaping the linewidth [27] and in dictating nonadiabatic dynamics [40]. We
measure the dephasing dynamics via the homogeneous linewidth extracted by
means of two-dimensional coherent excitation spectroscopy [34, 35]. In our
measurements, excitons generated coherently by a sequence of time-ordered and
phase-matched femtosecond pulses scatter from incoherent background excitons
and thereby undergo EID, which is perceived via changes of the homogeneous
linewidth. We find that EID affects the complex lineshape by mixing absorptive
and dispersive features in the real and imaginary spectral components; the
real component of the two-dimensional coherent spectrum initially displays a
dispersive lineshape that evolves into an absorptive over the timescale in
which EID couplings persist, and the imaginary component evolves in the
converse fashion. Furthermore, we find that the homogeneous contribution to
the spectral linewidth narrows with population time, indicating a dynamic
slowing down of the dephasing rate as the EID correlations active at early
time dissipate. We find that the dynamic line narrowing phenomenon is
reproduced by our stochastic scattering theory, which allows us to explore the
effect of dynamic Coulomb screening on EID quantum dynamics.

Figure 1: (a) Schematic representation of optical absorption of excitons and
exciton-exciton scattering with a background population, where the dispersion
relation is in the exciton representation and
$\vec{k}=\vec{k_{e}}+\vec{k_{h}}$ is the exciton wavevector. (b) Time
evolution of population $N(t)/N_{\mathrm{eq}}$ from an initial nonstationary
state produced by exciton injection. Individual trajectories are represented
by blue dots. Asymptotically, the function reaches a stationary state that
yields the Anderson-Kubo limit. (c) Crystal structure of (PEA)2PbI4 with
schematic representation of exciton-exciton elastic scattering interactions.
Reproduced with permission from ref. 41. Copyright 2020 American Institute of
Physics.

##  2 NONLINEAR SPECTROSCOPIC SIGNATURES OF EID

Figure 2: (a) Linear absorption spectrum of (PEA)2PbI4 at 5 K. Real part of
the corresponding rephasing (b) and non-rephasing (c) spectra at a population
time of $\tau_{p}=0$ fs and at 5 K. The bar to the right of the figure
displays the vertical false color scale in arbitrary units. Reproduced with
permission from ref. 41. Copyright 2020 American Institute of Physics.

Two-dimensional coherent electronic spectroscopies are very powerful
techniques to identify and quantify many-body effects in semiconductors, and
have been instrumental in the study EID of excitons in 2D materials [21, 22,
23, 24, 27, 41, 42]. Much of the early work was carried out by the group of
Cundiff on semiconductor nanostructures [21, 22, 23, 24], but here we focus on
our recent work a Ruddlesden-Popper metal halide, namely on (PEA)2PbI4 [27,
41, 42], (PEA = phenylethylamine), which is a 2D analogue of a lead-iodide
perovskite structure (see Fig. 1(c)). We highlight the peculiar signatures of
EID on the 2D exciton complex lineshape, the effect of exciton density on the
homogeneous linewidth, and finally on the time evolution of the spectral
lineshape.

We start this review by considering the complex 2D coherent excitation
spectrum of (PEA)2PbI4 to quantify the consequences of EID in the nonlinear
lineshape. The linear absorption spectrum (Fig. 2(a)) reveals a family of
exciton polarons with binding-energy offsets of $\sim 40$ meV [37, 39, 31]; we
label the dominant excitons as $X_{A}$ and $X_{B}$, and a shoulder at the blue
edge of $X_{A}$ as $X_{A^{\prime}}$. Shown in Figs. 2(b) and 2(c) are the real
parts of two different coherent excitation pathways; the time-ordering of the
three optical pulses in the experiment and phase-matching conditions define
the specific excitation pathways, based on which rephasing [Fig. 2(b)] and
non-rephasing [Fig. 2(c)] spectra are obtained [43]. In the rephasing
experiment, the pulse sequence is such that the phase evolution of the
polarization after the first pulse and the third pulse are of opposite sign,
while in the non-rephasing experiment, they are of the same sign (see equation
52 and Fig. 10). Both measurements shown in Fig. 2 are taken at a population
waiting time $\tau_{p}=0$ fs and an excitation fluence of 40 nJ/cm-2, which
corresponds to an exciton density in which we have identified effects of
elastic exciton-exciton scattering [27]. [] Rephasing vs. NonrephasingPhoton
momentum conservation requires the output signal wave vector, $k_{s}$ be the
sum of the input wave vectors. These can be understood by inspection of
Fig.10. (c.f.[44]) Rephasing
$\vec{k}_{s}=-\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3}$ Non-
rephasing$\vec{k}_{s}=\vec{k}_{1}-\vec{k}_{2}+\vec{k}_{3}$ $\tau_{p}$Time
between pulses 2 and 3 in which system is in a density-matrix population
state: $|n\rangle\langle n|$. Corresponding diagonal spectral features at the
energies of $X_{A}$, $X_{A^{\prime}}$ and $X_{B}$ (indicated by the magenta
vertical dotted-lines in Fig. 2) are observed, both in rephasing and non-
rephasing spectra. Apart from these diagonal peaks, we observe an off-diagonal
excited-state absorption feature (opposite phase with respect to the diagonal
features) corresponding to a correlation between the absorption energy of
$X_{A}$ and emission energy $\sim 2.3$ eV, which has no corresponding diagonal
signal. We have assigned this cross-peak to a biexciton resonance [38].

Figure 3: Fluence dependence of the exciton dephasing rates. Dephasing
parameters $\gamma$ of $X_{A}$, $X_{A^{\prime}}$ and $X_{B}$ are obtained from
the simultaneous fitting of diagonal and anti-diagonal cuts of the norm of the
zero-time rephasing spectrum, plotted as a function of excitation fluence.
Squares represent the experimental linewidths and lines are the best fit to
the model described in the main text. Error bars on the data are contained
within the markers. The sample temperature was maintained at 5 K. Reproduced
with permission from ref. 27. Copyright 2019 American Physical Society.

From the norm of the rephasing spectrum at zero time [not shown in Fig. 2 but
shown below in Fig. 4(i)], one can extract the homeogeneous and inhomogeneous
linewidths via a global analysis of the diagonal and the anti-diagonal
lineshape [34, 35]. To assess the contribution of many-body interactions on
the dephasing dynamics of the different excitons, we acquired 2D coherent
excitation spectra for a wide range of excitation fluences and sample
temperatures (the raw data is presented in Ref. 27). The monotonic rise of
$\gamma$ with the excitation fluence at 5 K is shown in Fig. 3. Such a
dependence on exciton density $n$ is a consequence of broadening induced by
exciton-exciton elastic scattering mediated by long-range Coulomb
interactions:

$\gamma_{\mathrm{EID}}(n)=\gamma_{0}+\Delta\cdot n.$ (1)

Here, $\gamma_{0}$ is the density-independent dephasing rate and $\Delta$ is
the exciton-exciton interaction parameter [25, 34, 35]. Excitons in
Ruddlesden-Popper metal halides are confined to one of the inorganic quantum
wells and are electronically isolated from the others due to the large inter-
layer distance [38] ($\sim 8$ Å) imposed between them by the long organic
cations. However, the sample itself, 40-nm thick, is composed of tens of these
quantum wells, leading to a highly anisotropic exciton-exciton interaction. To
quantify EID effects, we report the exciton-exciton interaction parameter,
$\Delta$ in units of energy per area. The associated fits and the fit
parameters are displayed in Fig. 3(a), (b) and (c). While $\gamma_{0}$ is
approximately 2 meV with modest variation across the three probed excitonic
transitions, $\Delta$ varies more substantially. It is $2.7\times 10^{-12}$
$\mu$eV cm2 for $X_{B}$, increases to $4.3\times 10^{-12}$ $\mu$eV cm2 for
$X_{A}$, and to $8.6\times 10^{-12}$ $\mu$eV cm2 for $X_{A{{}^{\prime}}}$.
Furthermore, for all excitons, $\gamma$ is consistently smaller with linearly
polarized excitation than with circularly polarized pulses [27]. Moreover,
except for exciton $X_{A^{\prime}}$, it is independent on the helicity of the
exciting pulses within our experimental uncertainty [27]. When exciting the
sample with linearly polarized pulses, excitons can scatter on both left and
right circularly polarized excitons. However, they can only scatter with
excitons of the same polarization when excited with circularly polarized
pulses due to conservation of angular momentum [45].

It is not straight-forward to compare these values with those of other
materials due to the ambiguities over the relevant values of the permittivity
function and thus the Bohr radii. However, given the two-dimensional nature of
the exciton and comparable exciton binding energies, monolayers of transition
metal-dichalchogenides (TMDCs) provide a realistic benchmark. Intriguingly,
previous 2D coherent excitation measurements on unencapsulated WSe2 [25] and
encapsulated MoSe2 [26] revealed $\Delta=2.7\times 10^{-12}$ meV cm2 and
$4\times 10^{-13}$ meV cm2, respectively, three and two orders of magnitude
higher than the value obtained here for (PEA)2PbI4. This and the linearity of
the dephasing rates over a wide range of excitation densities [45] highlights
the substantial screening of the exciton-exciton interactions in these
Ruddlesden-Popper metal halides. This is especially surprising given the high
biexciton binding energy [38], another characteristic that (PEA)2PbI4 shares
with TMDC monolayers [46, 47].

Figure 4: (a)–(d): Real parts of experimentally measured rephasing spectra at
population times $\tau_{p}$ indicated at the top of each panel, measured at 5
K. (e)–(h): Corresponding imaginary parts of the spectrum. (i)–(l): The norm
(absolute value) of the optical response. All spectra components are plotted
in the same relative vertical color scale to facilitate comparison of the
time-dependent signal. Reproduced with permission from ref. 41. Copyright 2020
American Institute of Physics.

[] TDMCtransition metal dichalchogenides Absorptive lineshape a diagonal slice
of a 2D spectrum directly maps onto the linear absorption spectra of the
system and has even reflection symmetry along the off-diagonal axis as shown
in Fig. 7 for a system of J-aggregates which do not exhibit the EID effect.
Dispersive lineshape the 2D spectrum has odd reflection symmetry along the
off-diagonal axis. as can be seen in Fig. 6 for GaAs multiple quantum wells
which show strong EID and EIS effects.

We now return to the complex zero-time spectral lineshape displayed in Fig. 2.
Upon close inspection, we notice that the real part of the spectrum displays
dispersive shape, i.e. derivative shape about the peak energy, both for
diagonal and off-diagonal resonances, in both the rephasing and non-rephasing
spectrum. Note the sign-flip for the off-diagonal feature, which is consistent
with its assignment to the excited state absorption to the biexcitonic state
[38]. Similarly, the imaginary part of the spectra (not shown in Fig. 2 but
shown in Fig. 4) display an absorptive lineshape. We have demonstrated that
such dispersive lineshapes are a consequence of many-body correlations [41],
consistent with the analysis of similar measurements in semiconductor quantum
wells [22]. These lineshapes are unexpected in the absence of many-body
correlations; the real part of the spectrum should be absorptive while the
imaginary part dispersive. The spectra in Fig. 2 therefore reveal phase mixing
due to many-body Coulomb correlations responsible for EID. In fact, Fig. 3
indicates that the EID dominates the non-linear response in the employed pump
fluence range.

The evolution of the rephasing lineshape shown in Fig. 2(b) with population
waiting time $\tau_{p}$ is displayed in Fig. 4. The top row displays the real
part of the spectrum at different values of $\tau_{p}$, the middle row the
imaginary component, and the bottom row the norm (absolute value) of the
complex spectrum. We observe that the phase scrambling phenomenon displayed in
the $\tau_{p}=0$ fs spectrum [Fig. 2(b)] dissipates within $\tau_{p}\leq 240$
fs: the real component of the spectrum evolves from an initially dispersive
[Fig. 4(a)] to absorptive [Fig. 4(d)] lineshape, while that of the imaginary
part evolves from absorptive [Fig. 4(e)] to dispersive [Fig. 4(h)] character.
We note that although the evolution of the real and imaginary components of
the complex lineshape is substantial over this ultrafast time window, the
population decay of the diagonal features for $X_{A}$ and $X_{B}$ is weak,
observed via the modest evolution of the total intensity in Fig. 4(i)–(l). The
decay of the $X_{A^{\prime}}$ diagonal peak and the biexciton cross peak
appears more substantial.

We also highlight the reduction in the total linewidth of the each diagonal
exciton resonance in the absolute value of the response shown in Fig. 4(i)–(l)
with population time. Inspection of these spectra reveal dynamic narrowing of
$X_{A}$ and $X_{B}$, primarily along the anti-diagonal spectral axis. It is
more difficult to visually ascertain the linewidth evolution of
$X_{A^{\prime}}$ and the biexciton cross peak given the non-negligible decay
over this time period.

Figure 5: Homogenous linewidths obtained from the lineshape analysis of the
absolute value of the rephasing spectra (see reference 27) plotted as a
function of the population time for (a) $X_{A}$, (b) $X_{A^{\prime}}$ and (c)
$X_{B}$ exciton lines shown in Fig. 2(a). Reproduced with permission from ref.
41. Copyright 2020 American Institute of Physics.

To quantify the measured dynamic line narrowing, we display in Fig. 5 the
homogeneous linewidth as extracted as in reference 27 as a function of
population time $\tau_{p}$. By this analysis, Fig. 5 shows that the linewidth
of $X_{A}$ reduces most drastically, but that of $X_{B}$ also reduces over a
typical time window, while $X_{A^{\prime}}$ displays no line narrowing. We
note that in reference 27 and in Fig. 3, we reported that $X_{A}$ has a
stronger density dependence of EID than $X_{B}$, which is consistent with the
observation derived from Figs. 5(a) and 5(c). We have found $X_{B}$ to be more
strongly displaced along phonon coordinates involving octahedral twist in the
plane of the inorganic layer, and out of plane scissoring of the Pb—I—Pb apex
[39]. The stronger exciton-phonon coupling implies that $X_{B}$ is more
susceptible to dynamic screening than $X_{A}$, which is consistent with the
data in Fig. 5 and Fig. 3. Finally, we point out that the asymptotic value of
the homogeneous linewidth for $X_{A}$, $X_{A^{\prime}}$, and $X_{B}$ tends
towards the low-exciton-density linewidths that we reported in reference 27.

The linewidth of $X_{A^{\prime}}$ remains relatively constant over the probed
population time. While this might initially suggest that this resonance is
immune to EID effects, we note that the real part of the rephasing spectrum
associated to this particular transition exhibits a dispersive lineshape at
all population times, consistent with the inital lineshapes of $X_{A}$ and
$X_{B}$. This indicates the clear presence of EID effects, as also confirmed
by the density dependent linewidth previously published in Ref 27. The trend
shown in Fig 5(b), on the other hand, suggests that the inter-exciton
scattering does not evolve with the population time, at least within the
probed time range. Inspection of the the lineshape, however, suggests that the
dispersive shape of the real part is preserved at all population times,
suggesting $X_{A^{\prime}}$ is subjected to EID over a much longer period of
time than the other two resonances. Following the arguments developed by the
theoretical work described below in the review, this implies the presence of a
background exciton population that contributes to the scattering of
$X_{A^{\prime}}$ and whose stochastic evolution is that of the background of
the other two resonances. This reiterates our assignment of the multiple
resonances within the spectral structure to excitonic states of distinct
character and possibly specific origin [31].

Figure 6: Linear absorption, excitation pulse spectrum (a) and experimental
real spectra for the rephasing (b) and nonrephasing (c) pulse sequences.
Reproduced with permission from Ref. 22.

The early-time complex lineshape indicative of EID many-body correlations,
shown in Figs. 2 and 4, have been observed previously in GaAs quantum wells
[22]. Fig. 6 displays the real part of the zero-population-time rephasing and
non-rephasing spectra at excitation densities in which the signature of many-
body interaction is clearly observed, and we identify a dispersive lineshape
akin to that observed in Fig. 2 for (PEA)2PbI4. In contrast, Fig. 7 displays
typical 2D spectra for a conjugated polymeric J-aggregate material (H2TPPS)
taken in solution at room temperature. In this case, exciton/exciton
interaction is essentially non-existent and the resulting 2D spectral maps
show absorptive line-shapes in both the rephasing and non-rephasing signals.
This latter case is well-described using current theoretical treatments of
non-linear spectroscopy assuming that the line-shape results from stationary
fluctuations about an average transition frequency.

Figure 7:  (a) Examples of experimental (upper line) and simulated (lower
line) 2D maps obtained in the rephasing (R), non rephasing (NR), and double-
quantum (2Q) configurations for H2TPPS aggregates in solution at room
temperature. The R and NR maps report the real part of the signal at $t_{2}=0$
fs. The 2Q maps refer to t1 = 0. All of the maps are normalized to their
maximum. The energies of relevant states and the coordinates where the traces
shown in panels d–f are extracted are also pinpointed in the maps. Reproduced
from Ref. 48.

We point out that inorganic and hybrid inorganic-organic semiconductors, such
as GaAs quantum wells and (PEA)2PbI4 discussed in this section, Coulomb
exciton-exciton interactions are strong and EID effects are evident in the 2D
coherent spectral lineshape. Organic semiconductors display a contrasting
situation, exemplified by measurements of 2D coherent spectra in molecular J
aggregates (Fig. 7) [48]. In this work, many-body excitonic effects are
identified and analyzed, but the distinct lineshape effects in Figs. 2, 4, and
6 are not observed. This points to the highly localized nature of excitons in
organic systems, where interactions with local vibration are dominant.

In the next section, we present an overview of our recent work in developing
spectroscopic models in which the energy gaps evolve in concert with a bath of
background excitations that are both non-equilibrium and non-stationary as the
result of being incoherently pumped by a series of laser pulses. Here, we
develop the theory starting from a microscopic/many-body description of a
system of excitons interacting via long-range Coulomb interactions and coupled
to a dissipative environment. The theory is developed by deriving stochastic
Langevin equations for the excitons and then using these derive spectral
responses in the mean-field limit. The resulting model reduces to the well-
known Anderson-Kubo model in the limit that the excitonic dynamics are
stationary. The model provides a microscopic origin for the EID and EIS
effects in semiconductor systems. As part of our review, we work through many
of the technical details of the theory and our use of stochastic calculus to
derive analytical expressions for the spectral responses.

##  3 NONLINEAR COHERENT SPECTROSCOPY OF NONSTATIONARY SYSTEMS

###  3.1 Optical Bloch Equations

The many-body exciton scattering signatures in the 2D coherent linshape,
reported in Fig. 6 for GaAs quantum wells, were rationalized by numerical
simulation based on modified optical Bloch equations (OBE) [22]. For a two-
level system that includes both excitation-induced dephasing (EID) and
excitation-induced shift (EIS) effects, the off-diagonal term of the density
matrix, $\rho_{12}$, follows the following equation of motion:

$\displaystyle\dot{\rho}_{12}$
$\displaystyle=-\left[(\gamma_{0}+\gamma^{\prime}N\rho_{22})-i(\omega_{0}+\omega^{\prime}N\rho_{22})\right]\rho_{12}+\frac{i}{\hbar}\vec{\mu}_{12}\cdot\vec{E}(t)(\rho_{22}-\rho_{11})$
$\displaystyle=i\left[\left(\omega_{0}+i\gamma_{0}\right)+\left(\omega^{\prime}+i\gamma^{\prime}\right)N\rho_{22}\right]\rho_{12}+\frac{i}{\hbar}\vec{\mu}_{12}\cdot\vec{E}(t)(\rho_{22}-\rho_{11}),$
(2)

where $\gamma_{0}$ is the natural dephasing rate, $\rho_{11}$ and $\rho_{22}$
are the ground and excited state populations connected by coherence term
$\rho_{12}$, $N$ is the number density of chromophores and $\gamma^{\prime}$
and $\omega^{\prime}$ characterize the collision rate and collective
interactions within the excited state population. The last term corresponds to
the driving field of the laser and dipole coupling between the ground and
excited states. As the OBE system evolves, the total dephasing rate,
$(\gamma_{0}+\gamma^{\prime}N\rho_{22}(t))$, and the phase oscillation
frequency, $(\omega_{0}+\omega^{\prime}N\rho_{22}(t))$, both depend upon the
fraction of chromophores in the excited state at time $t$.

Figure 8: Calculated real spectra for the rephasing sequence. (a) is based on
a simple V system without any many-body interactions, while excitation-induced
dephasing is included in (b). Reproduced with permission from Ref. 22

Fig. 8 shows the predicted four-wave mixing signals for a two-level system
with and without the term attributed to the EID component. While the simple
OBE approach does capture the narrowing, shift, and asymmetry of the
lineshape, it fails to capture the phase scrambling that is clearly observed
in the experimental signals in Fig. 6. The OBE real signals are clearly
absorptive, rather than dispersive.

However, the decomposition of a 2D spectrum into its real and imaginary parts
depends on the techniques applied e.g., through the comparison to an
independent spectrally resolved differential transmission measurement. In
other words, the feature of dispersive lineshape may be hidden in the
imaginary part of the spectrum. On the other hand, from the Green’s function
approach, the EID and EIS are attributed respectively to the real and
imaginary parts of the exciton self-energy renormalization. Therefore, the
term of $\omega^{\prime}+i\gamma^{\prime}$ in Eq. 2 cannot be treated
separately in a quantum mechanical theory. As we will illustrate later in our
model, including the many-body interaction in the Hamiltonian leads to the
excitation-induced dephasing and frequency shift effect in the spectral
signals simultaneously.

What is then desired is an approach that incorporates the many-body dynamics
of a dark, non-optical population that co-evolves with the optical signals. In
this section, we accomplish this via the use of a stochastic line-shape
approach, taken in the limit that the non-optical population contributes non-
equilibrium and non-stationary contributions to the fluctuations of the
optical energy gap.

###  3.2 Stochastic Many-body Processes

Our model is initiated by assuming that at $t=0$ a non-stationary population
of background excitations is created by a broad-band laser excitation. This
physical picture is sketched in Fig. 1. In the current work, excitation occurs
with a sequence of phase-matched and time-ordered femtosecond pulses used to
measure a coherent nonlinear excitation spectrum, and the excitons produced
and measured via a well defined coherent pathway (see Fig. 10 for the relevant
ones in this work) are assumed to scatter elastically with their incoherent
counterparts — excitons that are produced by the pulse sequence but have no
phase relationship to those that produce signal in our experiments. The
initial background population can be characterized by an average population
$N_{0}$ and variance $\sigma_{N_{0}}^{2}$ both of which depend upon the
excitation pulse as well as the density of states of the material. Optical
excitations at $k=0$ evolve in concert with a non-stationary ($k\neq 0$)
background of excitations in which the interaction determined by a screened
Coulomb potential giving rise to a noisy driving term that effectively
modulates the exciton energy gap.

Here we consider the case where we have an ensemble of bosonic excitons
described by a Hamiltonian written in second-quantized form as

$\displaystyle
H=\sum_{k}\hbar\omega_{k}a_{k}^{\dagger}a_{k}+\frac{1}{2}\sum_{kk^{\prime}q}V_{q}a^{\dagger}_{k+q}a^{\dagger}_{k^{\prime}-q}a_{k^{\prime}}a_{k},$
(3)

where $V=L^{3}$ is the unit volume and $V_{q}$

$\displaystyle V_{q}=\frac{1}{(2\pi)^{3}}\int V({\bf r})e^{-i{\bf q}\cdot{\bf
r}}d{\bf r}$ (4)

is the Fourier component of the many-body interaction potential. We now
collect the $k\neq 0$ terms by keeping those interacting with $k=0$ excitons
and involving no more than two $k\neq 0$ states.

$\displaystyle H$
$\displaystyle=\hbar\omega_{0}a^{\dagger}_{0}a_{0}+\sum_{q\neq
0}\hbar\omega_{q}a_{q}^{\dagger}a_{q}+\frac{V_{0}}{2}a^{\dagger}_{0}a^{\dagger}_{0}a_{0}a_{0}$
$\displaystyle+a^{\dagger}_{0}a_{0}\left[{2V_{0}}\sum_{q\neq
0}(a^{\dagger}_{q}a_{q})\right]$
$\displaystyle+a^{\dagger}_{0}a^{\dagger}_{0}\left[\frac{V_{0}}{2}\sum_{q\neq
0}a_{q}a_{-q}\right]+a_{0}a_{0}\left[\frac{V_{0}}{2}\sum_{q\neq
0}a_{q}^{\dagger}a_{-q}^{\dagger}\right]$ (5)

and focus only the $k=0$ term

$\displaystyle H/\hbar$
$\displaystyle=\omega_{0}a^{\dagger}_{0}a_{0}+\Omega\hat{A}^{\dagger}\hat{A}$
$\displaystyle+\frac{\gamma_{1}}{2}\left(a^{\dagger}_{0}a^{\dagger}_{0}a_{0}a_{0}+4a^{\dagger}_{0}a_{0}\hat{N}+a^{\dagger}_{0}a^{\dagger}_{0}\hat{A}\hat{A}+a_{0}a_{0}\hat{A}^{\dagger}\hat{A}^{\dagger}\right)$
(6)

where the $\hat{A}$ , $\hat{A}^{\dagger}$, and $\hat{N}$ operators are
collective bath operators defined by inspection of Eq. 5.
$\hbar\gamma_{1}=V_{0}$ is the exciton-exciton interaction, which we obtain
from the $s$-wave scattering length $a$ and reduced mass $\mu$ within the Born
approximation [49]

$\displaystyle\gamma_{1}=\frac{4\pi\hbar a}{\mu}.$ (7)

This assumption does not rely upon the specific form of the exciton-exciton
interaction, only that it be of finite range. In the current context, this
interaction will be due to Coulomb-mediated exciton-exciton scattering that
gives rise to EID [27]. However, it is possible that each distinct exciton
within the family of the 2D perovskite system considered here [31] have a
distinct and unique value of $\gamma_{1}$, as we reported in ref. 27, where we
demonstrated distinct Coulomb screening of different exciton polarons. For
purposes of our theoretical model, we assume that the system has a single
exciton species that is susceptible to many-body scattering and therefore EID
mediated via $\gamma_{1}$.

###  3.3 Exciton scattering contributions in the mean-field (Hartree) limit
and quantum Langevin equations.

The term involving $\hat{N}$ introduces mean-field between the $k=0$ excitons
and the net population of the $k\neq 0$ excitons. This term introduces energy
fluctuation simply due to scattering of the $k\neq 0$ population from the
$k=0$ population. The other two terms give rise to fluctuations/dissipation
due to exciton pair creation/annihilation. For the moment, we shall neglect
these terms, but will return to discussing them in a later section.

We now assume that the collective operators for the background (dark) excitons

$\gamma_{1}\hat{A}^{\dagger}\hat{A}=\sum_{q\neq
0}(V_{0}+V_{q})a^{\dagger}_{q}a_{q}$

are coupled to an ensemble of otherwise unspecified bath oscillators
$b_{i}/b_{i}^{\dagger}$ as described by an auxiliary Hamiltonian of the form:

$\displaystyle H_{aux}=$
$\displaystyle\hbar\Omega\left(\hat{A}^{\dagger}\hat{A}+\frac{1}{2}\right)+\sum_{i}\hbar\omega_{i}\left(\hat{b}_{i}^{\dagger}\hat{b}_{i}+\frac{1}{2}\right)$
$\displaystyle+\sum_{i}\left(g_{i}\hat{A}^{\dagger}\hat{b}_{i}+g_{i}^{*}\hat{A}\hat{b}_{i}^{\dagger}\right).$
(8)

For convenience we rotate the operators in the Heisenberg representation
(denoted by subscript H) according to their respective frequency

$\displaystyle\hat{A}(t)$ $\displaystyle=$ $\displaystyle\hat{A}_{\rm
H}(t)e^{-i\Omega t}$ (9) $\displaystyle\hat{b}_{i}(t)$ $\displaystyle=$
$\displaystyle\hat{b}_{i,{\rm H}}(t)e^{-i\omega_{i}t}$ (10)

so that $\hat{A}$ and $\hat{b}_{i}$ evolve only under the influence of the
interaction. Then we have the Heisenberg equation of motion for the rotated
exciton operator

$\frac{\rm d}{{\rm d}t}\hat{A}(t)=-\int_{0}^{t-t_{0}}{\rm
d}\tau\kappa(\tau)\hat{A}(t-\tau)+\hat{F}(t),$ (11)

where

$\displaystyle\kappa(\tau)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar^{2}}\sum_{i}|g_{i}|^{2}e^{i(\Omega-\omega_{i})\tau}$
(12) $\displaystyle\hat{F}(t)$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\sum_{i}g_{i}\hat{b}_{i}(t_{0})e^{i(\Omega-\omega_{i})t}.$
(13)

Considering that bath oscillator frequencies $\omega_{i}$ cover a wide range
and $|g_{i}|^{2}$ may vary slowly with $\omega_{i}$, the oscillating
exponentials in $\kappa(\tau)$ interfere destructively for $\tau>0$ and
$\kappa(\tau)$ becomes negligible when $\tau\gg\tau_{c}$, where $\tau_{c}$ is
the correlation time of the bath. In addition, the time scale (damping time)
of the excitons is often much greater than $\tau_{c}$, so that one may replace
$\hat{A}(t-\tau)$ by $\hat{A}(t)$ from the integral. Eq.(11) can be rewritten
as

$\frac{\rm d}{{\rm
d}t}\hat{A}(t)=-\left(\frac{\Gamma}{2}+i\Delta\right)\hat{A}(t)+\hat{F}(t),$
(14)

in which

$\displaystyle\Gamma$ $\displaystyle=$
$\displaystyle\frac{2\pi}{\hbar^{2}}\sum_{i}|g_{i}|^{2}\delta(\Omega-\omega_{i}),$
(15) $\displaystyle\Delta$ $\displaystyle=$
$\displaystyle{\cal{P}}\sum_{i}\frac{|g_{i}|^{2}}{\hbar^{2}\left(\Omega-\omega_{i}\right)}$
(16)

describe the spontaneous emission rate and the spontaneous radiative shift of
the system operator, respectively.

The bath operator $\hat{F}(t)$ has the following properties

$\displaystyle\left<\hat{F}(t)\right>$ $\displaystyle=$ $\displaystyle{\rm
Tr}\left[\sigma_{A}\sigma_{B}\hat{F}(t)\right]=0,$ (17)
$\displaystyle\left<\hat{F}(t^{\prime})\hat{F}(t)\right>$ $\displaystyle=$
$\displaystyle\left<\hat{F}^{\dagger}(t^{\prime})\hat{F}^{\dagger}(t)\right>=0,$
(18) $\displaystyle\left<\hat{F}^{\dagger}(t^{\prime})\hat{F}(t)\right>$
$\displaystyle=$ $\displaystyle\sum_{i}\frac{1}{\hbar^{2}}|g_{i}|^{2}\langle
n_{i}\rangle e^{i(\Omega-\omega_{i})(t-t^{\prime})},$ (19)
$\displaystyle\left<\hat{F}(t^{\prime})\hat{F}^{\dagger}(t)\right>$
$\displaystyle=$
$\displaystyle\sum_{i}\frac{1}{\hbar^{2}}|g_{i}|^{2}\left(\langle
n_{i}\rangle+1\right)e^{i(\Omega-\omega_{i})(t-t^{\prime})},$ (20)

where $n_{i}=b_{i}^{\dagger}(t_{0})b_{i}(t_{0})$. Therefore $\hat{F}(t)$ can
be considered as a Langevin force fluctuating around its zero average value,
and Eq.(14) is a Langevin equation. The diffusion coefficients are defined as

$\displaystyle 2D_{N}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}{\rm
d}\tau\left<\hat{F}^{\dagger}(t-\tau)\hat{F}(t)\right>$ (21) $\displaystyle
2D_{A}$ $\displaystyle=$ $\displaystyle\int_{-\infty}^{+\infty}{\rm
d}\tau\left<\hat{F}(t)\hat{F}^{\dagger}(t-\tau)\right>,$ (22)

with subscripts $N$ and $A$ denote for the normal and anti-normal order of
$\hat{F}^{\dagger}$ and $\hat{F}$. Furthermore,

$\displaystyle 2D_{N}=$ $\displaystyle\Gamma^{\prime}$
$\displaystyle=\Gamma\left<n(\Omega)\right>$ (23) $\displaystyle 2D_{A}=$
$\displaystyle\Gamma^{\prime}+\Gamma$
$\displaystyle=\Gamma\left(1+\left<n(\Omega)\right>\right),$ (24)

where $\langle n(\Omega)\rangle$ is the average number of quanta of the bath
modes having the same frequency $\Omega$ as the system.

In the mean-field Hamiltonian Eq.(6), the background exciton population
operator $\hat{N}(t)=\hat{A}^{\dagger}(t)\hat{A}(t)$ is of our interest. Using
Eq.(14) we deduce the Heisenberg equation of motion

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\left[\hat{A}^{\dagger}(t)\hat{A}(t)\right]=-\left[\Gamma+i(\Delta-\Delta^{*})\right]\hat{A}^{\dagger}(t)\hat{A}(t)+\hat{F}^{\dagger}(t)\hat{A}(t)+\hat{A}^{\dagger}(t)\hat{F}(t).$
(25)

Integrating the differential equation Eq.(14) from $t_{0}$ to $t^{\prime}$

$\hat{A}(t^{\prime})=\hat{A}(t_{0})e^{-(\Gamma/2+i\Delta)(t^{\prime}-t_{0})}+\int_{t_{0}}^{t^{\prime}}{\rm
d}t^{\prime\prime}\hat{F}(t^{\prime\prime})e^{-(\Gamma/2+i\Delta)(t^{\prime}-t^{\prime\prime})}.$
(26)

For $t_{0}\rightarrow-\infty$, the first term is negligible and we have, by
multiplying both sides by $\hat{F}^{\dagger}(t)$,

$\left<\hat{F}^{\dagger}(t)\hat{A}(t^{\prime})\right>=\int_{t_{0}}^{t^{\prime}}{\rm
d}t^{\prime\prime}\left<\hat{F}^{\dagger}(t)\hat{F}(t^{\prime\prime})\right>e^{-(\Gamma/2+i\Delta)(t^{\prime}-t^{\prime\prime})}.$
(27)

$\hat{F}^{\dagger}(t)$ is only correlated with $\hat{A}(t^{\prime})$ in an
interval of $t\in(t^{\prime},t^{\prime}-\tau_{c})$ and
$\left<\hat{F}^{\dagger}(t)\hat{F}(t^{\prime})\right>$ varies much more
rapidly with $t-t^{\prime}$ than the exponential, therefore

$\displaystyle\left<\hat{F}^{\dagger}(t)\hat{A}(t)\right>$
$\displaystyle\approx$ $\displaystyle\int_{t_{0}}^{t}{\rm
d}t^{\prime\prime}\left<\hat{F}^{\dagger}(t)\hat{F}(t^{\prime\prime})\right>$
(28) $\displaystyle\left<\hat{A}^{\dagger}(t)\hat{F}(t)\right>$
$\displaystyle\approx$ $\displaystyle\int_{t_{0}}^{t}{\rm
d}t^{\prime\prime}\left<\hat{F}^{\dagger}(t^{\prime\prime})\hat{F}(t)\right>.$
(29)

Let $\tau=t-t^{\prime\prime}$, hence we have

$\left<\hat{F}^{\dagger}(t)\hat{A}(t)\right>+\left<\hat{A}^{\dagger}(t)\hat{F}(t)\right>=\int_{-(t-t_{0})}^{+(t-t_{0})}{\rm
d}\tau\left<\hat{F}^{\dagger}(t-\tau)\hat{F}(t)\right>=2D_{N}.$ (30)

The background exciton population operator
$\hat{N}(t)=\hat{A}^{\dagger}(t)\hat{A}(t)$ satisfies

$\displaystyle\frac{\rm d}{{\rm
d}t}\left<\hat{A}^{\dagger}(t)\hat{A}(t)\right>=$
$\displaystyle\left<\frac{\rm d}{{\rm
d}t}\left[\hat{A}^{\dagger}(t)\right]\hat{A}(t)\right>+\left<\hat{A}^{\dagger}(t)\frac{\rm
d}{{\rm d}t}\left[\hat{A}(t)\right]\right>$ $\displaystyle=$
$\displaystyle-\Gamma\left<\hat{A}^{\dagger}(t)\hat{A}(t)\right>+\left<\hat{F}^{\dagger}(t)\hat{A}(t)\right>+\left<\hat{A}^{\dagger}(t)\hat{F}(t)\right>$
$\displaystyle=$
$\displaystyle-\Gamma\left<\hat{A}^{\dagger}(t)\hat{A}(t)\right>+2D_{N}.$ (31)

The nonzero average value of
$\langle\hat{F}^{\dagger}(t)\hat{A}(t)+\hat{A}^{\dagger}(t)\hat{F}(t)\rangle$
prevents us from treating it as a Langevin force. However, the equation of
motion Eq.(31) can still be described by a generalized Ornstein-Uhlenbeck
equation with a drift $2D_{N}/\Gamma$

${\rm d}N(t)=-\Gamma\left[N(t)-\frac{2D_{N}}{\Gamma}\right]{\rm d}t+\sigma{\rm
d}W(t),$ (32)

where ${\rm d}W(t)$ represents a Wiener process. Because the background
population operator $\hat{N}(t)$ and $k=0$ exciton operator
$a_{0}/a_{0}^{\dagger}$ evolves independently, hereafter we treat
$N(t)=\Tr[\hat{N}(t)]$ as a function rather than an operator. The formal
solution is

$N(t)=N(0)e^{-\Gamma t}+\frac{2D_{N}}{\Gamma}\left(1-e^{-\Gamma
t}\right)+\sigma\int_{0}^{t}e^{-\Gamma(t-s)}{\rm d}W_{s}.$ (33)

The covariance function does not depend on the drift,

$\displaystyle\mathrm{Cov}[N(s),N(t)]=\sigma_{N_{0}}^{2}e^{-\Gamma(s+t)}+\frac{\sigma^{2}}{2\Gamma}\left[e^{-\Gamma|t-s|}-e^{-\Gamma(t+s)}\right].$

However, the average value of the stationary state no longer vanishes but is
determined by the drift

$\left<N(t)\right>_{\rm ss}=\lim_{t\rightarrow\infty}\left[N(0)e^{-\Gamma
t}+\frac{2D_{N}}{\Gamma}\left(1-e^{-\Gamma
t}\right)\right]=\frac{2D_{N}}{\Gamma}.$ (34)

The drift term can also be approximated by $2D_{N}/\Gamma=\langle
n(\Omega)\rangle$, where $\langle n(\Omega)\rangle$ is the average number of
quanta of the bath modes having the same frequency $\Omega$ as the system.
This indicates the equilibrium between the system and the bath.

We can write the damping velocity as an operation denoted by ${\cal
D}\left(\hat{A}(t)\right)$

$\displaystyle{\cal D}\left(\hat{A}(t)\right)$ $\displaystyle=$
$\displaystyle-\left(\frac{\Gamma}{2}+i\Delta\right)\hat{A}(t)$ (35)
$\displaystyle{\cal D}\left(\hat{A}^{\dagger}(t)\right)$ $\displaystyle=$
$\displaystyle-\left(\frac{\Gamma}{2}-i\Delta\right)\hat{A}^{\dagger}(t)$ (36)
$\displaystyle{\cal D}\left(\hat{A}^{\dagger}(t)\hat{A}(t)\right)$
$\displaystyle=$
$\displaystyle-\Gamma\hat{A}^{\dagger}(t)\hat{A}(t)+\Gamma^{\prime}.$ (37)

Then we have

$2D_{N}=\left<{\cal D}\left(\hat{A}^{\dagger}\hat{A}\right)-{\cal
D}\left(\hat{A}^{\dagger}\right)\hat{A}-\hat{A}^{\dagger}{\cal
D}\left(\hat{A}\right)\right>=\Gamma^{\prime}=\Gamma\left<n(\Omega)\right>,$
(38)

which can be considered as a generalization of the Einstein relation
$(D=Mk_{B}T\gamma)$ between the diffusion coefficient $D$ and the decay rate
$\Gamma$.

Similarly, the diffusion coefficient $D_{A}$ is related to
$\hat{A}(t)\hat{A}^{\dagger}(t)$ by

$\frac{\rm d}{{\rm
d}t}\left<\hat{A}(t)\hat{A}^{\dagger}(t)\right>=-\Gamma\left<\hat{A}(t)\hat{A}^{\dagger}(t)\right>+2D_{A}.$
(39)

Since $\left<\left[\hat{A}(t),\hat{A}^{\dagger}(t)\right]\right>=1$, we
conclude

$D_{A}-D_{N}=\frac{\Gamma}{2}.$ (40)

Note that $D_{N}$ and $D_{A}$ are originated from the _non-Hamiltonian_ part
of the damping/relaxation. The radiative shift term of $\Delta$ can be
absorbed into an effective Hamiltonian thus it is not included in the
equations of motion for $A^{\dagger}(t)A(t)$ and $A(t)A^{\dagger}(t)$.

For a stationary background population, i.e. $\langle N(t)\rangle=0$ the
covariance evolves according to

$\left<N(t)N(s)\right>=\left<N(t-s)N(0)\right>=\frac{\sigma^{2}}{2\Gamma}\exp\left(-\Gamma|t-s|\right).$

In this limit, our model reduces to the Anderson-Kubo model in which the
frequency fluctuates about a stationary average according to an Ornstein-
Uhlenbeck process. In this case, the population relaxation time in our model
is equivalent to the correlation time in Anderson-Kubo and this gives the rate
at which the environment relaxes back to its stationary average given a small
push. Moreover, the fluctuation amplitude, $\Delta\omega^{2}$, in Anderson-
Kubo is equivalent to $\sigma^{2}/2\Gamma$ in our model. As we shall show,
what appears at first to be a simple modification to the dynamics of a system
has significant implications in terms of the non-linear spectral response of
the system.

At time $t=0$, we push the background population significantly away from the
steady-state distribution to an initial value of $\langle N(0)\rangle=N_{0}$,
the population evolves as

$\displaystyle N(t)=N(0)e^{-\Gamma t}+\sigma\int_{0}^{t}e^{-\Gamma(t-s)}{\rm
d}W(s).$ (41)

and

$\displaystyle\langle N(t)\rangle=e^{-\Gamma t}N_{0},$ (42)

where $N_{0}$ is the mean number of background excitations present at time
$t=0$. [] Itô lemma $(dW(t))^{2}=dt$ where $W(t)$ is the Wiener process.

Normal calculus Following the usual rules of calculus, we would write

$\displaystyle df$ $\displaystyle=f[x(t)+dx(t)]$ $\displaystyle-f[x(t)]$
$\displaystyle=f^{\prime}dx.$

taking $(dx)^{2}=0$.

Itô identityHowever, if $f(x(t))$ is a function of a stochastic variable with
$dx=adt+bdW$ then $f$ satisfies the SDE:

$\displaystyle df$ $\displaystyle=f[x(t)+dx(t)]$ $\displaystyle-f[x(t)]$
$\displaystyle=(af^{\prime}+\frac{1}{2}bf^{\prime\prime})dt$
$\displaystyle+bf^{\prime}dW.$

In principle, there will be a distribution about this mean characterized by a
variance $\sigma^{2}_{N_{0}}$. As a result, we break reversibility and the
time symmetry of the correlation functions. Mathematically, this means that
$\langle N(t)N(s)\rangle\neq\langle N(t-s)N(0)\rangle$ since the choice of
initial time is no longer arbitrary.

In Ref. [50, 41] we used Itô calculus to evaluate these correlation functions.
From a practical point of view, the Itô calculus is a tool for manipulating
stochastic processes that are closely related to Brownian motion and Itô’s
lemma allows us to easily perform noise-averaged interactions. For the model
at hand, the covariance of $N(s)$ and $N(t)$ is given by

$\displaystyle\mathrm{Cov}\left[N(s),N(t)\right]=$
$\displaystyle\left<(N(s)-\langle N(s)\rangle)(N(t)-\langle
N(t)\rangle)\right>$
$\displaystyle=\frac{\sigma^{2}}{2\Gamma}\left[e^{-\Gamma|t-s|}-e^{-\Gamma(t+s)}\right]+\sigma_{N_{0}}^{2}e^{-\Gamma(s+t)},$
(43)

with $\sigma_{N_{0}}^{2}$ being the variance of $N(0)$. Similarly, the
variance

$\displaystyle\mathrm{Var}[N(t)]=\left(\sigma_{N_{0}}^{2}-\frac{\sigma^{2}}{2\Gamma}\right)e^{-2\Gamma
t}+\frac{\sigma^{2}}{2\Gamma}$ (44)

also depends upon the initial fluctuation in the background population.
Mathematically, the Fourier transform of the kernel of the integral in Eq. 41
provides the spectral density of the noisy process. In fact, a trivial
modification of the approach would be to replace the kernel in Eq. 41 with
another kernel reflecting a more complex spectral density. The resulting
expressions for the responses will be more complex indeed. However, Itô’s
lemma provides a tractable route for computing the necessary response
functions.

###  3.4 Predictions from the stochastic model

Having established the mathematical model, let us briefly recapitulate some of
its features. First, we started by assuming that the background population
dynamics give rise to a stochastic process $N(t)$ that enters into the
Heisenberg equations of motion for the system operators. In particular, we
assumed that $N(t)$ corresponds to an overdamped Brownian oscillator and that
at time $t=0$ there is a non-stationary population of background excitations.
These two mathematical assumptions can be relaxed to some extent if one has a
more detailed description of the spectral density of the background process
and the initial background population. Secondly, we assume that averages over
exponential terms can be evaluated using the cumulant expansion. What then
follows are the mathematical consequences as expressed in terms of the
spectral responses of the model.

####  3.4.1 Linear response

The linear response for optical excitation is given by

$\displaystyle S^{(1)}(t)$
$\displaystyle=\frac{i}{\hbar}\left<\hat{\mu}(t)[\hat{\mu}(0),\rho(-\infty)]\right>,$
(45)

where $\hat{\mu}(t)=\mu(\hat{a}_{0}^{\dagger}(t)+\hat{a}_{0}(t))$ is the
excitonic transition dipole operator and $\rho(-\infty)$ is the initial
density matrix. The absorption spectrum is obtained by Fourier transformation.

Averaging over the fluctuations generates terms involving cumulants of the
background noise, which result in terms such as

$\displaystyle\left\langle\exp\left[i2\gamma_{1}\int_{0}^{t}N(\tau){\rm
d}\tau\right]\right\rangle\approx
e^{i2\gamma_{1}g_{1}(t)-2\gamma_{1}^{2}g_{2}(t)},$ (46)

where $\langle\cdots\rangle$ denotes averaging over noise. Note that the
exciton interaction strength is $2\gamma_{1}$ in Eq. 6. Here, the first
cumulant $g_{1}(t)$ gives rise to a characteristic frequency shift as the
background population decays:

$\displaystyle g_{1}(t)=\int_{0}^{t}\langle N(\tau)\rangle{\rm
d}\tau=\frac{N_{0}}{\Gamma}\left(1-e^{-\Gamma t}\right),$ (47)

and

$\displaystyle
g_{2}(t,t^{\prime})=\int_{0}^{t}\int_{0}^{t^{\prime}}\mathrm{Cov}$
$\displaystyle\left[N(\tau),N(\tau^{\prime})\right]{\rm d}\tau^{\prime}{\rm
d}\tau=$
$\displaystyle\frac{\sigma^{2}}{2\Gamma^{3}}\left[2\Gamma\mathrm{min}(t,t^{\prime})+2e^{-\Gamma
t}+2e^{-\Gamma
t^{\prime}}-e^{-\Gamma|t^{\prime}-t|}-e^{-\Gamma(t^{\prime}+t)}-2\right]$ (48)
$\displaystyle+\frac{\sigma_{N_{0}}^{2}}{\Gamma^{2}}\left[e^{-\Gamma(t+t^{\prime})}-e^{-\Gamma
t}-e^{-\Gamma t^{\prime}}+1\right].$

When the two time limits are the same, this reduces to

$\displaystyle g_{2}(t)$ $\displaystyle=\int_{0}^{t}\int_{0}^{t}{\rm
Cov}[N(\tau),N(\tau^{\prime})]{\rm d}\tau{\rm d}\tau^{\prime}$
$\displaystyle=\frac{\sigma^{2}}{2\Gamma^{3}}\left(2\Gamma t+4e^{-\Gamma
t}-e^{-2\Gamma
t}-3\right)+\frac{\sigma_{N_{0}}^{2}}{\Gamma^{2}}\left(1-e^{-\Gamma
t}\right)^{2}.$ (49)

In Fig. 9 we highlight some of the key physical effects that can appear in the
linear absorption spectra based upon our model. These effects are consistent
with experimental observations and theoretical models of 2D semiconductors and
transition metal dichalcogenides [27, 38, 29]. Fig. 9(a) displays the effect
of a non-stationary background on the linear absorption spectrum of a system.
The notable feature is the tail that extends to higher absorption energies.
The character of this tail depends most strongly upon the initial choice of
$N_{0}$ and is attributable to the $g_{1}(t)$ term in our response function
which is the time-integral over the evolving background population. This term,
as it appears in Eq. 45, produces an evolving frequency shift reflecting the
dynamical relaxation of the background. In the $S^{(1)}$ response, the
background evolution is manifest as a tail extending out to the blue.

Figure 9: The linear response function with (a) increasing background
population density $N_{0}$, and (b) different relaxation rate $\Gamma$, from
the homogeneous limit of $\Gamma=50$ meV to the inhomogeneous limit of
$\Gamma=2$ meV. (From Ref.50)

[Spectral effects that can be attributed to the non-stationary evolution.]

  1. 1.

Blocking: Increasing the initial background exciton density suppresses the
peak absorption intensity.

  2. 2.

Biexciton formation: The peak is split by $\gamma_{1}/2$ corresponding
biexciton interactions [38].

  3. 3.

(1D) Energy shift: The peak position shifts to the blue with increasing
background population due to increased Coulombic interactions.

  4. 4.

(1D) Broadening: The spectrum acquires a long tail extending to the blue due
to the dynamical evolution of the background.

  5. 5.

(2D) Phase scrambling: appears in the 2D coherent spectroscopy as an asymmetry
along the absorption axis and as phase scrambling in the rephasing and non-
rephasing signals.

  6. 6.

(2D) Excitation-induced shift: systematic shift of peak position that evolves
as the background population decays.

  7. 7.

(2D) Excitation-induced dephasing: Transient narrowing along the off-diagonal
due to decreasing rate of exciton/exciton scattering.

Fig. 9(b) shows how the linear spectra is affected by the background
relaxation rate $\gamma$ for fixed values of $N_{0}=4$. In the case of fast
background relaxation ($\Gamma=50$meV) the exciton and bi-exciton spitting is
clearly resolved and the lineshapes are Lorenzian about each peak. Decreasing
the relaxation rate $\Gamma$ produces a systematic shift towards the blue due
to the mean-field interaction between the exciton and the background
population. This shift saturates when the peak is fully shifted by
$2V_{0}N_{0}$ and the spectral peak acquires a Gaussian form reflecting mean
$N_{0}$ and variance $\sigma_{N_{0}}^{2}$ of the initial background. In this
slow-relaxation limit, the “bright” state is simply swamped and suppressed by
the background excitation.

####  3.4.2 Two-dimensional coherent spectroscopy

In Ref. 41 we discussed the linear response of our model and its relation to
the Anderson-Kubo model. Here we shall focus solely on the higher-order
responses that reveal the dynamic evolution of the two-dimensional coherent
excitation line-shape. The third-order response involves phase-matched
interactions of the system with a sequence of three laser pulses:

$\displaystyle
S^{(3)}(\tau_{3},\tau_{2},\tau_{1})=\left(\frac{i}{\hbar}\right)^{3}\left<\mu(\tau_{3})\left[\mu(\tau_{2}),[\mu(\tau_{1}),[\mu(0),\rho(-\infty)]]\right]\right>.$
(50)

The times $0<\tau_{1}<\tau_{2}<\tau_{3}$ define the sequence of the time-
ordered interactions in Fig. 10. The expressions for these can evaluated using
the standard rules for double-sided Feynman diagrams (Fig. 10, c.f. Ref. 44)
representing various optical paths that for a given pathway take the form

$\displaystyle R_{\alpha}(\tau_{1},\tau_{2},\tau_{3})=$
$\displaystyle\left(\frac{i}{\hbar}\right)^{3}\mu^{4}(n_{0}+1)^{2}\exp\left[i(\omega_{0}+n_{0}\gamma_{1})\sum_{j=1}^{3}(\pm)_{j}\tau_{j}\right]\left<\exp\left[i2\gamma_{1}\sum_{j=1}^{3}(\pm)_{j}\int_{0}^{\tau_{j}}N(s){\rm
d}s\right]\right>$ (51) $\displaystyle=$
$\displaystyle\left(\frac{i}{\hbar}\right)^{3}\mu^{4}(n_{0}+1)^{2}\exp\left[i(\omega_{0}+n_{0}\gamma_{1})\sum_{j=1}^{3}(\pm)_{j}\tau_{j}\right]$
$\displaystyle\times\exp\left[i2\gamma_{1}\sum_{j=1}^{3}(\pm)_{j}g_{1}(\tau_{j})\right]\exp\left[-2\gamma_{1}^{2}\sum_{i,j=1}^{3}(\pm)_{i}(\pm)_{j}g_{2}(\tau_{i},\tau_{j})\right].$
(52)

The sign function $(\pm)_{j}$ takes “$+$” and “$-$” depending upon whether or
not the time-step involves an excitation or de-excitation of the system. The
prefactor $(n_{0}+1)^{2}$ is for the pathways involving only single excitation
manifold (distinguished by subscript a), it is $(n_{0}+1)(n_{0}+2)$ when
double excitation (subscript b) is involved. Fig. 10 shows the most relevant
diagrams for the rephasing and non-rephasing optical response.

Figure 10: Double-sided Feynman Diagrams for coherent response functions
(equation 52) with rephasing phase matching (top): (a) $R_{2a}$, (b) $R_{3a}$,
(c) $R_{1b}^{*}$, non-rephasing phase matching (bottom): (d) $R_{1a}$, (e)
$R_{4a}$, (f) $R_{2b}^{*}$, and double quantum phase matching: (g)
$R_{3b}^{*}$, (h) $R_{4b}$.

It is important to notice that the the exciton-exciton interaction term
$\gamma_{1}$, and hence the screening due to exciton-lattice interactions,
appears in three distinct places in the third-order responses. First, as a
frequency shift due to self-interactions between the bright excitons. Second,
as a frequency shift due to interactions of bright excitons with the evolving
background population density. Third, as the leading contribution to the
lineshape. In addition, the third term involving $g_{2}(t)$ carries the
influence of the initial conditions (via $\sigma_{N_{0}}$). The effect of
many-body exciton-exciton scattering thus leads to time-evolving EID
processes. Given these observations, we expect that the homogeneous linewidth
will evolve with population time, dictated by the evolution of $g_{2}(t)$.

Figure 11: Theoretical real and imaginary spectra, respectively, of rephasing
[(a), (b)] and nonrephasing [(c), (d)] phase matching and at population
waiting time $\tau_{p}=0$ fs. The vertical false color scale indicated to the
right if the figure is in arbitrary units. Figure 12:  (a)–(d): Real parts of
theoretical rephasing spectra at population times $\tau_{p}$ indicated at the
top of each panel. (e)–(h): Corresponding imaginary parts of the spectrum.
(i)–(l): The norm (absolute value) of the optical response. Exciton 2D
coherent lineshape contour at half-maximum intensity as a function of
population waiting time derived from the theoretical rephasing absolute
spectral evolution in Fig. 12. The center of mass and one of the principle
axes are shown for each contour.

Figs. 11 and 12 correspond to the rephasing and non-rephasing behavior of
theoretical model as parametrized to approximate the excitons in the 2D metal-
halide perovskite system studied in the experimental investigations, which we
shall describe later in this section. The parameters used to produce these
spectra are given in Table 1. The two pairs of gray dashed lines correspond to
the bare exciton energy at $\hbar\omega_{0}=2.35$ eV and the dressed exciton
energy at $\hbar\omega_{0}+\gamma_{1}=2.36$ eV. Fig. 11 gives the rephasing
(a,b) and non-rephasing (c,d) spectra computed at $\tau_{p}=0$. Two features
highlighted above are immediately striking in the modelled 2D spectra. Both
the asymmetry of the signals as well as the lineshape inversion of the real
and imaginary spectral components can be traced specifically to terms within
the response functions in Eq. 52 that depend upon the transient background
relaxation and exciton self-interactions.

Both the phasing and asymmetry evolve with increasing population time as shown
in Fig. 12(a-l). Importantly, the rephasing signal evolves being dispersive at
$\tau_{p}=0$ to absorptive at longer times. The non-rephasing signal [Fig.
12(e-h)] has complementary behavior, evolving from absorptive to dispersive.
Figs. 12(i-l) give the absolute value of the total response as it evolves over
$\tau_{p}$. The peak is displaced from the diagonal and its position as well
as the linewidth evolves over $\tau_{p}$.

In Fig. 12 we extract the contour corresponding to the half-maximum intensity
at various indicated $\tau_{p}$ population times. Superimposed over each
contour is one of the principal axes of the contour scaled according to its
magnitude. The central points are the geometric centers of contours. This
analysis clearly shows that the peak systematically narrows, rotates, and
distorts as the exciton co-evolves with the background population. Moreover,
the center peak shifts by about 10 meV towards the red in both absorption and
emission spectral dimensions as Coulombic interactions with the evolving
background are diminished. [10].

Within the stochastic-model the frequency peak evolution is due to the first
cumulant $g_{1}(\pm\tau_{1},\pm\tau_{2},\pm\tau_{3})$ which introduced a
phase-shift that depends the background population evolution as well as the
background displacements following each interaction with the laser field. The
early-time blue shift as well as more rapid dephasing arise from many-body
effects contained within $g_{1}$. As the background population decays, the
scattering effects are diminished. We note that in equation 52, if we set
$g_{1}=0$ the coherent response functions reduce to a stationary background,
and the frequency peak evolution does not occur.

Table 1: Parameters used in the theoretical model to produce Figs. 11 and 12.

Description | Symbol | Value
---|---|---
bare exciton energy | $\hbar\omega_{0}$ | 2.35 eV
noise variance | $\sigma^{2}$ | 0.0025 fs-1
relaxation rate | $\Gamma$ | 0.01 fs-1
exciton/exciton interaction | $\gamma_{1}$ | 10 meV
avg. init. background density | $N_{0}$ | 2 per unit volume
init. background variance | $\sigma_{N_{0}}$ | 0.35 per unit volume

We assume here that the initial background excitation is broad compared to its
fluctuations about a stationary state and choose parameters to best represent
the experimental conditions of an ultrafast experiment with $\sim 20$ fs
pulses.

###  3.5 Exciton/polaron formation dynamics due to exchange interactions

In deriving this model, we also assumed that an additional term corresponding
to pair creation/annihilation could be dropped from consideration. That term
takes the form

$\displaystyle H_{pair}=\sum_{q\neq
0}\gamma_{q}(a_{0}^{\dagger}a_{0}^{\dagger}a_{q}a_{-q}+a_{q}^{\dagger}a_{-q}^{\dagger}a_{0}a_{0})$
(53)

and corresponds to the Boson exchange interaction whereby momentum is
transferred within the background population ($q\neq 0$) as the result of
interaction with the optical ($q=0$) exciton. However, such exchange terms may
give important and interesting contributions to the spectral lineshape,
especially in systems in which excitons are formed near the Fermi energy. In
such systems, the exciton becomes dressed by virtual electron/hole
fluctuations about the Fermi sea producing spectral shifts and broadening of
the spectral lineshape. Such states are best described as exciton/polarons
whose wave function consists of the bare exciton/hole excitation dressed by
electron/hole fluctuations.

Ordinarily, as in the Bogoliubov treatment of Bose-Einstein condensates [51],
one makes the semiclassical approximation that the condensate population can
be taken as macroscopic and as a result, one can replace the $a_{0}$ and
$a_{0}^{\dagger}$ operators with c-number $\sqrt{N_{0}}$. In our case, we
shall continue to treat the background within the semiclassical limit and
replace $a_{q}$ and $a_{q}^{\dagger}$ with $\sqrt{N}_{q}$ and write the
coupling

$\displaystyle\gamma(t)=\sum_{q\neq 0}N_{q}\gamma_{q}\approx\gamma_{pair}N(t)$
(54)

where $\gamma_{pair}$ is the exchange coupling constant and $N(t)$ the net
background population at time $t$.

To pursue the effect of the exciton/polaron formation, we start with the basic
form of the Hamiltonian

$\displaystyle
H=\hbar\omega_{0}(a^{\dagger}a+1/2)+\hbar\gamma(t)(a^{\dagger}a^{\dagger}+aa)/2$
(55)

where $\gamma(t)$ is the coupling which we take to be an unspecified
stochastic process. One can bring $H$ into a diagonal form by unitary
transformation

$\displaystyle\tilde{H}=e^{-S}He^{S}=\hbar\tilde{\omega}_{0}(t)\left(\tilde{a}^{\dagger}\tilde{a}+1/2\right).$
(56)

with

$\displaystyle\tilde{\omega}_{0}(t)$
$\displaystyle=\sqrt{\omega_{0}^{2}-\gamma(t)^{2}}.$ (57)

However, since $\gamma(t)$ is a stochastic process, we need to use the Itô
identity to properly derive the underlying SDE for the renormalized harmonic
frequency, $\tilde{\omega}_{0}(t)$, in order to compute correlation functions.

In the regime of weak pair-excitation interaction, $\gamma/\omega_{0}\ll 1$,
the eigen-frequency can be approximated as

$\displaystyle\tilde{\omega}_{0}(t)$
$\displaystyle=\omega_{0}\sqrt{1-(\gamma/\omega_{0})^{2}}$
$\displaystyle\approx\omega_{0}\left(1-z(t)/2\right),$ (58)

where $z(t)=\gamma(t)^{2}/\omega_{0}^{2}$. Therefore, $\sqrt{z}$ represents
the coupling strength of the pair-excitation relative to the excitation
frequency. For the moment, we leave the stochastic variable unspecified and
find the linear response function

$\displaystyle S^{(1)}(t)=$
$\displaystyle\frac{i}{\hbar}\langle[\hat{\mu}(t),\hat{\mu}(0)]\rho(-\infty)\rangle$
$\displaystyle=$
$\displaystyle\frac{i}{\hbar}\mu^{2}\left<\left[\tilde{a}^{\dagger}(t),\tilde{a}_{0}\right]\rho(-\infty)-{\rm
c.c}\right>$ $\displaystyle=$
$\displaystyle\frac{2\mu^{2}}{\hbar}\imaginary\left<\exp(i\omega_{0}t)\exp\left[-\frac{i\omega_{0}}{2}\int_{0}^{t}z(\tau){\rm
d}\tau\right]\right>$ (59) $\displaystyle=$
$\displaystyle\frac{2\mu^{2}}{\hbar}\imaginary\left\\{\exp(i\omega_{0}t)\exp\left[\sum_{n=1}^{\infty}\frac{(-i\omega_{0}/2)^{n}}{n!}\left<\left(\int_{0}^{t}z(\tau){\rm
d}\tau\right)^{n}\right>_{\rm c}\right]\right\\}$ (60)

in the form of cumulant expansion, where $\left<x^{n}\right>_{\rm c}$ denotes
the $n$-th cumulant. According to the theorem of Marcinkiewicz,[52, 53] the
cumulant generating function is a polynomial of degree no greater than two to
maintain the positive definiteness of the probability distribution function.
Therefore, we truncate the cumulant expansion to the second order and write
the spectral line-shape functions $g_{1}(t)$ and $g_{2}(t)$ from the first
cumulant and second cumulants,

$\displaystyle g_{1}(t)=$ $\displaystyle\int_{0}^{t}\langle z(\tau)\rangle{\rm
d}\tau$ (61)

and

$\displaystyle g_{2}(t)=$
$\displaystyle\int_{0}^{t}\int_{0}^{t}\left<z(\tau),z(\tau^{\prime})\right>{\rm
d}\tau{\rm d}\tau^{\prime},$ (62)

respectively.

We now make the simplifying assumption that $\gamma(t)$ satisfies the
Ornstein-Uhlenbeck process, corresponding to vacuum fluctuations about bare
exciton state.

$\displaystyle{\rm d}\gamma_{t}=-\theta\gamma_{t}{\rm d}t+\sigma{\rm d}W_{t}.$
(63)

We should emphasize that this is not properly in the regime of quantum
fluctuations since we have not enforced the bosonic commutation relation
within the background in making the semi-classical ansatz. This is clearly an
avenue for future exploration. Applying the Itô identity, we arrive at a SDE
for the exciton frequency,

$\displaystyle{\rm d}z_{t}$
$\displaystyle=2\theta\left(\frac{\sigma^{2}}{2\theta\omega_{0}^{2}}-z_{t}\right){\rm
d}t+\frac{2\sigma}{\omega_{0}}\sqrt{z_{t}}{\rm d}W_{t},$ (64)

in which the relaxation rate is $2\theta$, and the drift term
$\sigma^{2}/2\theta\omega_{0}^{2}$ corresponds to the mean value of the
stationary state. The formal solution, analogous to $\gamma(t)$ as the
solution of the Ornstein-Uhlenbeck SDE, is

$\displaystyle z(t)^{1/2}$ $\displaystyle=\left[z(0)\right]^{1/2}e^{-\theta
t}+\frac{\sigma}{\omega_{0}}\int_{0}^{t}e^{-\theta(t-s)}{\rm d}W_{s}$ (65)
$\displaystyle\gamma(t)$ $\displaystyle=\gamma(0)e^{-\theta
t}+\sigma\int_{0}^{t}e^{-\theta(t-s)}{\rm d}W_{s}.$

Using Itô isometry we find the mean value

$\displaystyle\langle z(t)\rangle=z_{0}e^{-2\theta
t}+\frac{\sigma^{2}}{2\omega_{0}^{2}\theta}\left(1-e^{-2\theta t}\right),$
(66)

and the correlation function

$\displaystyle\langle
z(t),z(s)\rangle=\sigma_{z_{o}}^{2}e^{-2\theta(t+s)}+\frac{\sigma^{4}}{2\theta^{2}\omega_{0}^{4}}\left[e^{-\theta|t-s|}-e^{-\theta(t+s)}\right]^{2}+\frac{2\sigma^{2}}{\theta\omega_{0}^{2}}z_{0}e^{-\theta(t+s)}\left[e^{-\theta|t-s|}-e^{-\theta(t+s)}\right],$
(67)

which can be used to construct the $g_{1}$ and $g_{2}$ cumulants. Here the
initial distribution of $z$ is determined by the mean $z_{0}=\langle
z(0)\rangle$, and the variance
$\sigma_{z_{o}}^{2}=\left<\left(z(0)-z_{0}\right)^{2}\right>$.

$\displaystyle g_{1}(t)=$ $\displaystyle\int_{0}^{t}\langle z(\tau)\rangle{\rm
d}\tau$ $\displaystyle=$
$\displaystyle\frac{\sigma^{2}t}{2\theta\omega_{0}^{2}}+\frac{1}{2\theta}\left(z_{0}-\frac{\sigma^{2}}{2\theta\omega_{0}^{2}}\right)\left(1-e^{-2\theta
t}\right),$ (68)

The first cumulant of the model produces a red shift that more complex than
the counterpart in the simpler model where interactions between paired-
excitations are neglected. The initial frequency shift $z_{0}\omega_{0}/2$
agrees with the Anderson-Kubo theory, but converges to
$\sigma^{2}/4\theta\omega_{0}$ rather than decaying to zero. The first term
then can be considered as a correction term that accounts for the interaction
of exchange terms and leads to a constant red shift of
$\sigma^{2}/4\theta\omega_{0}$.

[Implications from the exciton/exciton exchange term]

  1. 1.

The model produces the lineshape function given by the Anderson-Kubo model in
the stationary limit, albeit with twice the coherence time;

  2. 2.

The model captures the formation of exciton/polarons as the steady-state/long-
time limit and gives an exciton/polaron reorganization energy of
$\sigma^{2}/4\theta\omega_{0}$ that reflects the coupling and spectral density
of the background.

The second cumulant, $g_{2}(t)$, evaluates to

$\displaystyle g_{2}(t)=$
$\displaystyle\int_{0}^{t}\int_{0}^{t}\left<z(\tau),z(\tau^{\prime})\right>{\rm
d}\tau{\rm d}\tau^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\sigma_{z_{o}}^{2}}{4\theta^{2}}\left(1-e^{-2\theta
t}\right)^{2}+\frac{\sigma^{4}}{8\theta^{4}\omega_{0}^{4}}\left(e^{-4\theta
t}+8\theta te^{-2\theta t}+4e^{-2\theta t}+4\theta t-5\right)$
$\displaystyle+$
$\displaystyle\frac{\sigma^{2}}{2\theta^{3}\omega_{0}^{2}}z_{0}\left(1-4\theta
te^{-2\theta t}-e^{-4\theta t}\right)$ (69)

If we compare this to the spectral function derived above

$\displaystyle g_{2}^{\rm
EID}(t)=\frac{\sigma_{\gamma_{o}}^{2}}{\theta^{2}}\left(1-e^{-\theta
t}\right)^{2}+\frac{\sigma^{2}}{2\theta^{3}}\left(2\theta t+4e^{-\theta
t}-e^{-2\theta t}-3\right).$ (70)

the first term is recovered; however, the present model provides a more
sophisticated description of the dependency on the initial average $z_{0}$. In
the limiting case of stationary state where
$\sigma_{\gamma_{o}}^{2}=\sigma^{2}/2\theta$, the second cumulant in the
present model turns into

$\displaystyle
g_{2}(t)=\frac{\sigma^{4}}{4\theta^{4}\omega_{0}^{4}}\left(e^{-2\theta
t}+2\theta
t-1\right)+\frac{\sigma^{2}\gamma_{0}^{2}}{\theta^{3}\omega_{0}^{4}}(e^{2\theta
t}-2\theta t-1)e^{-2\theta t},$ (71)

in which the first term reproduces the Anderson-Kubo lineshape but with half
correlation time $\tau_{c}=(2\theta)^{-1}$ compared to that of the Anderson-
Kubo theory $\theta^{-1}$. Furthermore, the second term gives the line
broadening due to the initial average of the background exciton population,
$\gamma_{0}^{2}$, which only results in a frequency shift in our previous
model.

Figure 13: Real spectra of rephasing (a) and non-rephasing (b) signal computed
by the exciton/polaron formation model.

####  3.5.1 Effect on 2D spectroscopy

The inhomogeneous and homogeneous contributions to the lineshape can be
separated using 2D coherent spectroscopic methods. [42, 54, 43, 55, 35] In
most molecular applications of 2D spectroscopy, the evolving background plays
little to no role in the spectral dynamics. However, evolving background does
affect the spectral lineshape by mixing absorptive and dispersive features in
the real and imaginary spectral components. Generally speaking, systems
lacking background dynamics exhibit absorptive line-shapes and dispersive
lineshapes are a consequence of many-body correlations [41], consistent with
the analysis of similar measurements in semiconductor quantum wells [22].
Furthermore it is useful to compare the model presented here, which pertains
to the exciton/exciton exchange coupling, versus our previous model which did
not include this term and only considered the direct (Hartree) interaction.
For this, we compute the third-order response $S^{(3)}$ from Eq. 50 under the
impulsive/rotating-wave approximation. One easily finds the responses for the
various Liouville-space pathways take the form

$\displaystyle
R_{\alpha}(\tau_{3},\tau_{2},\tau_{1})=\left(\frac{i}{\hbar}\right)^{3}\mu^{4}\left\langle\exp\left[i\sum_{j=1}^{3}(\pm)_{j}\int_{0}^{\tau_{j}}\tilde{\omega}_{0}(\tau){\rm
d}\tau\right]\right\rangle$ (72)

where the angular brackets denote averaging over the stochastic noise term and
the $(\pm)_{j}$ corresponds to whether or not the time-step involves an
excitation (+) or de-excitation (-) of the system. The time-ordering of the
three optical pulses in the experiment and phase-matching conditions define
the specific excitation pathways, based on which photon echo
($k_{s}=-k_{1}+k_{2}+k_{3}$) and virtual echo ($k_{s}=+k_{1}-k_{2}+k_{3}$)
signals can be obtained by heterodyne detection (the fourth pulse) [43].
Equivalently, in the experiments using co-linear phase-modulated pulses,
rephasing [$-(\phi_{43}-\phi_{21})$] and non-rephasing
[$-(\phi_{43}+\phi_{21})$] signals can be measured. In the rephasing
experiment, the pulse sequence is such that the phase evolution of the
polarization after the first pulse and the third pulse are of opposite sign,
while in the non-rephasing experiment, they are of the same sign. Eq.(72) can
be evaluated by cumulant expansion and the full expressions are given in the
Appendix. Since the $\tilde{\omega}_{0}(\tau)$ corresponds to a non-stationary
process, both the lineshape functions $g_{1}$ and $g_{2}$ contribute to the
output signal.

Fig. 14 presents the 2D rephasing and non-rephasing spectra corresponding to a
single quantum state dressed by the pair-excitation terms. Focusing on the
effect of interactions of paired excitations, rather than that of the initial
condition, we set $\sigma_{\gamma_{o}}^{2}=\sigma^{2}/(2\theta)$ so that the
initial fluctuation is the same as that of the Wiener process[50]. The initial
distribution of $z(0)$ is given by that of $\gamma(0)$

$\displaystyle\omega_{0}^{2}z_{0}$
$\displaystyle=\sigma_{\gamma_{o}}^{2}+\gamma_{0}^{2},$ (73a)
$\displaystyle\omega_{0}^{4}\sigma_{z_{o}}^{2}$
$\displaystyle=2\sigma_{\gamma_{o}}^{4}+4\sigma_{\gamma_{o}}^{2}\gamma_{0}^{2}.$
(73b)

The “dispersive” lineshape is observed in the real spectra for both rephasing
and non-rephasing pulse sequences, which is a clear indication of the EID. The
center of the peak deviates from the bare exciton energy
$\hbar\omega_{0}=2.35~{}{\rm eV}$ (black dashed lines) due to the coupling
between exciton pairs. Both the absorption and emission energies shift to red
because $z(t)$ is positive by definition Eq.(58). Although the Hamiltonian is
diagonal after the exciton/polaron transformation using matrix $S$, the
diagonal peaks are off the diagonal. Noting Eqs.(58) and (66), we find that
the emission frequency shift from $\omega_{0}$ by
$-\sigma^{2}/(4\theta\omega_{0})$ (red dashed line), as long as the time scale
of the experiment is greater than the relaxation time $(2\theta)^{-1}$.
Indeed, this energy discrepancy attributed to stationary state of $z(t)$ can
be considered as the exciton/polaron dressing energy. Regarding the absorption
frequency measured by the first two pulses, because the system may not have
sufficient time to relax, we can estimate, from Eq.(66), that the shift ranges
between $\omega_{0}z_{0}/2$ and $\sigma^{2}/(4\theta\omega_{0})$. The median
$\omega_{0}z_{0}/4+\sigma^{2}/(8\theta\omega_{0})$ is shown as red dashed line
for absorption.

Figure 14: Rephasing (top) and non-rephasing (bottom) spectra at population
time 100 fs based on the SDE in Eq.(63). (a) and (d) are the real, (b) and (e)
are the imaginary, (c) and (f) are the norm of the spectra. The parameters
used in the simulation: $\hbar\omega_{0}=2.35$ eV, $\sigma=0.05~{}{\rm
fs}^{-3/2}$, $\theta=0.01~{}{\rm fs}^{-1}$, and $\gamma_{0}=0.5~{}{\rm
fs}^{-1}$.

####  3.5.2 Comparison to the Anderson-Kubo model and our previous excitation-
induced dephasing (EID) theory

The well-known Anderson-Kubo theory describes the line shape broadening with
regard to the stationary state of a random variable (usually the frequency
fluctuation, $z_{t}$ here) characterized by an Ornstein-Uhlenbeck process. The
expansion of the linear optical response function leads to the first cumulant
$g_{1}^{\rm AK}(t)=\beta t$, in which $\beta$ is the drift term, i.e., the
long-term mean value, and the second cumulant

$\displaystyle g_{2}^{\rm
AK}(t)=\frac{\sigma^{2}}{2\theta^{3}}\left(e^{-\theta t}+\theta t-1\right).$
(74)

In the short time limit $\theta t\ll 1$, the first cumulant in Eq. (68) turns
to $g_{1}(t)\approx z_{0}t$, which has the same linear form as $g_{1}^{\rm
AK}(t)$. It also agrees with our previous model in Section. 3.4.1 Eq. 47 [50,
41]

$\displaystyle g_{1}^{\rm EID}(t)=\frac{z_{0}}{2\theta}\left(1-e^{-2\theta
t}\right)\approx z_{0}t.$ (75)

It is worth noting that the relaxation rate here is $2\theta$ instead of
$\theta$, because $\gamma(t)^{2}$ is the stochastic process of interest rather
than $\gamma(t)$ characterized by the rate $\theta$. For deterministic initial
condition, $z_{0}=\gamma_{0}^{2}/\omega_{0}^{2}$. The first cumulant
$g_{1}(t)$ results in a red shift of $z_{0}\omega_{0}/2$ in the linear
spectrum in the short time limit, which is determined by the initial average
of the stochastic process $z(t)$.

When the initial fluctuation of $\gamma(t)$ obeys the same Ornstein-Uhlenbeck
process, we conclude $\sigma_{\gamma_{o}}^{2}=\sigma^{2}/2\theta$ from the
stationary state corresponding to the long-time limit where Eq. (66) turns
into ${\rm Var}[\gamma(t)]=\sigma^{2}/2\theta$. Considering Eq. (73a), we
recast Eq. 68 into

$\displaystyle
g_{1}(t)=\frac{\sigma^{2}}{2\theta\omega_{0}^{2}}t+\frac{\gamma_{0}^{2}}{2\theta\omega_{0}^{2}}\left(1-e^{-2\theta
t}\right),$ (76)

in which the second term looks similar to the $g_{1}^{\rm EID}(t)$ function in
Eq. (75). However, $z_{0}=\gamma_{0}^{2}/\omega_{0}^{2}$ is true only for the
deterministic initial condition, which is not the case in the above equation.
Eq. (75) given in our previous model leads to a time-dependent red shift that
eventually vanishes after sufficiently long time. The first term then can be
considered as a correction term that accounts for the interaction of paired-
excitation and leads to a constant red shift of
$\sigma^{2}/4\theta\omega_{0}$.

Therefore, the first cumulant of the present model produces the red shift
similar to but more complex than the counterpart in our previous model, where
interactions between paired-excitations are neglected. The initial frequency
shift $z_{0}\omega_{0}/2$ agrees with the Anderson-Kubo theory, but converges
to $\sigma^{2}/4\theta\omega_{0}$ rather than decaying to zero as in the
previous model.

Regarding the second cumulant, $g_{2}(t)$, Compared the result from the
previous model Eq. 49 to Eq. 69, the second term is recovered; however, the
present model provides a more sophisticated description of the dependency on
the initial average $z_{0}$.

In the limiting case of stationary state where
$\sigma_{\gamma_{o}}^{2}=\sigma^{2}/2\theta$, the second cumulant in the
present model turns into

$\displaystyle
g_{2}(t)=\frac{\sigma^{4}}{4\theta^{4}\omega_{0}^{4}}\left(e^{-2\theta
t}+2\theta
t-1\right)+\frac{\sigma^{2}\gamma_{0}^{2}}{\theta^{3}\omega_{0}^{4}}(e^{2\theta
t}-2\theta t-1)e^{-2\theta t},$ (77)

in which the first term reproduces the Anderson-Kubo lineshape but with half
correlation time $\tau_{c}=(2\theta)^{-1}$ compared to that of the Anderson-
Kubo theory $\theta^{-1}$. Furthermore, the second term gives the line
broadening due to the initial average of the background exciton population,
$\gamma_{0}^{2}$, which only results in a frequency shift in our previous
model.

##  4 Perspective

Stochastic models have a long and important history in the field of chemical
physics since they allow one to incorporate a trajectory-based viewpoint
directly into the dynamics. The Anderson/Kubo model was an early attempt at
providing a physical rationalization of how frequency fluctuations contribute
to the absorption and emission lineshapes of molecules in contact with a
thermal environment. [56, 33, 32] The models have been continuously improved
upon over the years, notably including more detailed descriptions of the bath
and the actual coupling mechanisms between the system and the environment that
lead to the frequency fluctuations. [57, 58, 59, 44, 60] As we have repeatedly
pointed out, such approaches assume that the environment is in a stationary
(e.g. thermal) state at time $t=0$ and does not interact with any external
stimulus over the course of the dynamics of the system. We argue here, that in
many cases one can not ignore the fact that the broad-band excitation pulses
used in contemporary ultrafast experiments can create a background gas of
excitons that can interact with an optical bright state–leading to fluency-
dependent dynamics that can be manifest in terms of spectral shifts and tails
even in linear absorption spectra. We conclude that these details can be
further revealed through 2D coherent spectroscopy, especially when paired with
a theoretical approach that accounts for the non-stationary evolution of the
background.

Here we have reviewed our approach based upon a stochastic many-body treatment
of the background and have provided a number of principal results and some
technical details of our theoretical models. Throughout, we have used the Itô
stochastic calculus approach when integrating over stochastic variables. This
mathematical method provides us with a powerful avenue for obtaining
analytical expressions for the various cumulants and correlation functions
needed to compute the spectral responses. Fortunately, the current release of
Mathematica ($v>10)$ has a powerful stochastic calculus module that can be
harnessed to evaluate both formally and numerically the cumulants and
correlation functions for simple and complex transformed processes.

[Open questions] We conclude this review by posing a number of open questions
that we are currently addressing:

  1. 1.

Our current model assumes a mean-field/semiclassical treatment. Can this be
extended to include correlations between the optical and non-optical degrees
of freedom?

  2. 2.

How does one include multiple bands or a more accurate description of the
dark-state density of states?

  3. 3.

Can a similar model be developed for fermionic degrees of freedom allowing
separate treatment of electron and hole degrees of freedom?

  4. 4.

Can the spectral response be “inverted” to reveal an underlying stochastic
model for the background spectral density?

  5. 5.

While the EID and EIS effects appear to be largely present in semiconductor
systems and quantum dots, can similar effects be observed in molecular-based
systems?

  6. 6.

In ionic semiconductors such as the Ruddlesden-Popper metal halides described
in this review, how do the EID quantum dynamics reflect the polaronic nature
of excitons? What is the effect of metal, halide, and organic cation
substitution on the homogeneousl linewidth?

  7. 7.

Can this approach be used to model the effects of dark-states present in
microcavity polaritonic systems?

  8. 8.

Can the model be extended towards the strongly quantum limit to model
interacting quantum photons?

With regards to this last open question, we recently have demonstrated how
frequency correlations between emitters can be detected using quantum photons.
In particular, we show that pairs of photons, originating in a common Fock-
state, can become entangled via interactions with quantum emitters that are
correlated only through their interaction with a mutual environment. Our
analysis indicates that the cross-correlation can be detected by observing the
change in the entanglement entropy of scattered bi-photon states. [61, 36]
This latter example is probably best taken as a thought experiment since an
experimental realization of this will certainly be non-trivial due to
contemporary difficulties in preparing and detecting high-quality bi-photon
states at a sufficient flux to produce sufficient output signal. However, the
theoretical results are highly tantalizing and we look forward to future
developments along these lines.

## DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or
financial holdings that might be perceived as affecting the objectivity of
this review.

## ACKNOWLEDGMENTS

The work at the University of Houston was funded in part by the National
Science Foundation (CHE-2102506) and the Robert A. Welch Foundation (E-1337).
The work at LANL was funded by Laboratory Directed Research and Development
(LDRD) program, 20220047DR. The work at Georgia Tech was funded by the
National Science Foundation (DMR-1904293).

## References

  * [1] Mysyrowicz A, Grun J, Levy R, Bivas A, Nikitine S. 1968. Excitonic molecule in CuC1. Phys. Lett. A 26(12):615–616
  * [2] Magde D, Mahr H. 1970. Exciton-exciton interaction in CdS, CdSe, and ZnO. Phys. Rev. Lett. 24(16):890
  * [3] Grun J, Nikitine S, Bivas A, Levy R. 1970. Luminescence of copper halides excited by a high power laser. J. Lumin. 1:241–253
  * [4] Miller R, Kleinman D, Gossard A, Munteanu O. 1982. Biexcitons in GaAs quantum wells. Phys. Rev. B 25(10):6545
  * [5] Kleinman D. 1983. Binding energy of biexcitons and bound excitons in quantum wells. Phys. Rev. B 28(2):871
  * [6] Hu Y, Koch SW, Lindberg M, Peyghambarian N, Pollock E, Abraham FF. 1990. Biexcitons in semiconductor quantum dots. Phys. Rev. Lett. 64(15):1805
  * [7] Brunner K, Abstreiter G, Böhm G, Tränkle G, Weimann G. 1994. Sharp-line photoluminescence and two-photon absorption of zero-dimensional biexcitons in a GaAs/AlGaAs structure. Phys. Rev. Lett. 73(8):1138
  * [8] Albrecht T, Bott K, Meier T, Schulze A, Koch M, et al. 1996. Disorder mediated biexcitonic beats in semiconductor quantum wells. Phys. Rev. B 54(7):4436
  * [9] Stone KW, Gundogdu K, Turner DB, Li X, Cundiff ST, Nelson KA. 2009. Two-quantum 2d ft electronic spectroscopy of biexcitons in gaas quantum wells. Science 324(5931):1169–1173
  * [10] Karaiskaj D, Bristow AD, Yang L, Dai X, Mirin RP, et al. 2010. Two-quantum many-body coherences in two-dimensional fourier-transform spectra of exciton resonances in semiconductor quantum wells. Phys. Rev. Lett. 104(11):117401
  * [11] Turner DB, Nelson KA. 2010. Coherent measurements of high-order electronic correlations in quantum wells. Nature 466(7310):1089–1092
  * [12] Schultheis L, Kuhl J, Honold A, Tu CW. 1986. Ultrafast phase relaxation of excitons via exciton-exciton and exciton-electron collisions. Phys. Rev. Lett. 57(13):1635–1638
  * [13] Honold A, Schultheis L, Kuhl J, Tu CW. 1989. Collision broadening of two-dimensional excitons in a GaAs single quantum well. Phys. Rev. B 40(9):6442–6445
  * [14] Wang H, Ferrio K, Steel DG, Hu YZ, Binder R, Koch SW. 1993. Transient nonlinear optical response from excitation induced dephasing in GaAs. Phys. Rev. Lett. 71(8):1261–1264
  * [15] Wang H, Ferrio KB, Steel DG, Berman PR, Hu YZ, et al. 1994. Transient four-wave-mixing line shapes: Effects of excitation-induced dephasing. Phys. Rev. A 49(3):R1551–R1554
  * [16] Hu YZ, Binder R, Koch SW, Cundiff ST, Wang H, Steel DG. 1994. Excitation and polarization effects in semiconductor four-wave-mixing spectroscopy. Phys. Rev. B 49(20):14382–14386
  * [17] Rappen T, Peter UG, Wegener M, Schäfer W. 1994. Polarization dependence of dephasing processes: A probe for many-body effects. Phys. Rev. B 49(15):10774–10777
  * [18] Wagner HP, Schätz A, Maier R, Langbein W, Hvam JM. 1997. Coherent optical nonlinearities and phase relaxation of quasi-three-dimensional and quasi-two-dimensional excitons in ZnSxSe1-x/ZnSe structures. Phys. Rev. B 56(19):12581–12588
  * [19] Wagner HP, Schätz A, Langbein W, Hvam JM, Smirl AL. 1999. Interaction-induced effects in the nonlinear coherent response of quantum-well excitons. Phys. Rev. B 60(7):4454–4457
  * [20] Shacklette JM, Cundiff ST. 2002. Role of excitation-induced shift in the coherent optical response of semiconductors. Phys. Rev. B 66(4):045309
  * [21] Shacklette JM, Cundiff ST. 2003. Nonperturbative transient four-wave-mixing line shapes due to excitation-induced shift and excitation-induced dephasing. J. Opt. Soc. Am. B 20(4):764–769
  * [22] Li X, Zhang T, Borca CN, Cundiff ST. 2006. Many-body interactions in semiconductors probed by optical two-dimensional fourier transform spectroscopy. Phys. Rev. Lett. 96(5):057406
  * [23] Moody G, Siemens ME, Bristow AD, Dai X, Karaiskaj D, et al. 2011. Exciton-exciton and exciton-phonon interactions in an interfacial GaAs quantum dot ensemble. Phys. Rev. B 83(11):115324
  * [24] Nardin G, Moody G, Singh R, Autry TM, Li H, et al. 2014. Coherent excitonic coupling in an asymmetric double ingaas quantum well arises from many-body effects. Phys. Rev. Lett. 112(4):046402
  * [25] Moody G, Dass CK, Hao K, Chen CH, Li LJ, et al. 2015. Intrinsic homogeneous linewidth and broadening mechanisms of excitons in monolayer transition metal dichalcogenides. Nat. Commun. 6:8315
  * [26] Martin EW, Horng J, Ruth HG, Paik E, Wentzel MH, et al. 2018. Encapsulation narrows excitonic homogeneous linewidth of exfoliated MoSe2 monolayer. ArXiv:1810.09834 [cond-mat.mtrl-sci]
  * [27] Thouin F, Cortecchia D, Petrozza A, Srimath Kandada AR, Silva C. 2019a. Enhanced screening and spectral diversity in many-body elastic scattering of excitons in two-dimensional hybrid metal-halide perovskites. Phys. Rev. Res. 1:032032
  * [28] Karki KJ, Widom JR, Seibt J, Moody I, Lonergan MC, et al. 2014. Coherent two-dimensional photocurrent spectroscopy in a pbs quantum dot photocell. Nature Communications 5(1):5869
  * [29] Katsch F, Selig M, Knorr A. 2020. Exciton-scattering-induced dephasing in two-dimensional semiconductors. Phys. Rev. Lett. 124(25):257402
  * [30] Erkensten D, Brem S, Malic E. 2020. Excitation-induced dephasing in 2D materials and van der Waals heterostructures. ArXiv:2006.08392 [cond-mat.mtrl-sci]
  * [31] Srimath Kandada AR, Silva C. 2020. Exciton polarons in two-dimensional hybrid metal-halide perovskites. J. Phys. Chem. Lett. 11(9):3173–3184
  * [32] W. Anderson P. 1954. A mathematical model for the narrowing of spectral lines by exchange or motion. Journal of the Physical Society of Japan 9(3):316–339
  * [33] Kubo R. 1954. Note on the stochastic theory of resonance absorption. Journal of the Physical Society of Japan 9(6):935–944
  * [34] Siemens ME, Moody G, Li H, Bristow AD, Cundiff ST. 2010. Resonance lineshapes in two-dimensional Fourier transform spectroscopy. Optics Express 18(17):17699–17708
  * [35] Bristow AD, Zhang T, Siemens ME, Cundiff ST, Mirin R. 2011. Separating homogeneous and inhomogeneous line widths of heavy-and light-hole excitons in weakly disordered semiconductor quantum wells. J. Phys. Chem. B 115(18):5365–5371
  * [36] Li H, Piryatinski A, Srimath Kandada AR, Silva C, Bittner ER. 2019. Photon entanglement entropy as a probe of many-body correlations and fluctuations. The Journal of Chemical Physics 150(18):184106
  * [37] Neutzner S, Thouin F, Cortecchia D, Petrozza A, Silva C, Srimath Kandada AR. 2018\. Exciton-polaron spectral structures in two dimensional hybrid lead-halide perovskites. Phys. Rev. Mater. 2(6):064605
  * [38] Thouin F, Neutzner S, Cortecchia D, Dragomir VA, Soci C, et al. 2018. Stable biexcitons in two-dimensional metal-halide perovskites with strong dynamic lattice disorder. Phys. Rev. Mater. 2(3):034001
  * [39] Thouin F, Valverde-Chávez DA, Quarti C, Cortecchia D, Bargigia I, et al. 2019b. Phonon coherences reveal the polaronic character of excitons in two-dimensional lead halide perovskites. Nat. Mater. 18:349–356
  * [40] Thouin F, Srimath Kandada AR, Valverde-Chávez DA, Cortecchia D, Bargigia I, et al. 2019c. Electron-phonon couplings inherent in polarons drive exciton dynamics in two-dimensional metal-halide perovskites. Chem. Mater. 31:7085–7091
  * [41] Srimath Kandada AR, Li H, Thouin F, Bittner ER, Silva C. 2020. Stochastic scattering theory for excitation-induced dephasing: Time-dependent nonlinear coherent exciton lineshapes. The Journal of Chemical Physics 153(16):164706
  * [42] Srimath Kandada AR, Li H, Bittner ER, Silva-Acuña C. 2022. Homogeneous optical line widths in hybrid ruddlesden–popper metal halides can only be measured using nonlinear spectroscopy. The Journal of Physical Chemistry C 126(12):5378–5387
  * [43] Cho M. 2008. Coherent two-dimensional optical spectroscopy. Chem. Rev. 108(4):1331–1418
  * [44] Mukamel S. 1995. Principles of Nonlinear Optics and Spectroscopy. Oxford University Press
  * [45] Ciuti C, Savona V, Piermarocchi C, Quattropani A, Schwendimann P. 1998. Role of the exchange of carriers in elastic exciton-exciton scattering in quantum wells. Physical Review B 58(12):7926–7933
  * [46] You Y, Zhang XX, Berkelbach TC, Hybertsen MS, Reichman DR, Heinz TF. 2015. Observation of biexcitons in monolayer WSe2. Nature Physics 11:477–482
  * [47] Kylänpää I, Komsa HP. 2015. Binding energies of exciton complexes in transition metal dichalcogenide monolayers and effect of dielectric environment. Physical Review B 92(20):205418
  * [48] Bolzonello L, Fassioli F, Collini E. 2016. Correlated fluctuations and intraband dynamics of j-aggregates revealed by combination of 2des schemes. The Journal of Physical Chemistry Letters 7(24):4996–5001PMID: 27973862
  * [49] Born M. 1926. Quantenmechanik der stoßvorgänge. Zeitschrift für Physik 38(11-12):803–827
  * [50] Li H, Srimath Kandada AR, Silva C, Bittner ER. 2020. Stochastic scattering theory for excitation-induced dephasing: Comparison to the anderson–kubo lineshape. The Journal of Chemical Physics 153(15):154115
  * [51] de Gennes PG. 1999. Superconductivity Of Metals And Alloys. CRC PRess
  * [52] Marcinkiewicz J. 1939. Sur une propriété de la loi de Gauß. Mathematische Zeitschrift 44(1):612–618
  * [53] Rajagopal AK, Sudarshan ECG. 1974. Some generalizations of the marcinkiewicz theorem and its implications to certain approximation schemes in many-particle physics. Phys. Rev. A 10(5):1852–1857
  * [54] Fuller FD, Ogilvie JP. 2015. Experimental implementations of two-dimensional fourier transform electronic spectroscopy. Annu. Rev. Phys. Chem. 66:667–690
  * [55] Tokmakoff A. 2000. Two-dimensional line shapes derived from coherent third-order nonlinear spectroscopy. J. Phys. Chem. A 104(18):4247–4255
  * [56] Kubo R. 1969. A stochastic theory of line shape. John Wiley & Sons, Ltd, 101–127
  * [57] Reichman D, Silbey RJ, Suárez A. 1996. On the nonperturbative theory of pure dephasing in condensed phases at low temperatures. The Journal of Chemical Physics 105(23):10500–10506
  * [58] Hsu D, Skinner JL. 1984. On the thermal broadening of zero-phonon impurity lines in absorption and fluorescence spectra. The Journal of Chemical Physics 81(4):1604–1613
  * [59] Skinner JL, Hsu D. 1986. Pure dephasing of a two-level system. The Journal of Physical Chemistry 90(21):4931–4938
  * [60] Mukamel S. 1984. Stochastic theory of resonance Raman line shapes of polyatomic molecules in condensed phases. The Journal of Chemical Physics
  * [61] Bittner ER, Li H, Piryatinski A, Srimath Kandada AR, Silva C. 2020. Probing exciton/exciton interactions with entangled photons: Theory. The Journal of Chemical Physics 152(7):071101