# Aspects of two-photon absorption of squeezed light: the CW limit

C. Drago<EMAIL_ADDRESS>J. E. Sipe  Department of Physics,
University of Toronto, 60 St. George Street, Toronto, Ontario, Canada, M5S 1A7

###### Abstract

We present a theoretical analysis of two-photon absorption of classical and
squeezed light valid when one-photon absorption to an intermediate state is
either resonant or far-detuned from resonance, and in both the low and high
intensity regimes. In this paper we concentrate on continuous-wave excitation,
although the approach we develop is more general. We calculate the energy
removed from an incident field for typical experimental parameters and
consider the limiting cases when the photon pairs are narrowband or broadband
compared to the molecular linewidths. We find an enhancement of the two-photon
absorption due to resonant contributions from the large squeezed light
bandwidth and due to photon bunching in the low intensity regime. However, in
both cases, for the parameters we choose, the one-photon absorption is the
dominant process in the region of parameter space where a large enhancement of
the two-photon absorption is possible.

††preprint: APS/123-QED

##  I introduction

With sources of non-classical light becoming widely available, the study of
advantages such light might offer in spectroscopic and sensing methods is an
active area of research [1, 2, 3, 4, 5]. One current research topic is whether
or not the use of broadband photon pairs can lead to enhanced rates for two-
photon absorption, thereby increasing the effectiveness of applications that
rely on classical two-photon absorption, such as two-photon microscopy,
photodynamic cancer therapy, and 3D photopolymerization [6, 7, 8].

Photon pairs can be generated such that their correlation time is much less
than molecular decay times, while also being anti-correlated in frequency such
that pairs of photons have the correct total energy to lead to two-photon
absorption. It was first argued in the 1980s that correlated pairs of photons
would lead to a scaling of two-photon absorption linear with the incident
intensity [9, 10, 11]. This was shortly confirmed with experiments on trapped
atomic cesium [12], and followed by many theoretical studies of two-photon
absorption [1, 2, 3, 13, 2, 14, 15, 16] and its potential advantage in
applications.

Other experiments also followed on two-photon absorption [17, 18, 19, 20] and
sum frequency generation [21, 22], and some researchers have advocated for the
development of two-photon microscopy [23]. In each of these studies it is
claimed that two-photon absorption can be achieved with the use of photon
pairs at photon fluxes many orders of magnitude lower than would be required
were classical light used. However, these results have been called into
question by recent work, where no enhancement due to photon pairs is observed
[24, 25, 26], and an alternate explanation for the experimental data can be
provided [27, 28].

A new theoretical initiative has been led by Raymer et al. [29, 30, 31, 32],
where the low flux (“isolated pair”) regime was considered. This was later
generalized to the high flux regime [33]. In both studies they considered the
“far detuned limit,” i.e., for the incident light used it is assumed that
there are no one-photon resonances involving the ground state and intermediate
states of the molecule. Their conclusion is that in the low flux regime the
predicted rate of two-photon absorption should be undetectable with current
technologies, while in the high flux regime, where the squeezed light
bandwidth is much broader than molecular linewidths, squeezed light does not
provide any enhancement. These new theoretical results are in agreement with
the results of recent experiments [24, 25, 26, 27].

In the far detuned limit, the low flux regime can be treated by usual
perturbation approaches. For the high flux regime, there are, broadly
speaking, two strategies that can be employed. The first is to “begin at the
beginning,” with a nonlinear interaction Hamiltonian governing the generation
of broadband squeezed light generated via a narrowband pump [13]. This was the
approach of Dayan, who found that in the high flux regime there are two
quadratic dependencies on intensity, labeled as “coherent” and “incoherent”.
Alternatively, if the incident field is squeezed light generated by an
ultrashort pump pulse, one can employ a Schmidt decomposition of the squeezed
state, allowing for the calculation of the time-frequency correlation
functions governing the absorption process [15]. Here one must sum over the
Schmidt modes, but the structure of the calculation is close in form to the
structure of the calculation following from the approach of Dayan. Indeed, in
their recent generalization to the high flux regime Raymer et al. [33] used a
method closely related to Dayan’s with similar results.

What has generally not been addressed is the “resonance limit,” by which we
mean that the incident light is on or close to a one-photon resonance between
the ground state and one of the excited states of the molecule. In this
article we present a new description of two-photon absorption of squeezed
light (and classical light) valid in both the far-detuned and resonance
limits, in both the low and high intensity regimes, and for both CW and pulsed
excitation.

Our formalism differs from past treatments in that we work with Heisenberg
operators instead of employing a density operator approach [34, 35, 32].
Within our formalism we include the effects of both non-radiative decay and
dephasing in our calculations, the first of which is treated dynamically in
the Markov limit, but could be generalized to beyond that limit. Further, we
calculate the total energy removed from the incident field instead of the
final state populations. This has the advantage that it is directly related to
what is measured by experiment, although a simple distinction between one- and
two-photon absorption is not always apparent.

A central quantity in squeezed light two-photon absorption is the photon
correlation functions. Our calculation of the squeezed light correlation
functions differs from those of Dayan [13] and Raymer et al. [33] in that we
begin with a joint spectral amplitude for the incident squeezed light, and
derive a field operator transformation that is used to evaluate the time-
frequency correlation functions leading to the absorption. We use a simple
model of the joint spectral amplitude of the squeezed light to derive explicit
functions and scalings with intensity, and consider them in situations where
the photons are narrowband compared to the molecular linewidths, and in
situations where they are broadband compared to molecular linewidths. However,
more complicated models for the squeezed light could also be implemented using
our approach.

In this first communication we focus on CW excitation, and compare the energy
absorption of classical light with that of squeezed light. In agreement with
the results of earlier studies, we find the two scalings with intensity
previously labeled as “coherent” and “incoherent”. We find that when the
incident squeezed light is narrow band with respect to molecular broadening,
the squeezed light cross-sections are equivalent to those of classical light
and the ratio of absorption is given by the normalized second order
correlation function $g^{(2)}(0,0)$, in agreement with recent work [33, 36,
37, 38, 39].

We then turn to a consideration of how the classical and broadband squeezed
light cross-sections vary with the detuning from one-photon resonances for
typical experimental parameters, extending past treatments [1, 2, 3, 13, 2,
14, 15, 16]. We show that the near resonant behaviour is nontrivial, and there
is a significant enhancement of the two-photon absorption cross-section of
broadband squeezed light when the magnitude of detuning is less than half the
squeezed light bandwidth. However, for the typical parameters we adopt and in
the region of parameter space where this enhancement is visible, the one-
photon absorption cross-section is also significantly enhanced, and overpowers
the two-photon absorption.

In the far detuned limit, which for squeezed light is when the magnitude of
detuning is greater than half the squeezed light bandwidth, and when the
squeezed light bandwidth is greater than molecular linewidths, the ratio of
absorption is _approximately_ given by $g^{(2)}(0,0)$. Our results agree with
the calculation made by Raymer et al. [33] when the squeezed light bandwidth
is _much_ greater than molecular linewidths. Further, for broadband squeezed
light in the low flux regime, where $g^{(2)}(0,0)\gg 1$ we do find that the
squeezed light could provide an advantage. However, here again we find that
the one-photon absorption dominates.

The outline of the paper is as follows: Working in the rotating-wave-
approximation (RWA), in Section II we present a model of a multilevel molecule
(see Fig. 1) that allows one to explore a wide range of the parameter space.
The incident field can be either far-detuned from, or near a one-photon
resonance with, the excited states.

To properly model a molecule in typical experimental settings we include two
sources of broadening, the nonradiative decay of excited states and dephasing.
For large molecules with many degrees of freedom there are many ways in which
absorbed energy can decay. To model these processes we couple the molecule to
a quantum reservoir for each nonradiative decay pathway. For typical
experiments, molecules are in solution and thus subjected to broadening due to
molecular collisions. To include this broadening affect we include in our
Hamiltonian a stochastic fluctuation of the energy levels of the molecule.

In many previous treatments, the calculations are referred to as addressing
the “probability of two-photon absorption.” However, long after an exciting
pulse has passed the molecule it returns to its ground electronic state, and
some fraction of the energy extracted from the exciting pulse has been re-
emitted as fluorescence. In perhaps the most proper meaning of the word
“absorption,” that fraction of the incident energy has not been “absorbed,”
since it has been returned to the electromagnetic field. In Section III we
derive a quantum mechanical equation for the amount of energy removed from the
electromagnetic field, perhaps most properly the _absorption_. However, in
this first communication we treat the fluorescence phenomenologically, by
extending the use of the reservoir to include the effects of radiation from
fluorescent molecules with high quantum efficiencies that are used for current
squeezed light two-photon absorption experiments. Thus the “absorption” we
then calculate has its more colloquial meaning as the energy removed from the
incident field.

In Sections IV-VIII we develop the interaction picture in which we will be
working, derive the equations of motion, and calculate the absorption for any
incident field (pulsed or CW) near or far from resonance. Then in Section IX
we specialize to a three level model for which we can construct closed form
results and set up the calculations to follow.

Finally, with the formalism developed and general expressions for the
absorption presented, we specify to different incident fields. In Sections X
and XI we provide an analysis of CW coherent and CW squeezed light in various
regimes. Technical details are presented in appendices.

##  II model Hamiltonian

We begin with a simple model of a molecule which has a ground state $\ket{g}$,
a set of states $\\{\ket{e}\\}$ that can serve as intermediate states in two-
photon absorption, and a set of states $\\{\ket{f}\\}$ that can serve as final
states. For a review of quantum light-matter interactions, see, e.g. Gerry and
Knight [34], and Mukamel [35]. The level diagram for such a molecule is shown
in Fig. 1. The free molecule Hamiltonian, $H_{\text{M}}$, is then

$H_{\text{M}}=\hbar\omega_{g}\sigma_{gg}+\sum_{e}\hbar\omega_{e}\sigma_{ee}+\sum_{f}\hbar\omega_{f}\sigma_{ff}$
(II.1)

where $\sigma_{ij}=\ket{i}\bra{j}$. We assume these states can be approximated
as a complete basis,

$\hat{1}=\sigma_{gg}+\sum_{e}\sigma_{ee}+\sum_{f}\sigma_{ff}.$ (II.2)

From here on we do not explicitly indicate sums over intermediate $(e)$ or
final $(f)$ states, and use the convention that whenever such state indices
appear on the right-hand side of an equation, but not on the left-hand side,
they are summed over.

The free EM field Hamiltonian, $H_{\text{EM}}$, is given by

$H_{\text{EM}}=\sum_{I}\int
d\bm{k}\hbar\omega_{k}a_{I}^{\dagger}(\bm{k})a_{I}(\bm{k}),$ (II.3)

with $k=|\bm{k}|$, where

$[a_{I}(\bm{k}),a_{J}^{\dagger}(\bm{k}^{\prime})]=\delta(\bm{k}-\bm{k}^{\prime})\delta_{IJ},$
(II.4)

and where the subscript labels the polarization of the field. We treat the
coupling between the molecule and the field within the electric dipole
approximation, neglecting certain terms that lead to renormalization effects
[40], by using an interaction Hamiltonian

$H_{\text{M-EM}}=-\bm{\mu}\cdot\bm{E},$ (II.5)

where $\bm{\mu}$ is the dipole moment operator and $\bm{E}$ is the electric
field operator at the position of the molecule. These two operators commute,
so we are free to order them in any way.

We expand each operator in terms of its positive and negative frequency
components $\bm{E}=\bm{E_{+}}+\bm{E_{-}}$ and
$\bm{\mu}=\bm{\mu_{+}}+\bm{\mu_{-}}$, with
$(\bm{E_{+}})^{\dagger}=\bm{E_{-}}$, $(\bm{\mu_{+}})^{\dagger}=\bm{\mu_{-}}$.
We assume the relevant incident field components are close enough to
resonances that the rotating-wave approximation (RWA) is valid; within this
approximation the interaction is given by

$H_{\text{M-EM}}=-\bm{\mu_{-}}\cdot\bm{E_{+}}-\bm{E_{-}}\cdot\bm{\mu_{+}},$
(II.6)

where we have chosen the normally ordered form, which is the simpler to work
with. In Appendix A we confirm that the results are independent of the order
chosen.

We assume that the only dipole allowed transitions of interest are
$\ket{g}\leftrightarrow\\{\ket{e}\\}$ and
$\\{\ket{e}\\}\leftrightarrow\\{\ket{f}\\}$. The transitions are sketched in
Fig. 1, where we indicate one scenario where there is resonance with an
intermediate state and one where there is not. For atomic systems, where
parity is a good quantum number, if the transitions
$\ket{g}\leftrightarrow\\{\ket{e}\\}$ and
$\\{\ket{e}\\}\leftrightarrow\\{\ket{f}\\}$ are allowed then the transition
$\ket{g}\leftrightarrow\\{\ket{f}\\}$ would be forbidden; for large florescent
molecules the situation is more complicated. In any case, in this paper we
treat any direct transition $\\{\ket{f}\\}\leftrightarrow\ket{g}$, or indeed
any transition from $\ket{f}$ to other states that might exist with energies
close to that of $\ket{g}$ but not shown in Fig. 1, in a phenomenological way,
if they are relevant, as we discuss in Section III below. The expansion of the
dipole moment $\bm{\mu_{+}}$ using the molecular basis states is then

$\bm{\mu_{+}}=\bm{\mu}_{ge}\sigma_{ge}+\bm{\mu}_{ef}\sigma_{ef},$ (II.7)

where $\bm{\mu}_{ij}=\bra{i}\bm{\mu}\ket{j}$.

Figure 1: Level diagram for the molecular system we are considering. We make
no assumptions to where the intermediate set of levels $\\{\ket{e}\\}$ and
final levels $\\{\ket{f}\\}$ are, they can be near or far from resonance
(within the RWA). We include the resonant and non-resonant diagram. We assume
the only allowed transitions from the external field are between
$\ket{g}\leftrightarrow\\{\ket{e}\\}$ and
$\\{\ket{e}\\}\leftrightarrow\\{\ket{f}\\}$.

In many two-photon absorption experiments [18, 19, 23, 24, 25, 26, 27], large
molecules are put in solution and then irradiated. Such molecules have many
degrees of freedom that can be involved in the non-radiative decay of
optically excited electronic states, leading to a loss of energy from the
electromagnetic field. As well, the molecule suffers fluctuations in the local
environment that can lead to fluctuations in its energy levels, and thus to
dephasing of induced dipole moments. We include both effects in the following
sections.

###  II.1 Non-radiative decay

To model the first of these we couple the molecule to a quantum reservoir.
Such a reservoir is often taken to be a set of harmonic oscillators, with the
coupling to the system of interest treated in the Markov limit; this leads to
a Lindblad equation for the reduced density operator of the system [41].
However, since here we want to keep the details of correlations between the
molecule and the electromagnetic field in our equations, only the reservoir
would be “traced over,” and the reduced density operator of the system of
interest – molecule and electromagnetic field – would act on the combined
Hilbert space of the molecule _and_ the electromagnetic field. To avoid such a
complicated entity we work instead in the Heisenberg picture, and use the
result of Fischer [42] that dynamics at the Markov level can be modeled by
coupling each transition of interest to the excitations in a formal 1D
waveguide, which is taken as a reservoir at zero temperature. In such a model,
quantum fluctuations in the waveguide propagate towards the molecule before
interacting with it, and then away from it after the interaction, assuring
that the molecule is effectively always interacting with the vacuum state of
the quantum reservoir. This idealizes the physical assumption that the
coherence generated between the electronic transitions of the molecule and the
other degrees of freedom, which is responsible for the loss of energy from the
electromagnetic field, quickly decays.

In our implementation of this approach, we assume there are independent non-
radiative decay pathways of three types, involving transitions between
$\ket{g}\leftrightarrow\\{\ket{e}\\}$,
$\\{\ket{e}\\}\leftrightarrow\\{\ket{f}\\}$ and
$\ket{g}\leftrightarrow\\{\ket{f}\\}$ respectively. We use the form for a 1D
waveguide Hamiltonian introduced earlier [43], where the propagation of light
in integrated photonic structures was discussed. For each decay process we
introduce the waveguide “reservoir operators” $\psi_{ij}(z)$ and
$\psi^{\dagger}_{ij}(z)$, which will lead to transitions between $\ket{i}$ and
$\ket{j}$, with $\omega_{i}>\omega_{j}$. For a single waveguide, the free
reservoir Hamiltonian is given by

$\begin{split}H_{\text{R}}^{ij}=\hbar\omega_{ij}&\int
dz\psi^{\dagger}_{ij}(z)\psi_{ij}(z)\\\ &+\frac{i\hbar v_{ij}}{2}\int
dz\left(\frac{d\psi^{\dagger}_{ij}(z)}{dz}\psi_{ij}(z)-\text{H.c.}\right),\end{split}$
(II.8)

where $\omega_{ij}=\omega_{i}-\omega_{j}$ is taken as a center frequency and
$v_{ij}$ is the group velocity of excitation propagation in the waveguide. We
take the waveguide operators to satisfy the usual commutation relation

$\displaystyle[\psi_{ij}(z),\psi_{mn}^{\dagger}(z^{\prime})]=\delta_{im}\delta_{jn}\delta(z-z^{\prime}),$
(II.9)

which ensures the dynamics of all the reservoir operators are independent. The
full reservoir Hamiltonian is then written as the sum of three contributions,

$H_{\text{R}}=H_{\text{R}}^{eg}+H_{\text{R}}^{fe}+H_{\text{R}}^{fg},$ (II.10)

where $H_{\text{R}}^{eg}$, $H_{\text{R}}^{fe}$ and $H_{\text{R}}^{fg}$ contain
reservoir operators that facilitate transitions between
$\ket{g}\leftrightarrow\\{\ket{e}\\}$,
$\\{\ket{e}\\}\leftrightarrow\\{\ket{f}\\}$ and
$\ket{g}\leftrightarrow\\{\ket{f}\\}$, respectively. Again, and below, there
is an implicit sum over the intermediate and final states.

For the interaction between a transition $\sigma_{ij}$ and its associated
reservoir, we assume a point coupling between the molecule and the waveguide
modeling the reservoir. After applying the RWA, the coupling responsible for
transitions between $\ket{i}\leftrightarrow\ket{j}$ due to the reservoir is
given by [43]

$H_{\text{M-R}}^{ij}=\hbar\eta_{ij}^{*}\sigma_{ij}\psi_{ij}(0)+\hbar\eta_{ij}\psi^{\dagger}_{ij}(0)\sigma_{ji},$
(II.11)

taking the “position” of the molecule to be the origin of the fictitious
waveguide, with the full coupling Hamiltonian then given by

$H_{\text{M-R}}=H_{\text{M-R}}^{eg}+H_{\text{M-R}}^{fe}+H_{\text{M-R}}^{fg}.$
(II.12)

In equation (II.8) we introduced a nominal reference frequency and group
velocity, and in Eq. (II.11) we introduced a nominal coupling constant
$\eta_{ij}$. In fact, the group velocity $v_{ij}$ and coupling $\eta_{ij}$
will lead to a new constant that will identify the non-radiative decay widths
associated with each transition $\ket{i}\leftrightarrow\ket{j}$; it will be
the physical parameter that can be identified by experiment, and there will be
one of these physical parameters for each coupling of a transition to a
waveguide reservoir.

###  II.2 Dephasing

We now turn to the effects of fluctuations in the energy levels of the
molecule that arise due to fluctuations in its environment. The dephasing that
results could be captured by a density operator formalism with a quantum
reservoir model such as that described above, but one that would involve terms
$H^{ij}_{\text{M-R}}$ with $i=j$ [44]. In the low temperature regime such a
treatment is necessary, but one can consider a simpler model of dephasing that
reproduces the fully quantum treatment at higher temperatures [44, 45, 46,
35]. In that simpler model the molecule is subject to classical stochastic
fluctuations in its energy eigenvalues, described with a contribution to the
Hamiltonian given by

$H_{\text{fluc}}(t)=\hbar\tilde{\omega}_{g}(t)\sigma_{gg}+\hbar\tilde{\omega}_{e}(t)\sigma_{ee}+\hbar\tilde{\omega}_{f}(t)\sigma_{ff},$
(II.13)

where the $\tilde{\omega}_{i}(t)$ are classical stochastic functions of time.
We make four assumptions concerning these random variables [44, 45]: i) They
have zero mean; ii) The correlation between state $\ket{i}$ and $\ket{j}$ is
characterized by $c_{ij}$; iii) the correlation between two variables at
different times follows a fast exponential decay with a time scale $\tau_{c}$;
and iv) the random variables follow a Gaussian distribution. We denote the
classical expectation value as $[\\![\cdot]\\!]$, and then the first three
assumptions are summarized by

$\displaystyle[\\![\tilde{\omega}_{i}(t)]\\!]=0$ (II.14a)
$\displaystyle[\\![\tilde{\omega}_{i}(t)\tilde{\omega}_{j}(t^{\prime})]\\!]=c_{ij}e^{-|t-t^{\prime}|/\tau_{c}}.$
(II.14b)

In the full problem we calculate all quantities of interest, and then finally
take the stochastic averages as given above to describe the effects of
dephasing.

Then the full Hamiltonian for our model is

$H(t)=H_{0}+V(t),$ (II.15)

where $H_{0}$ is the free Hamiltonian given by

$H_{0}=H_{\text{M}}+H_{\text{EM}}+H_{\text{R}},$ (II.16)

and $V(t)$ the time dependent interaction term given by

$V(t)=H_{\text{M-EM}}+H_{\text{M-R}}+H_{\text{fluc}}(t).$ (II.17)

We note that although $V(t)=V^{\dagger}(t)$ is time dependent it is still a
Schrödinger operator.

To write the interaction Hamiltonian in a compact form we define the operators

$F_{+}^{ij}=\frac{\bm{\mu}_{ij}}{\hbar}\cdot\bm{E_{+}}-\eta_{ij}^{*}\psi_{ij}(0),$
(II.18)

and their adjoints,

$(F_{+}^{ij})^{\dagger}=\frac{\bm{\mu}_{ji}}{\hbar}\cdot\bm{E_{-}}-\eta_{ij}\psi_{ij}^{\dagger}(0)\equiv
F_{-}^{ji}.$ (II.19)

The interaction $V(t)$ can then be written as

$\begin{split}V(t)=&-\hbar\sigma_{eg}F_{+}^{eg}-\hbar\sigma_{fe}F_{+}^{fe}-\hbar\sigma_{fg}F_{+}^{fg}+\text{H.c.}\\\
&+H_{\text{fluc}}(t).\end{split}$ (II.20)

We note that there are three $F_{+}^{ij}$ operators that contribute to the
interaction: $F_{+}^{eg}$, $F_{+}^{fe}$, and $F_{+}^{fg}$. However, the last
term $F_{+}^{fg}$ is _sui_ _generis_. It only contains reservoir operators, as
we have assumed there are no one-photon transition between
$\ket{g}\leftrightarrow\\{\ket{f}\\}$, i.e., $\bm{\mu}_{fg}=0$. The form of
the first line in Eq. (II.20) allows a simple interpretation: the energy of
the molecule is increased (decreased) by annihilating (creating) energy in the
external field or reservoir.

##  III absorption(s)

The term _absorption_ is used in a number of different ways in optics. In one
usage it refers to the removal of energy from the electromagnetic field, and
in our discussion here we begin with that meaning. In Appendix A we show that,
within the RWA, the absorption of light between $t_{I}$ and $t_{F}>t_{I}$ by a
molecule is given by

$\mathcal{A}=\int\limits_{t_{I}}^{t_{F}}dt\bra{\Psi(t_{0})}\bm{E_{-}}^{H}(t)\cdot\frac{d\bm{\mu_{+}}^{H}(t)}{dt}\ket{\Psi(t_{0})}+\text{c.c.},$
(III.1)

where $\ket{\Psi(t_{0})}$ is the initial ket ($t_{0}<t_{I}$) of the full
molecule-field-reservoir system, and the superscript $H$ indicates that the
operators evolve according to the full Hamiltonian. (For a discussion of the
special case where the state of the electromagnetic field is a single photon,
see Valente et al. [47]). Although the result (III.1) is exact we have yet to
solve for the Heisenberg operators $\bm{\mu_{\pm}}^{H}(t)$ and
$\bm{E_{\pm}}^{H}(t)$.

Besides renormalization contributions, which we neglect, the solution for
these quantities also involve radiation reaction effects [48]. These we also
neglect in this paper. For typical one-photon absorption processes in
molecules this is well-justified, for the natural linewidth is much less than
that due to non-radiative decay and dephasing. But for typical fluorescent
molecules used for two-photon absorption experiments [17, 18, 19], there can
be large fluorescence quantum efficiencies that would make this neglect
suspect. To compensate for this we expand the use the reservoir introduced in
Section II to also provide a phenomenological model of the decay of the states
$\\{\ket{f}\\}$ to the ground state $\ket{g}$ due to fluorescent processes; we
leave a more complete description with radiation reaction to later work

Then the expression for the absorption is given by

$\mathcal{A}\rightarrow\int\limits_{t_{I}}^{t_{F}}dt\bra{\Psi(t_{0})}\bm{E_{-}}^{H_{0}}(t)\cdot\frac{d\bm{{\mu}_{+}}^{H^{\prime}}(t)}{dt}\ket{\Psi(t_{0})}+c.c.,$
(III.2)

where by $\bm{\mu_{\pm}}^{H^{\prime}}(t)$ we mean the evolution of those
quantities according to the Heisenberg equations, but with
$\bm{E_{\pm}}^{H}(t)$ in those equations replaced by
$\bm{E_{\pm}}^{H_{0}}(t)$. Since this expression for the “absorption” only
_models_ the fluorescent contributions to the extinction as loss from the full
electromagnetic field, it includes effects that physically are not absorptive
in the strict sense introduced at the start of this section. Nonetheless, in
the more colloquial sense of “absorption” implicit in the phrase “two-photon
absorption,” this can be understood as describing the removal of energy from
the components of the electromagnetic field near the fundamental frequency.
For the rest of this paper we use Eq. (III.2) to calculate _this_
“absorption,” but to avoid clutter in the notation we henceforth drop the
primes on the $H^{\prime}$ in $\bm{\mu_{\pm}}^{H^{\prime}}(t)$, and also drop
the inverted commas around “absorption.”

We then proceed by using the expression for $\bm{\mu_{\pm}}^{H}(t)$. The
calculated absorption then has two contributions, the first due to the
coherence operator $\sigma_{ge}^{H}(t)$ and the second is due to
$\sigma_{ef}^{H}(t)$. Since within our approximations Eq. (III.2) identifies
the _full_ absorption, we will see that a simple division of the result into
“one-photon” and “two-photon” terms is not always possible, especially when
the fundamental frequency is on resonance with a transition to an intermediate
level.

We proceed by expanding each coherence operator $\sigma_{ij}^{H}(t)$ in a
perturbative series, using the interaction picture we present below.

##  IV Interaction Picture

In our full Hamiltonian (II.15), $V(t)$ can be written as a function of a set
of Schrödinger operators $\left\\{O_{\alpha}\right\\}$ and time,
$V\left(\left\\{O_{\alpha}\right\\};t\right)$. The time evolution of a ket
from $t_{a}$ to $t_{b}$ in the absence of
$V\left(\left\\{O_{\alpha}\right\\};t\right)$ is formally described by the
unitary evolution operator $\mathcal{U}_{0}(t_{b},t_{a})$, which satisfies the
Schrödinger equation

$i\hbar\frac{d}{dt_{b}}\mathcal{U}_{0}(t_{b},t_{a})=H_{0}\mathcal{U}_{0}(t_{b},t_{a}),$
(IV.1)

and has a solution

$\mathcal{U}_{0}(t_{b},t_{a})=e^{-iH_{0}(t_{b}-t_{a})/\hbar}\,.$ (IV.2)

Including $V\left(\left\\{O_{\alpha}\right\\};t\right)$, the full evolution
operator $\mathcal{U}(t_{b},t_{a})$ satisfies the equation

$i\hbar\frac{d}{dt_{b}}\mathcal{U}(t_{b},t_{a})=H(t_{b})\mathcal{U}(t_{b},t_{a}),$
(IV.3)

with the initial condition $\mathcal{U}(t_{a},t_{a})=\hat{1}$ for all $t_{a}$.

Now consider times $t_{\text{min}}$ and $t_{\text{max}}$ such that for $t\leq
t_{\text{min}}$ or $t\geq t_{\text{max}}$ the interaction between the
electromagnetic field and the molecule vanishes. In general, of course, such
times do not exist; even when there is no pulse of light, if initially a
molecule were in its ground state or any other eigenstate of $H_{\text{M}}$
(Eq. (II.1)), it would interact with the quantized electromagnetic field.
However, for the RWA form of the Hamiltonian we adopt, such “vacuum state”
interactions vanish; see Appendix A. So for times $t_{a},t_{b}\leq
t_{\text{min}}$ or $t_{a},t_{b}\geq t_{\text{max}}$ we have

$\mathcal{U}(t_{a},t_{b})=\mathcal{U}_{0}(t_{a},t_{b}).$ (IV.4)

It is then useful to introduce the operator

$U(t_{b},t_{a})=\mathcal{U}_{0}(0,t_{b})\mathcal{U}(t_{b},t_{a})\mathcal{U}_{0}(t_{a},0),$
(IV.5)

which satisfies all the properties of a unitary time evolution operator,
including the condition that $U(t_{a},t_{a})=\hat{1}$ for all $t_{a}$. We can
then write

$\left|\Psi(t)\right\rangle=\mathcal{U}_{0}(t,0)U(t,t_{0})\left|\Psi_{\text{in}}\right\rangle,$
(IV.6)

where the “asymptotic-in” ket

$\left|\Psi_{\text{in}}\right\rangle\equiv\mathcal{U}_{0}(0,t_{0})\left|\Psi(t_{0})\right\rangle$
(IV.7)

is what the ket would be at $t=0$ were there no interaction $V(t)$ and
$t_{0}<t_{I},t_{\text{min}}$ is the initial start time. Note that since
$t_{0}$ is long before any interaction occurs,
$\left|\Psi_{\text{in}}\right\rangle$ is independent of $t_{0}$. From the
definition (IV.5) we can show, using (IV.4), that for
$t_{a},t_{a}^{\prime}\leq t_{\text{min}}$ we have
$U(t,t_{a})=U(t,t_{a}^{\prime}$), and thus both must be equal to
$U(t,-\infty)$; if we put

$U(t)\equiv U(t,-\infty)$ (IV.8)

we can write (IV.6) as

$\left|\Psi(t)\right\rangle=\mathcal{U}_{0}(t,0)U(t)\left|\Psi_{\text{in}}\right\rangle.$
(IV.9)

We now move to the Heisenberg picture, and develop an interaction
representation related to this.

To begin we introduce an interaction representation of each Schrödinger
operator $O_{\alpha}$,

$\hat{O}_{\alpha}(t)\equiv\mathcal{U}_{0}(0,t)O_{\alpha}\mathcal{U}_{0}(t,0).$
(IV.10)

For our operators of interest we will have, schematically,

$\hat{O}_{\alpha}(t)=\sum_{\beta}g_{\alpha\beta}(t)O_{\beta},$ (IV.11)

where the $g_{\alpha\beta}(t)$ are (classical) functions of time. For example,
if $O_{\alpha}$ is $\sigma_{ij}$ or $a(\bm{k})$ then there will be only one
such $g_{\alpha\beta}(t)$ and we have

$\begin{gathered}\hat{\sigma}_{ij}(t)=\sigma_{ij}e^{-i\omega_{ji}t},\\\
\hat{a}_{I}(\bm{k},t)=a_{I}(\bm{k})e^{-i\omega_{k}t},\end{gathered}$ (IV.12)

but if $O_{\alpha}$ is $\psi_{ij}(z),$ we have

$\begin{split}\hat{\psi}_{ij}(z,t)&=\psi_{ij}(z-v_{ij}t)e^{-i\omega_{ij}t}\\\
&=\int g_{zz^{\prime}}(t)\psi_{ij}(z^{\prime})dz^{\prime},\end{split}$ (IV.13)

where $g_{zz^{\prime}}(t)=\delta(z^{\prime}-z+v_{ij}t)\exp(-i\omega_{ij}t)$.

With the expressions (IV.10,IV.11), and the definitions of
$\mathcal{U}_{0}(t_{a},t_{b})$, $\mathcal{U}(t_{a},t_{b})$, and $U(t)$, it is
easy to confirm that the expectation value of any operator $O_{\alpha}$ is

$\bra{\Psi(t)}O_{\alpha}\ket{\Psi(t)}=\bra{\Psi_{\text{in}}}\bar{O}_{\alpha}(t)\ket{\Psi_{\text{in}}},$
(IV.14)

where

$\begin{split}\bar{O}_{\alpha}(t)&=U^{\dagger}(t)\mathcal{U}_{0}(0,t)O_{\alpha}\mathcal{U}_{0}(t,0)U(t)\\\
&=\sum_{\beta}g_{\alpha\beta}(t)\check{O}_{\beta}(t),\end{split}$ (IV.15)

and

$\check{O}_{\beta}(t)\equiv U^{\dagger}(t)O_{\beta}U(t).$ (IV.16)

Since absorption only occurs for times $t$ between the “interaction region”
$t_{\text{min}}\leq t\leq t_{\text{max}}$, as long as
$t_{0}<t_{I}<t_{\text{min}}$ and $t_{F}>t_{\text{max}}$ we can extend the
limits of integration to infinity because the absorption is zero for those
times and write Eq. (III.2) as

$\mathcal{A}=\int\limits_{-\infty}^{\infty}dt\langle\Psi_{\text{in}}|\bm{\hat{E}}_{-}(t)\cdot\frac{d\bm{\bar{\mu}}_{+}(t)}{dt}|\Psi_{\text{in}}\rangle+c.c.,$
(IV.17)

where
$\bm{\hat{E}}_{-}(t)\equiv\bm{E}_{-}(\\{\hat{a}_{I}^{\dagger}(\bm{k},t)\\})=\bm{E}_{-}(\\{a_{I}^{\dagger}(\bm{k})e^{i\omega_{k}t}\\})$,
and since in calculating $\bm{\mu}_{+}^{H}(t)$ we are to replace
$\bm{E}_{+}^{H}(t)$ by $\bm{E}_{+}^{H_{0}}(t)$, in calculating
$\bar{\bm{\mu}}_{+}(t)$ we are to replace
$\bm{\bar{E}}_{+}(t)\equiv\bm{E}_{+}(\\{\bar{a}_{I}(\bm{k},t)\\})=\bm{E}_{+}(\\{\check{a}_{I}(\bm{k},t)e^{-i\omega_{k}t}\\})$
by
$\bm{\hat{E}}_{+}(t)\equiv\bm{E}_{+}(\\{\hat{a}_{I}(\bm{k},t)\\})=\bm{E}_{+}(\\{a_{I}(\bm{k})e^{-i\omega_{k}t}\\})$.
This we will do in the equations of motion we derive below.

##  V equations of motion

The operator $\bm{\bar{\mu}}_{+}(t)$ involves the operators
$\bar{\sigma}_{ge}(t)$ and $\bar{\sigma}_{ef}(t)$, and we will see below that
the dynamics for those operators involves other operators
$\bar{\sigma}_{ij}$(t), and so the $\bar{\sigma}_{ij}(t)$ will be our
operators $\bar{O}_{\gamma}(t)$ of interest. Now the complicated part of the
dynamics of an operator $\bar{O}_{\gamma}(t)$ is due to its dependence on the
$\check{O}_{\gamma}(t)$, and the evolution of an operator
$\check{O}_{\gamma}(t)$ follows from the evolution of $U(t).$ From its
definition (IV.8) the dynamical equation for $U(t)$ can be found,

$i\hbar\frac{d}{dt}U(t)=V\left(\left\\{\hat{O}_{\alpha}(t)\right\\};t\right)U(t),$
(V.1)

from which follows the dynamical equation for $\check{O}_{\gamma}(t),$

$i\hbar\frac{d}{dt}\check{O}_{\gamma}(t)=\left[\check{O}_{\gamma}(t),V\left(\left\\{\bar{O}_{\alpha}(t)\right\\};t\right)\right].$
(V.2)

In these equations the $\bar{O}_{\alpha}(t)$ are functions of the
$\check{O}_{\beta}(t)$, as given by (IV.15). Noting this relation is the same
as that between the $\hat{O}_{\alpha}(t)$ operators and the Schrödinger
operators $O_{\beta}$ (IV.11), we have

$\begin{gathered}\bar{\sigma}_{ij}(t)=\check{\sigma}_{ij}(t)e^{-i\omega_{ji}t},\\\
\bar{a}_{I}(\bm{k},t)=\check{a}_{I}(\bm{k},t)e^{-i\omega_{k}t},\\\
\bar{\psi}_{ij}(z,t)=\check{\psi}_{ij}(z-v_{ij}t,t)e^{-i\omega_{ij}t}.\end{gathered}$
(V.3)

Since the relation between $\bar{\sigma}_{ij}(t)$ and $\check{\sigma}_{ij}(t)$
is so simple, and it is the operators $\bar{O}_{\alpha}(t)$ appearing in
$V\left(\left\\{\bar{O}_{\alpha}(t)\right\\};t\right)$ in (V.2), we begin by
constructing dynamical equations for the $\bar{\sigma}_{ij}(t)$ using
(V.2,V.3),

$i\hbar\left(\frac{d}{dt}+i\omega_{ji}\right)\bar{\sigma}_{ij}(t)=\left[\bar{\sigma}_{ij}(t),V\left(\left\\{\bar{O}_{\alpha}(t)\right\\};t\right)\right].$
(V.4)

They involve the quantities

$\begin{gathered}\bar{F}_{+}^{ij}(t)\equiv\frac{\bm{\mu}_{ij}}{\hbar}\cdot\bm{E}_{+}(\\{\bar{a}_{I}(\bm{k},t)\\})-\eta_{ij}^{*}\bar{\psi}_{ij}(0,t),\end{gathered}$
(V.5)

for $\omega_{i}>\omega_{j}$, and their adjoints (see Eq. (II.18)). Note that
the dynamical equations (V.4) for the $\bar{\sigma}_{ij}(t)$ take the same
structure as the Heisenberg equations $\sigma_{ij}^{H}(t)$. However, while the
ket relevant for expectation values of the $\sigma_{ij}^{H}(t)$ is
$\left|\Psi(t_{0})\right\rangle$, that relevant for expectation values of the
$\bar{\sigma}_{ij}(t)$ is $\left|\Psi_{\text{in}}\right\rangle$. And while the
$\sigma_{ij}^{H}(t)$ equal the corresponding Schrödinger operators
$\sigma_{ij}$ at $t_{0}$, it is the operators $\check{\sigma}_{ij}(t)$ that
equal the corresponding Schrödinger operators $\sigma_{ij}$ at a special time
given by $t=-\infty$. We return to this point below.

Since the relation between $\bar{\psi}_{ij}(z,t)$ and $\check{\psi}_{ij}(z,t)$
is not as simple as that between $\bar{\sigma}_{ij}(t)$ and
$\check{\sigma}_{ij}(t)$ (see (V.3)), we simply construct the dynamical
equation for $\check{\psi}_{ij}(z,t)$ directly from (V.2); again, since an
$\check{O}_{\alpha}(t)$ is related to the corresponding Schrödinger operator
$O_{\alpha}$ by a unitary transformation (IV.16), the commutation relations
are preserved, and we find

$\frac{\partial}{\partial
t}\check{\psi}_{ij}(z,t)=-i\eta_{ij}\check{\sigma}_{ji}(t)\delta(z+v_{ij}t).$
(V.6)

The equations (V.4,V.6) are the set of dynamical equations we wish to solve;
note that (V.4) could easily be rewritten as equations for the
$\check{\sigma}_{ij}(t)$, and we would be solving for the set
$\left\\{\check{O}_{\alpha}(t)\right\\}$ of operators.

##  VI Formal solutions

The equation (V.6) that describes how the $\check{\sigma}_{ji}(t)$ affect the
$\check{\psi}_{ij}(z,t)$ can be solved immediately. Integrating from
$t=-\infty$ to $t$, we recall that
$\check{\psi}_{ij}(z,-\infty)=\psi_{ij}(z),$ the Schrödinger operator, and we
find

$\check{\psi}_{ij}(-v_{ij}t,t)=\psi_{ij}(-v_{ij}t)-\frac{i\eta_{ij}}{2v_{ij}}\check{\sigma}_{ji}(t).$
(VI.1)

Multiplying by $\exp(-i\omega_{ij}t)$, using (V.3) we can also write this as

$\begin{split}&\bar{\psi}_{ij}(0,t)=\hat{\psi}_{ij}(0,t)-\frac{i\eta_{ij}}{2v_{ij}}\bar{\sigma}_{ji}(t),\end{split}$
(VI.2)

where we have also used (IV.13). Note that $\hat{\psi}_{ij}(0,t)$ can be
written in terms of Schrödinger operators, and as such can be thought of as
the “input” of the reservoir field at the site of the molecule;
$\bar{\psi}_{ij}(0,t)$ can then be thought of as the “output,” which contains
a contribution from the molecule.

With this in hand we can simplify the expression (V.5) for the
$\bar{F}^{ij}_{+}(t)$. For $\omega_{i}>\omega_{j}$ we have

$\begin{split}\bar{F}_{+}^{ij}(t)=&\frac{\bm{\mu}_{ij}}{\hbar}\cdot\bm{E}_{+}[\\{\check{a}_{I}(\bm{k},t)e^{-i\omega_{k}t}\\}]-\eta_{ij}^{*}\hat{\psi}_{ij}(0,t)\\\
&+i\frac{\left|\eta_{ij}\right|^{2}}{2v_{ij}}\bar{\sigma}_{ji}(t)\end{split}$
(VI.3)

Neglecting the radiation reaction (and renormalization) effects on the
electric field operator, as discussed in section III, we put
$\check{a}_{I}(\bm{k},t)\rightarrow a_{I}(\bm{k})$ (see discussion after (Eq.
(III.1)), and so

$\bm{E}_{+}[\\{\check{a}_{I}(\bm{k},t)e^{-i\omega_{k}t}\\}]\to\bm{E}_{+}[\\{a_{I}(\bm{k})e^{-i\omega_{k}t}\\}]=\bm{\hat{E}}_{+}(t),$
(VI.4)

and we can write

$\bar{F}_{+}^{ij}(t)=\hat{F}_{+}^{ij}(t)+i\Gamma_{ij}\bar{\sigma}_{ji}(t),$
(VI.5)

where we have put

$\hat{F}_{+}^{ij}(t)\equiv\frac{\bm{\mu}_{ij}}{\hbar}\cdot\bm{\hat{E}}_{+}\left(t\right)-\eta_{ij}^{*}\hat{\psi}_{ij}(0,t),$
(VI.6)

and

$\Gamma_{ij}=\frac{\left|\eta_{ij}\right|^{2}}{2v_{ij}}.$ (VI.7)

The reservoir parameters that contribute to calculation are now combined into
the constant $\Gamma_{ij}$, and we have written $\bar{F}_{+}^{ij}(t)$ as the
sum of two contributions. The first is from the operator
$\hat{F}_{+}^{ij}(t)$, which is a Schrödinger operator, and therefore a
“known” quantity which drives transitions. The second contribution depends on
the molecule operator $\bar{\sigma}_{ji}(t)$ which is yet to be determined.
Finally, the constant $\Gamma_{ij}$ which has units of inverse time, will turn
out to be the decay rate associated with transitions between states
$\ket{i}\leftrightarrow\ket{j}$.

Unlike the subscripts on $\sigma_{ij}$ and $\bm{\mu}_{ij}$, the subscript of
$\Gamma_{ij}$ follows from the label we gave to each reservoir operator
$\psi_{ij}(z)$, which was only defined when $\omega_{i}>\omega_{j}$; thus
$\Gamma_{ij}$ is also only defined then. However, in the following it will be
convenient to define $\Gamma_{ji}$ when $\omega_{i}>\omega_{j}$ to be
$\Gamma_{ji}\equiv\Gamma_{ij}$, i.e. the decay constant $\Gamma_{ij}$ is
defined to be symmetric. Although $\Gamma_{ij}$ is defined to be symmetric, it
always represents the decay rate of energy of the molecule transitioning from
a state of high energy to one of lower energy. The definition is convenient in
introducing the Green function below.

Using the solution (VI.5) for each $\bar{F}_{+}^{ij}(t)$ in the equation (V.4)
for the corresponding $\bar{\sigma}_{ij}(t)$, damping terms $\Gamma_{ij}$ are
introduced which we group on the left-hand side of each equation and define
the total relaxation constants due to the reservoir for each molecule operator
by

$\displaystyle\bar{\Gamma}_{eg}\equiv\Gamma_{eg},$ (VI.8a)
$\displaystyle\bar{\Gamma}_{gg}\equiv 2\Gamma_{eg},$ (VI.8b)
$\displaystyle\bar{\Gamma}_{e^{\prime}e}\equiv\Gamma_{e^{\prime}g}+\Gamma_{eg},$
(VI.8c)
$\displaystyle\bar{\Gamma}_{fe}\equiv\Gamma_{eg}+\sum_{e^{\prime}}\Gamma_{fe^{\prime}}+\Gamma_{fg},$
(VI.8d)
$\displaystyle\bar{\Gamma}_{fg}\equiv\sum_{e^{\prime}}\Gamma_{fe^{\prime}}+\Gamma_{fg},$
(VI.8e)

where each barred decay constant $\bar{\Gamma}_{ij}$ ($=\bar{\Gamma}_{ji}$) is
associated with $\bar{\sigma}_{ji}(t)$. We find

$\displaystyle\left(\frac{d}{dt}{+}\bar{\Gamma}_{eg}{+}i\omega_{eg}{+}i\tilde{\omega}_{eg}(t)\right)\bar{\sigma}_{ge}(t)=i[\bar{\sigma}_{gg}(t)\hat{F}_{+}^{eg}(t){-}\bar{\sigma}_{e^{\prime}e}(t)\hat{F}_{+}^{e^{\prime}g}(t){+}\hat{F}_{-}^{ef}(t)\bar{\sigma}_{gf}(t)],$
(VI.9a)
$\displaystyle\left(\frac{d}{dt}{+}\bar{\Gamma}_{fe}{+}i\omega_{fe}{+}i\tilde{\omega}_{fe}(t)\right)\bar{\sigma}_{ef}(t)=i[\bar{\sigma}_{ee^{\prime}}(t)\hat{F}_{+}^{fe^{\prime}}(t){-}\hat{F}_{-}^{ge}(t)\bar{\sigma}_{gf}(t){-}\bar{\sigma}_{f^{\prime}f}(t)\hat{F}_{+}^{f^{\prime}e}(t)],$
(VI.9b)
$\displaystyle\left(\frac{d}{dt}{+}\bar{\Gamma}_{fg}{+}i\omega_{fg}{+}i\tilde{\omega}_{fg}(t)\right)\bar{\sigma}_{gf}(t)=i[\bar{\sigma}_{ge}(t)\hat{F}_{+}^{fe}(t){-}\bar{\sigma}_{ef}(t)\hat{F}_{+}^{eg}(t)],$
(VI.9c)
$\displaystyle\left(\frac{d}{dt}{+}\bar{\Gamma}_{e^{\prime}e}{+}i\omega_{e^{\prime}e}{+}i\tilde{\omega}_{e^{\prime}e}(t)\right)\bar{\sigma}_{ee^{\prime}}(t)=i[\bar{\sigma}_{eg}(t)\hat{F}_{+}^{e^{\prime}g}(t){-}\hat{F}_{-}^{ge}(t)\bar{\sigma}_{ge^{\prime}}(t){-}\bar{\sigma}_{fe^{\prime}}(t)\hat{F}_{+}^{fe}(t){+}\hat{F}_{-}^{e^{\prime}f}(t)\bar{\sigma}_{ef}(t){-}2i\Gamma_{fe}\bar{\sigma}_{ff}(t)\delta_{ee^{\prime}}],$
(VI.9d)

where $\tilde{\omega}_{ij}(t)=\tilde{\omega}_{i}(t)-\tilde{\omega}_{j}(t)$.
Further, we have dropped terms involving $\hat{F}_{+}^{fg}(t)$, which can be
done with impunity: Since the equations here are normally ordered, when we
solve for each molecule operator perturbatively they will remain normally
ordered; $\hat{F}_{+}^{fg}(t)$ is a known quantity that only involves
reservoir operators, so it will always annihilate the ket
$\ket{\Psi_{\text{in}}}$, since we take the reservoir to be initially the
vacuum state. Thus the dropped terms will lead to no contribution to the
absorption. For the same reason, at this point we could replace each
$\hat{F}_{+}^{ij}(t)\to\bm{\mu}_{ij}/\hbar\cdot\bm{E}_{+}(t)$ for $ij\neq fg$,
but for now we choose to stay in the $\hat{F}_{+}^{ij}(t)$ notation to keep
the equations more condensed. We note that Eq. (VI.9), with the replacement of
$\hat{F}_{+}^{ij}(t)\to\bm{\mu}_{ij}/\hbar\cdot\bm{E}_{+}(t)$, is the usual
Heisenberg time evolution for the coherence operators; see, for example, Gerry
and Knight [34].

The equations in (VI.9) and their adjoints, together with the corresponding
equations for $\bar{\sigma}_{gg}(t)$ and $\bar{\sigma}_{f^{\prime}f}(t)$ – not
written down because they are not necessary to calculate the absorption in a
perturbative calculation – form a closed set of coupled differential equations
for the molecule operators $\bar{\sigma}_{ij}(t)$, since each
$\hat{F}_{\pm}^{ij}(t)$ is a Schrödinger operator which is known.

We now consider the general form of each equation in (VI.9), which is given by

$\left(\frac{d}{dt}+\bar{\Gamma}_{ij}+i\omega_{ij}+i\tilde{\omega}_{ij}(t)\right)\bar{\sigma}_{ji}(t)=K_{ij}(t),$
(VI.10)

or

$\left(\frac{d}{dt}+\bar{\Gamma}_{ij}+i\tilde{\omega}_{ij}(t)\right)\check{\sigma}_{ji}(t)=e^{i\omega_{ij}t}K_{ij}(t).$
(VI.11)

The second form has a formal solution that includes a homogeneous solution
satisfying the initial condition at $t=-\infty$, where each
$\check{\sigma}_{ji}(t)$ is equal to the corresponding Schrödinger operator
$\sigma_{ji}$. Because of the damping term $\bar{\Gamma}_{ij}$, that
homogeneous solution will vanish for finite times. Correspondingly, in a
formal solution of (VI.10) the homogeneous solution will vanish as well, and
introducing a Green function

$G_{ij}(t,t_{1})\equiv
ie^{-(\bar{\Gamma}_{ij}+i\omega_{ij})(t-t_{1})}e^{-i\int\limits_{t_{1}}^{t}dt^{\prime}\tilde{\omega}_{ij}(t^{\prime})},$
(VI.12)

the exact solution to Eq. (VI.10) and each operator in Eq. (VI.9) is

$\bar{\sigma}_{ji}(t)=\int\limits_{-\infty}^{t}dt_{1}G_{ij}(t,t_{1})K_{ij}(t_{1}),$
(VI.13)

for the appropriate $K_{ij}(t)$ on the right-hand side. We note that the Green
function for each equation satisfies the property that

$G_{ij}^{*}(t,t_{1})=-G_{ji}(t,t_{1}),$ (VI.14)

where we used the property that $\bar{\Gamma}_{ij}$ is symmetric.

To process the terms further we need to evaluate the time integral in the
definition of $G_{ij}(t,t_{1})$. The time integral can be evaluated by taking
the classical average of each coherence operator

$[\\![\bar{\sigma}_{ji}(t)]\\!]=\int\limits_{-\infty}^{t}dt_{1}[\\![G_{ij}(t,t_{1})K_{ij}(t_{1})]\\!],$
(VI.15)

where the classical average must include $K_{ij}(t_{1})$ since it is a
function of molecule operators, which also vary stochastically. For this
reason the expectation value in Eq. (VI.15) cannot be solved. Further, the
right-hand side involves $K_{ij}(t)$, which includes unknown molecule
operators. Therefore, to solve for the molecule operators where the right-hand
side involves “known” quantities we need to resort to perturbation theory. It
is at this point where each term will only depend on $\tilde{\omega}_{ij}(t)$
explicitly, allowing us to evaluate the expectation value.

Before resorting to perturbation theory we consider the exact formal
expressions for $\bar{\sigma}_{ge}(t)$ and $\bar{\sigma}_{ef}(t)$ that will be
used to calculate the absorption. Using the definition of $G_{ij}(t,t_{1})$,
exact expressions for $\bar{\sigma}_{ge}(t)$ and $\bar{\sigma}_{ef}(t)$ are

$\displaystyle\bar{\sigma}_{ge}(t)=\int\limits_{-\infty}^{t}dt_{1}[G_{eg}(t,t_{1})\bar{\sigma}_{gg}(t_{1})\hat{F}_{+}^{eg}(t_{1})$
(VI.16a) $\displaystyle\hskip
5.69054pt+G_{ge}^{*}(t,t_{1})\bar{\sigma}_{e^{\prime}e}(t_{1})\hat{F}_{+}^{e^{\prime}g}(t_{1})+G_{eg}(t,t_{1})\hat{F}_{-}^{ef}(t_{1})\bar{\sigma}_{gf}(t_{1})],$
$\displaystyle\bar{\sigma}_{ef}(t)=\int\limits_{-\infty}^{t}dt_{1}[G_{fe}(t,t_{1})\bar{\sigma}_{ee^{\prime}}(t_{1})\hat{F}_{+}^{fe^{\prime}}(t_{1})$
(VI.16b) $\displaystyle\hskip
5.69054pt+G^{*}_{ef}(t,t_{1})\hat{F}_{-}^{ge}(t_{1})\bar{\sigma}_{gf}(t_{1})+G^{*}_{ef}(t,t_{1})\bar{\sigma}_{f^{\prime}f}(t_{1})\hat{F}_{+}^{f^{\prime}e}(t_{1})],$

where we used Eq. (VI.14) and which both have three contributions on the
right-hand side, each corresponding to a distinct physical process that we
will discuss.

Figure 2: Exact diagrams for $\bar{\sigma}_{ge}(t)$ and
$\bar{\sigma}_{ef}(t)$. Time flows upwards such that $t\geq t_{1}$. Each
diagram starts with the exact form of an operator in a bubble at $t_{1}$ and
interacts with the field. The time evolution from $t_{1}\to t$ evolves
according to the Green function $G_{ij}(t,t_{1})$. Diagram I and IV contribute
to one- and two-photon absorption respectively. Diagrams II and V are
saturation terms of the one- and two-photon absorption respectively. The
physics of diagrams III and VI depend on the detuning, but far from resonance
the diagrams contribute to two-photon absorption.

It will be useful to introduce diagrammatic representations of the exact
expressions in Eqs. (VI.16a) and (VI.16b); see Fig. 2. These diagrams are
similar to the double-sided Feynman diagrams considered in many studies [2, 3,
49] and they follow similar rules. However, in our approach we work in the
Heisenberg picture as apposed to a density operator formalism, and include
decay processes from the beginning; this sets each of the molecule operators
in Eq. (VI.9) to zero until the interaction with the electromagnetic field is
introduced (see discussion under Eq. (VI.11).) When we begin perturbation
theory, it follows that the zeroth order solution in our perturbation
expansion is _not_ given by $\bar{\sigma}^{(0)}_{gg}\neq\ket{g}\bra{g}$, as we
discuss in detail below. Thus, the starting point of the perturbative
expansion of each diagram does not have a “$\ket{g}\bra{g}$” term. Our
diagrams involve the exact form of each coherence operator at time $t$ given
by $\bar{\sigma}_{ij}(t)$, and later when we begin perturbation theory we will
include a superscript to denote the order of perturbation theory, with
$\bar{\sigma}_{ij}^{(k)}(t)$, for example, denoting the $k$’th order in
perturbation theory.

To represent the fact that the coherence operators on the left- and right-hand
side of Eq. (VI.16a) and (VI.16b) are formal solutions we put them in a
“bubble”. When we begin to solve for them perturbatively we will provide a
perturbative expansion of each “bubble” diagram.

The rules of each diagram are as follows. Time increases upward. Inward
(outward) arrows on the right- and left-hand side represent absorption
(emission). Inward arrows on the bra (ket) side raise the bra (ket) side of
$\bar{\sigma}_{ij}(t_{1})$ and outward arrows on the bra (ket) side lower the
bra (ket) side of $\bar{\sigma}_{mn}(t_{1})$. Each inward/outward arrow is
connected to a vertical arrow and a Green function $G_{pq}(t,t_{1})$. If the
vertical arrow is on the right-hand (left-hand) side we include a factor of
$G_{qp}(t,t_{1})$ ($G^{*}_{pq}(t,t_{1})$) and the subscript is opposite (same)
of the $\bar{\sigma}_{qp}(t_{1})$ that it connects to. We note that because
the equations are normally ordered and the initial ket of the reservoir is the
vacuum state we always take
$\hat{F}^{ij}_{\pm}(t_{1})\to\bm{\mu}_{ij}\cdot\hat{\bm{E}}_{\bm{\pm}}(t_{1})$
in each diagram and read the dipole matrix elements on the bra (ket) side from
right to left (left to right).

With the aid of each diagram in Fig. 2 we can describe each contribution on
the right-hand side of $\bar{\sigma}_{ge}(t)$ and $\bar{\sigma}_{ef}(t)$ in
Eq. (VI.16); in later sections we will expand these simple descriptions with
more details.

First we consider the coherence operator $\bar{\sigma}_{ge}(t)$ (Eq. (VI.16a))
which has three contributions on the right-hand side. The first contribution
to the dynamics (I in Fig. 2) is given by the term
$\bar{\sigma}_{gg}(t)\hat{F}_{+}^{eg}(t)$; it depends on the population in the
ground state and an annihilation operator $\hat{F}_{+}^{eg}(t)$, which
facilitates transitions between $g$ and $e$. Taking the expectation value, we
understand the annihilation operator $\hat{F}_{+}^{eg}(t)$ as acting on the
right to annihilate a photon from the field, which raises the bra side of
$\bar{\sigma}_{gg}(t)\to\bar{\sigma}_{ge}(t)$. At lowest order in perturbation
theory this term leads to one-photon absorption and its higher order
contribution will be saturation.

The second contribution to the dynamics of $\bar{\sigma}_{ge}(t)$ (II in Fig.
2) is given by $\bar{\sigma}_{e^{\prime}e}(t)\hat{F}_{+}^{e^{\prime}g}(t)$; it
depends on coherences between intermediate states and the annihilation
operator $\hat{F}_{+}^{e^{\prime}g}(t)$, which facilitates transitions between
$g$ and $e^{\prime}$. The annihilation operator $\hat{F}_{+}^{e^{\prime}g}(t)$
acts on the left creating an excitation in the field-reservoir system, which
lowers the ket side of $\bar{\sigma}_{e^{\prime}e}(t)\to\bar{\sigma}_{ge}(t)$.
At lowest order this term describes a saturation of the one-photon absorption.

The third and final contribution to the dynamics of $\bar{\sigma}_{ge}(t)$
(III in Fig. 2) is given by $\hat{F}_{-}^{ef}(t)\bar{\sigma}_{gf}(t)$; it
depends on the field-reservoir creation operator $\hat{F}_{-}^{ef}(t)$, which
facilitates transitions between $e$ and $f$, and on the coherence between $g$
and $f$. Coherence between $g$ and $e$ is generated from coherence between $g$
and $f$ and the creation operator $\hat{F}_{-}^{ef}(t)$ acting on the right,
which creates an excitation in the field-reservoir system and lowers the bra
side of $\bar{\sigma}_{gf}(t)\to\bar{\sigma}_{ge}(t)$. In later sections we
will find that this term’s contribution will depend largely on the detunings
from resonances, and a simple characterization of “one-photon absorption” or
“two-photon absorption” seems not possible. However, when the incident field
is far detuned from resonance with an intermediate state, this term can be
identified as a contribution to two-photon absorption.

Next we consider the coherence operator $\bar{\sigma}_{ef}(t)$ (Eq. (VI.16b))
which has three contributions on the right-hand side. The first contribution
to the dynamics (IV in Fig. 2) is given by the term
$\bar{\sigma}_{ee^{\prime}}(t)\hat{F}_{+}^{fe^{\prime}}(t)$; it depends on
coherence between excited states and the annihilation operator
$\hat{F}_{+}^{fe^{\prime}}(t)$, which facilitates transitions between
$e^{\prime}$ and $f$. Taking the expectation value, we understand the
annihilation operator $\hat{F}_{+}^{fe^{\prime}}(t)$ as acting on the right to
annihilate a photon from the field, which raises the bra side of
$\bar{\sigma}_{ee^{\prime}}(t)\to\bar{\sigma}_{ef}(t)$. At lowest order in
perturbation theory this term contributes to two-photon absorption.

The second contribution to the dynamics of $\bar{\sigma}_{ef}(t)$ (V in Fig.
2) is given by $\hat{F}_{-}^{ge}(t)\bar{\sigma}_{gf}(t)$; it depends on the
field-reservoir creation operator $\hat{F}_{-}^{ge}(t)$, which facilitates
transitions between $g$ and $e$, and on the coherence between $g$ and $f$.
Coherence is generated between $e$ and $f$ from coherence between $g$ and $f$
and the creation operator $\hat{F}_{-}^{ge}(t)$ acting on the left, which
annihilates a photon from the field and raises the ket side of
$\bar{\sigma}_{gf}(t)\to\bar{\sigma}_{ef}(t)$. Similar to the term shown in
III in Fig. 2, its contribution depends strongly on the detunings from
resonances, and does not allow a simple characterization. But like III, when
the incident field is far detuned from a resonance with an intermediate state
this term can be identified as a contribution two-photon absorption.

The third and final contribution to the dynamics of $\bar{\sigma}_{ef}(t)$ (VI
in Fig. 2) is given by
$\bar{\sigma}_{f^{\prime}f}(t)\hat{F}_{+}^{f^{\prime}e}(t)$; it depends on
coherences between final states and the annihilation operator
$\hat{F}_{+}^{f^{\prime}e}(t)$, which facilitates transitions between $e$ and
$f^{\prime}$. We understand this term as modifying the coherence
$\bar{\sigma}_{ef}(t)$ by the annihilation operator
$\hat{F}_{+}^{f^{\prime}e}(t)$ acting on the left to create an excitation in
the field-reservoir system, which lowers the ket side of
$\bar{\sigma}_{f^{\prime}f}(t)\to\bar{\sigma}_{ef}(t)$. At lowest order this
term describes a saturation of the two-photon absorption; however, at the
order of perturbation we will be considering, it will not contribute.

##  VII perturbative solution

To solve the system of equations (VI.9) we expand each $\bar{\sigma}_{ij}(t)$
in a perturbative expansion given by

$\bar{\sigma}_{ij}(t)=\bar{\sigma}_{ij}^{(0)}(t)+\lambda\bar{\sigma}_{ij}^{(1)}(t)+\lambda^{2}\bar{\sigma}_{ij}^{(2)}(t)+...$
(VII.1)

and take $\hat{\bm{E}}_{\bm{\pm}}(t)\to\lambda\hat{\bm{E}}_{\bm{\pm}}(t)$
where $\lambda$ characterizes the order of the perturbation. We use the
perturbative expansion (VII.1) with the exact formal solution (VI.13),
identify the zeroth order solution, and iterate order-by-order.

For the zeroth order solution, since the initial ket of the reservoir is the
vacuum state and the equations are normal ordered, no coherence is generated
and no population is moved out of the ground state. Then using the identity
relation (II.2), which is also satisfied in the interaction picture we work
in, we solve for the last unknown molecule operator, and for finite times we
have

$\bar{\sigma}_{gg}^{(0)}(t)=\hat{1}.$ (VII.2)

This occurs because the interaction between the molecule and reservoir is
still operative, and since we take the initial ket of the reservoir to be
vacuum, the molecule can only lose energy to the reservoir. So for
$\ket{\Psi_{\text{in}}}$ describing the molecule in any state, for finite
times $t>-\infty$ the molecule will move to the ground state, and Eq. (VII.2)
necessarily results.

To generate each higher order term we use the zeroth order term in the right-
hand side of each exact result and then continue order by order. In Appendix B
we begin with the exact solutions for each molecule operator and apply the
steps outlined above to third order in $\bar{\sigma}_{ge}(t)$ and
$\bar{\sigma}_{ef}(t)$. The nonzero terms that contribute to the absorption
are $\bar{\sigma}_{ge}^{(1)}(t),\bar{\sigma}_{ge}^{(3)}(t)$ and
$\bar{\sigma}_{ef}^{(3)}(t)$ and are given in Eq. (B.2), (B.4) and (B.5)
respectively.

Each coherence operator
$\bar{\sigma}_{ge}^{(1)}(t),\bar{\sigma}_{ge}^{(3)}(t)$ and
$\bar{\sigma}_{ef}^{(3)}(t)$ still involves the stochastic functions
$\tilde{\omega}_{ij}(t)$, which are included in the definition of
$G_{ij}(t,t_{1})$. To push forward we need to take the average over the
classical distributions. In Appendix C we begin with the simplest case and
work out the expectation value of $G_{ij}(t,t_{1})$ using the assumptions
after Eq. (II.13). Taking the “impact limit” in which correlations decay on
fast time scales compared to the coupling strength [35, 46],
$\tau_{c}^{2}c_{ij}\ll 1$, we find $[\\![G_{ij}(t,t_{1})]\\!]=G_{ij}(t-t_{1})$
where

$G_{ij}(t-t_{1})\equiv ie^{-(\gamma_{ij}+i\omega_{ij})(t-t_{1})},$ (VII.3)

where we defined the total decay width

$\gamma_{ij}\equiv\bar{\Gamma}_{ij}+\Lambda_{ij},$ (VII.4)

where $\Lambda_{ij}$ is the dephasing decay rate between state $i$ and $j$
defined in Appendix C and satisfies $\Lambda_{ij}=\Lambda_{ji}$ and
$\Lambda_{ii}=0$. Since each contribution to $\gamma_{ij}$ is symmetric so is
$\gamma_{ij}$. Then the new Green function defined in Eq. (VII.3) also
satisfies

$G^{*}_{ij}(t-t_{1})=-G_{ji}(t-t_{1}).$ (VII.5)

The result for $[\\![G_{ij}(t,t_{1})]\\!]$ can be directly applied to the
stochastic average for $\bar{\sigma}_{ge}^{(1)}(t)$ because there is only one
$G_{eg}(t,t_{1})$ present. However, for the higher order molecule operators
there are three Green functions present and we must take the expectation value
over the product of them which can become quite complicated. If we assume that
the correlation between random variables which decays on time scales given by
$\tau_{c}$ is much shorter than the dephasing rates $\Lambda_{ij}$, for all
$i$ and $j$, i.e. $\Lambda_{ij}\tau_{c}\ll 1$, then we can apply the
“factorization approximation” [46, 35] where the expectation value of a
product of green functions is equal to the product of expectation values, that
is

$[\\![G_{ij}(t,t_{n})\dots G_{pq}(t_{2},t_{1})]\\!]=G_{ij}(t-t_{n})\dots
G_{pq}(t_{2}-t_{1}).$ (VII.6)

Thus to take the average value of any coherence operator we replace each Green
function with Eq. (VII.3), which includes the total rate $\gamma_{ij}$.

After applying the above prescription and taking the stochastic average we
have solved for the perturbative expansion of $\bar{\sigma}_{ge}(t)$ and
$\bar{\sigma}_{ef}(t)$ up to third order and in doing so we also solved for
$\bar{\sigma}_{gg}^{(2)}(t)$, $\bar{\sigma}_{e^{\prime}e}^{(2)}(t)$ and
$\bar{\sigma}_{gf}^{(2)}(t)$.

In Fig. 2 we included bubbles to denote the formal solution of each coherence
operator. Using the perturbative solution of $\bar{\sigma}_{gg}(t)$,
$\bar{\sigma}_{e^{\prime}e}(t)$ and $\bar{\sigma}_{gf}(t)$ we can fill in each
bubble with a perturbative diagram shown in Fig. 3, which we now describe.

The first diagram is shown in Fig. 3a and is the perturbative expansion of the
operator $\bar{\sigma}_{gg}(t)$. The zeroth order term is Eq. (VII.2), while
the second order solution involves the sum of two contributions, absorption
and emission on the bra and ket side. Since the result of the second order
term involves an emission of a photon, the second order term
$\bar{\sigma}_{gg}^{(2)}(t)$ is a correction to the zeroth order term,
lowering the amplitude that the molecule is in the ground state. Thus when
considering their contribution to the dynamics of $\bar{\sigma}_{ge}(t)$, the
zeroth order term leads to one-photon absorption, and the second order term
contributes to its saturation.

The second diagram shown in in Fig. 3b is the perturbative expansion for
$\bar{\sigma}_{e^{\prime}e}(t)$, which is only nonzero at second order and has
two contributions. In the terminology of nonlinear spectroscopy the first term
leads to the “rephasing” contribution and the second to the “nonrephasing”
contribution [50, 51]. Both of these terms combine to give coherence between
intermediate states and when $e=e^{\prime}$ amplitude in the intermediate
state. Thus when considering the contribution to the dynamics of
$\bar{\sigma}_{ge}(t)$ we expect it to lead to saturation, and to that of
$\bar{\sigma}_{ef}(t)$ we expect it to lead to two-photon absorption.

The third and final diagram is for the perturbative expansion of
$\bar{\sigma}_{gf}(t)$ given by Fig. 3.c) and has one contribution which is
given by the absorption of two photons on the bra side. This terms leads to a
virtual two-photon transition to the excited state $f$, which has been
referred to as the “double quantum coherence” [32], and the effect that this
term has on $\bar{\sigma}_{ge}(t)$ and $\bar{\sigma}_{ef}(t)$ is complicated
in general, but far from resonance will lead to two-photon absorption.

##  VIII absorption expressions

We are now at the point where we are able to calculate the absorption. The
general expressions that result are given in Appendix D. Here we define the
terms involved and the steps necessary; in the next section we deal with
simplified limits of the general expressions in an example calculation.

We begin by taking the classical expectation value which is done by replacing
$G_{ij}(t,t_{1})\to G_{ij}(t-t_{1})$ and drop the average value symbol
$[\\![\cdot]\\!]$ to avoid extra notation. Then we input the coherence
operators $\bar{\sigma}_{ge}^{(1)}(t)$,$\bar{\sigma}_{ge}^{(3)}(t)$ and
$\bar{\sigma}_{ef}^{(3)}(t)$ given in Eq. (B.2), (B.4), (B.5) respectively
into Eq. (III.2) for the absorption. The result involves iterated time
integrals (Eq. (D.1)).

Figure 3: Perturbative expansion of each bubble diagram on the left.

To consider different incident states of light it is useful to move to
frequency space to evaluate this expression. For any Hermitian operator
$O(t)$, such as the electric field operator at the site of the molecule, we
write

$O(t)=\int\limits_{-\infty}^{\infty}d\bar{}\hskip 1.00006pt\omega
O(\omega)e^{-i\omega t},$ (VIII.1)

where we define $d\bar{}\hskip 1.00006pt\omega=d\omega/\sqrt{2\pi}$.
Introducing positive and negative frequency components in the usual way we
then write all frequency integrals to range from $0\to\infty$, and leave these
limits implicit in the expressions below.

Evaluating each time integral introduces resonant frequency denominator
contributions

$Q_{ij}(\omega)\equiv\omega_{ij}-\omega-i\gamma_{ij},$ (VIII.2)

which follow from the exponent of each Green function $G_{ij}(t-t_{1})$. And
the “nonrephasing” and “rephasing” contributions to
$\bar{\sigma}_{e^{\prime}e}(t)$ mentioned above can be combined into a single
term, involving a new function

$R_{ee^{\prime}}(\omega)\equiv\frac{\omega_{ee^{\prime}}-\omega-i\gamma_{eg}-i\gamma_{e^{\prime}g}}{\omega_{ee^{\prime}}-\omega-i\gamma_{ee^{\prime}}}.$
(VIII.3)

The term $R_{ee^{\prime}}(\omega)$ deviates from unity only because we
consider two broadening mechanisms, population decay and dephasing. In the
limit when dephasing is negligible ($\gamma_{eg}\to\bar{\Gamma}_{eg}$),
$R_{ee^{\prime}}(\omega)\to 1$. Physically, $R_{ee^{\prime}}(\omega)$
corresponds to contributions to the absorption that involve either a
population term or a coherence between excited states.

With the use of Eqs (VIII.2) and (VIII.3), our expression involving iterated
time integrals (Eq. (D.1)) reduces to one involving frequency integrals (Eq.
(D.2)), which we write in terms of $\hat{F}_{\pm}^{ij}(\omega)$ to keep the
notation simpler. But since the initial ket of the reservoir is the vacuum
state,

$\hat{F}_{+}^{ij}(\omega)\ket{\Psi_{\text{in}}}=\frac{\bm{\mu}_{ij}}{\hbar}\cdot\hat{\bm{E}}_{+}(\omega)\ket{\Psi_{\text{in}}},$
(VIII.4)

since each absorption term is normally ordered. Transitions are thus only made
by the external field, as expected.

This expression (Eq. (D.2)) for the absorption of a molecule with many
intermediate and final states; with possible resonant excitation to
intermediate state and with two sources of broadening, is a main result of
this paper. It is valid near resonance and for any incident state of light, in
general propagating in three dimensions with any polarization, and either
pulsed or CW. However, for our sample calculation in the following section we
adopt a usual approximation and consider a field essentially uniform over a
cross-section area $A$ in the $xy$ plane, and propagating in vacuum in the $z$
direction. The positive frequency components of the electric field operator
within the cross-section area $A$, which in our case is sometimes called the
“entanglement area” [52, 24, 29, 31, 1], are given by

$\hat{\bm{E}}_{\bm{+}}(z,\omega)=\bm{e}\sqrt{\frac{\hbar\omega}{2\epsilon_{0}cA}}a(\omega)e^{ik(\omega)z},$
(VIII.5)

where we assumed the incident field has a single polarization given by
$\bm{e}$, the wave-vector is given by $k(\omega)=\omega/c$, and the creation
and annihilation operators satisfy the commutation relations

$\left[a(\omega),a^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime}),$
(VIII.6)

see, e.g., Blow et al. [53]. In our sample calculations below, we consider a
single molecule at the origin, and assume that the incident fields we consider
have a centre frequency $\bar{\omega}$ and a bandwidth much smaller than their
centre frequency. Then the positive frequency component of the electric field
operator that appears in the dipole interaction can be taken as

$\bm{\hat{E}_{+}}(\omega)=\bm{e}Ea(\omega),$ (VIII.7)

where

$E\equiv\sqrt{\frac{\hbar\bar{\omega}}{2\epsilon_{0}cA}}.$ (VIII.8)

##  IX Three level molecule

As a sample calculation where the results can be given in closed form, we
consider a three level molecule with a ground state $g$, single intermediate
state $e$ and single final state $f$. As in Section VIII we consider the
incident field is given by a single polarization $\bm{e}$, and with that
assumption we define the two constants

$\displaystyle d_{eg}\equiv\frac{\bm{\mu}_{eg}\cdot\bm{e}}{\hbar},$ (IX.1a)
$\displaystyle d_{fe}\equiv\frac{\bm{\mu}_{fe}\cdot\bm{e}}{\hbar},$ (IX.1b)

which are the dipole transition elements scaled by $\hbar$.

Using these constants and our expression (VIII.7) for the field operator, the
absorption (Eq. (D.2)) can be written as the sum of five contributions,
$\mathcal{A}=\mathcal{A}^{1}+\mathcal{A}^{2}+\mathcal{A}^{3}+\mathcal{A}^{4}+\mathcal{A}^{5}$,
where

$\displaystyle\mathcal{A}^{1}=\int
d\omega\text{Im}\left[\mathcal{R}^{1}(\omega)\right]\frac{\hbar\bar{\omega}\langle
a^{\dagger}(\omega)a(\omega)\rangle}{A},$ (IX.2a)
$\displaystyle\mathcal{A}^{i}=2\pi\int d\bar{}\hskip
1.00006pt\omega_{1}d\bar{}\hskip 1.00006pt\omega_{2}d\bar{}\hskip
1.00006pt\omega_{3}d\bar{}\hskip
1.00006pt\omega_{4}\text{Im}\left[\mathcal{R}^{i}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})\frac{(\hbar\bar{\omega})^{2}\langle
a^{\dagger}(\omega_{4})a^{\dagger}(\omega_{3})a(\omega_{2})a(\omega_{1})\rangle}{A^{2}}\right]\delta(\omega_{1}+\omega_{2}-\omega_{3}-\omega_{4}),$
(IX.2b)

for $i=2,3,4,5$, with the “broadening functions” given by

$\displaystyle\mathcal{R}^{1}(\omega_{1})=\frac{\hbar\bar{\omega}}{\epsilon_{0}c}\frac{|d_{eg}|^{2}}{Q_{eg}(\omega_{1})},$
(IX.3a)
$\displaystyle\mathcal{R}^{2}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})=\frac{-\hbar\bar{\omega}|d_{eg}|^{4}R_{ee}(\omega_{2}-\omega_{3})}{(\epsilon_{0}c)^{2}Q_{eg}(\omega_{4})Q_{eg}^{*}(\omega_{3})Q_{eg}(\omega_{2})},$
(IX.3b)
$\displaystyle\mathcal{R}^{3}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})=\frac{\hbar\bar{\omega}|d_{fe}d_{eg}|^{2}}{2(\epsilon_{0}c)^{2}Q_{eg}(\omega_{4})Q_{fg}(\omega_{1}+\omega_{2})Q_{eg}(\omega_{2})},$
(IX.3c)
$\displaystyle\mathcal{R}^{4}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})=\frac{\hbar\bar{\omega}|d_{fe}d_{eg}|^{2}R_{ee}(\omega_{2}-\omega_{3})}{2(\epsilon_{0}c)^{2}Q_{fe}(\omega_{4})Q_{eg}^{*}(\omega_{3})Q_{eg}(\omega_{2})},$
(IX.3d)
$\displaystyle\mathcal{R}^{5}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})=\frac{-\hbar\bar{\omega}|d_{fe}d_{eg}|^{2}}{2(\epsilon_{0}c)^{2}Q_{fe}(\omega_{4})Q_{fg}(\omega_{1}+\omega_{2})Q_{eg}(\omega_{2})}.$
(IX.3e)

Each absorption term is then given as a product of two functions. The lowest
order term $\mathcal{A}^{1}$ consists of the imaginary part of the broadening
function $\mathcal{R}^{1}(\omega)$ multiplied by the energy density/area of
the incident field given by the one-photon correlation function $\langle
a^{\dagger}(\omega)a(\omega)\rangle$, containing all the one-photon
statistics. Each higher order term $\mathcal{A}^{i}$ ($i=2,3,4,5$) is given by
the imaginary part of the higher order broadening function
$\mathcal{R}^{i}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ multiplied by
the $\text{(energy density/area)}^{2}$ given by the four point correlation
function $\langle
a^{\dagger}(\omega_{4})a^{\dagger}(\omega_{3})a(\omega_{2})a(\omega_{1})\rangle$,
containing all the two-photon statistics of the incident light.

The absorption terms $\mathcal{A}^{1}$ and $\mathcal{A}^{2}$ are generated
from the perturbative expansion of the diagram in Fig. 2.I and 2.II and
contribute to the one-photon absorption, the second term describing the onset
of saturation. The fourth absorption term $\mathcal{A}^{4}$ is generated from
the perturbative expansion of the diagrams in Fig. 2.IV and contributes to
two-photon absorption. The third and fifth absorption terms, $\mathcal{A}^{3}$
and $\mathcal{A}^{5}$, are generated from the perturbative expansion of the
diagrams in Fig. 2.III and in Fig. 2.V, respectively, and are more
complicated. But in the far detuned limit they will also contribute to two-
photon absorption.

We now consider different states of light and calculate each term. Our
interest in this paper is in the CW limit.

##  X CW coherent light

We begin with coherent light, described by a state

$\ket{\alpha}=D(\alpha)\ket{\text{vac}},$ (X.1)

where $\alpha$ is a complex coefficient and the displacement operator

$D(\alpha)\equiv e^{\alpha\int
d\omega\varphi(\omega)a^{\dagger}(\omega)-\text{H.c.}}.$ (X.2)

The spectral profile of the state is given by $\varphi(\omega)$ and is
normalized to satisfy

$\int dt|\varphi(t)|^{2}=\int d\omega|\varphi(\omega)|^{2}=1.$ (X.3)

To model coherent light in the CW limit, we consider a field oscillating at
the frequency $\bar{\omega}$ for a length of time $T_{p}$. The normalized
distribution in time and frequency is then

$\displaystyle\varphi(t)=\frac{e^{-i\bar{\omega}t}}{\sqrt{T_{p}}},\hskip
14.22636pt-\frac{T_{p}}{2}\leq t\leq\frac{T_{p}}{2},$ (X.4a)
$\displaystyle\varphi(\omega)=\frac{1}{\sqrt{\Omega_{p}}}\operatorname{sinc}\left(\frac{(\omega-\bar{\omega})\pi}{\Omega_{p}}\right),$
(X.4b)

where we take $\Omega_{p}=2\pi/T_{p}$ to identify an effective bandwidth. The
CW limit is then $T_{p}\to\infty$, $\Omega_{p}\to 0$, yielding a spectral
profile $\varphi(\omega)$ strongly peaked at $\omega=\bar{\omega}$ with the
squared spectral profile given by [54]

$\lim_{\Omega_{p}\to 0}|\varphi(\omega)|^{2}=\delta(\omega-\bar{\omega}).$
(X.5)

Taking the CW limit will lead to infinite absorption, since the light is
always ‘on,’ and in the calculation we divide each absorption term by $T_{p}$
to get a finite result as $T_{p}\to\infty$ that can be identified as an
absorption rate.

Using the transformation property

$D^{\dagger}(\alpha)a(\omega)D(\alpha)=a(\omega)+\alpha\varphi(\omega),$ (X.6)

and the unitarity of the displacement operator, we can calculate the two
correlation functions we need. The one-photon correlation function is given by

$\bra{\alpha}a^{\dagger}(\omega)a(\omega)\ket{\alpha}=|\alpha|^{2}\varphi^{*}(\omega)\varphi(\omega),$
(X.7)

and the two-photon correlation function by

$\begin{split}&\bra{\alpha}a^{\dagger}(\omega_{4})a^{\dagger}(\omega_{3})a(\omega_{2})a(\omega_{1})\ket{\alpha}\\\
&\hskip
99.58464pt=|\alpha|^{4}\varphi^{*}(\omega_{4})\varphi^{*}(\omega_{3})\varphi(\omega_{2})\varphi(\omega_{1}).\end{split}$
(X.8)

In the CW limit both these quantities will peak and be significantly different
from zero only when their argument(s) are close to $\bar{\omega}$.

From the one-photon correlation function we can immediately identify the
energy of the light, which is given by

$E_{\text{coh}}=\hbar\bar{\omega}|\alpha|^{2},$ (X.9)

where as usual $|\alpha|^{2}$ is the average number of photons in the coherent
state. We will let this diverge as $T_{p}\to\infty$ so the intensity in the
pulse does not change in that limit, introducing a fixed parameter
$\alpha_{0}$ by setting $\alpha=\alpha_{0}\sqrt{T_{p}}$. Then
$F_{\text{coh}}\equiv|\alpha_{0}|^{2}/A$ is identified as the photon flux; the
intensity of the beam is then

$I_{\text{coh}}=\hbar\bar{\omega}F_{\text{coh}}.$ (X.10)

Using the result of the correlation functions in Eq. (X.7) and (X.8), and that
in the CW limit $\varphi(\omega)$ is strongly peaked at $\omega=\bar{\omega}$,
the rate of each absorption term for coherent light is given by

$\displaystyle\frac{\mathcal{A}^{1}_{\text{coh}}}{T_{p}}=I_{\text{coh}}\sigma^{1}_{\text{coh}}$
(X.11a)
$\displaystyle\frac{\mathcal{A}^{i}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\sigma^{i}_{\text{coh}},$
(X.11b)

for $i=2,3,4,5$, where each cross-section is

$\displaystyle\sigma^{1}_{\text{coh}}=\text{Im}\left[\mathcal{R}^{1}(\bar{\omega})\right]$
(X.12a)
$\displaystyle\sigma^{i}_{\text{coh}}=\text{Im}\left[\mathcal{R}^{i}(\bar{\omega},\bar{\omega},\bar{\omega},\bar{\omega})\right].$
(X.12b)

The one-photon term scales with intensity and each higher order two-photon
term scales with the intensity squared, as expected. We note that each
$\sigma^{i}_{\text{coh}}$ for $i=1,2,3,4,5$ is independent of the beam
parameters aside from frequency. The quantity $\sigma^{1}_{\text{coh}}$ is the
lowest order cross-section due to a one-photon transition, and has units of
area. Each higher order cross-section $\sigma^{i}_{\text{coh}}$ for
$i=2,3,4,5$ has units of $\text{area}/\text{intensity}$.

Finally, we calculate the rate of absorption in full by evaluating and
expanding each cross-section term in Eq. (X.11) with the appropriate
broadening functions in Eq. (IX.3), then each absorption term is given by

$\displaystyle\frac{\mathcal{A}^{1}_{\text{coh}}}{T_{p}}=I_{\text{coh}}\frac{|d_{eg}|^{2}\hbar\bar{\omega}}{\epsilon_{0}c}\frac{\gamma_{eg}}{(\omega_{eg}-\bar{\omega})^{2}+\gamma_{eg}^{2}},$
(X.13a)
$\displaystyle\frac{\mathcal{A}^{2}_{\text{coh}}}{T_{p}}=-I_{\text{coh}}^{2}\frac{|d_{eg}|^{4}}{(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}^{2}}{[(\omega_{eg}-\bar{\omega})^{2}+\gamma_{eg}^{2}]^{2}},$
(X.13b)
$\displaystyle\frac{\mathcal{A}^{3}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}\hbar\bar{\omega}}{2(\epsilon_{0}c)^{2}}\frac{(\omega_{eg}-\bar{\omega})^{2}\gamma_{fg}+2(\omega_{eg}-\bar{\omega})(\omega_{fg}-2\bar{\omega})\gamma_{eg}-\gamma_{fg}\gamma_{eg}^{2}}{[(\omega_{eg}-\bar{\omega})^{2}+\gamma_{eg}^{2}]^{2}[(\omega_{fg}-2\bar{\omega})^{2}+\gamma_{fg}^{2}]},$
(X.13c)
$\displaystyle\frac{\mathcal{A}^{4}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}}{(\omega_{eg}-\bar{\omega})^{2}+\gamma_{eg}^{2}}\frac{\gamma_{fe}}{(\omega_{fe}-\bar{\omega})^{2}+\gamma_{fe}^{2}},$
(X.13d)
$\displaystyle\frac{\mathcal{A}^{5}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}\hbar\bar{\omega}}{2(\epsilon_{0}c)^{2}}\frac{\gamma_{eg}\gamma_{fe}\gamma_{fg}-(\omega_{eg}-\bar{\omega})(\omega_{fe}-\bar{\omega})\gamma_{fg}-(\omega_{fg}-2\bar{\omega})[(\omega_{eg}-\bar{\omega})\gamma_{fe}+(\omega_{fe}-\bar{\omega})\gamma_{eg}]}{[(\omega_{eg}-\bar{\omega})^{2}+\gamma_{eg}^{2}][(\omega_{fe}-\bar{\omega})^{2}+\gamma_{fe}^{2}][(\omega_{fg}-2\bar{\omega})^{2}+\gamma_{fg}^{2}]}.$
(X.13e)

From equation (X.13) it is clear that the contribution of the terms
$\mathcal{A}^{3}_{\text{coh}}$ and $\mathcal{A}^{5}_{\text{coh}}$ depends
strongly on the relative magnitudes of the detunings from resonance and decay
parameters. In certain regimes $\mathcal{A}^{3}_{\text{coh}}$ and
$\mathcal{A}^{5}_{\text{coh}}$ can be either positive or negative, and thus
can contribute to either absorption or saturation.

###  X.1 Far detuned from intermediate state resonances

To get a better understanding of $\mathcal{A}^{3}_{\text{coh}}$ and
$\mathcal{A}^{5}_{\text{coh}}$, and to connect with the literature, we
consider the limit when the field is far detuned from any resonances involving
the intermediate level. That is, we assume
$|\omega_{eg}-\bar{\omega}|\gg\gamma_{eg}$ and
$|\omega_{fe}-\bar{\omega}|\gg\gamma_{fe}$; here one-photon absorption will be
small, and two-photon absorption will dominate for large intensities. Defining
$\delta_{eg}\equiv\omega_{eg}-\bar{\omega}$ and
$\delta_{fe}\equiv\omega_{fe}-\bar{\omega}$ to be detunings from the
resonances (see Fig. 4.a for molecule diagram), in the limit when
$|\delta_{eg}|\gg\gamma_{eg}$ and $|\delta_{fe}|\gg\gamma_{fe}$, to good
approximation the absorption terms are

$\displaystyle\frac{\mathcal{A}^{1}_{\text{coh}}}{T_{p}}=I_{\text{coh}}\frac{|d_{eg}|^{2}\hbar\bar{\omega}}{\epsilon_{0}c}\frac{\gamma_{eg}}{\delta_{eg}^{2}},$
(X.14a)
$\displaystyle\frac{\mathcal{A}^{2}_{\text{coh}}}{T_{p}}=-I_{\text{coh}}^{2}\frac{|d_{eg}|^{4}}{(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}^{2}}{\delta_{eg}^{4}},$
(X.14b)
$\displaystyle\frac{\mathcal{A}^{3}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\delta_{eg}^{2}}\frac{\gamma_{fg}}{(\omega_{fg}-2\bar{\omega})^{2}+\gamma_{fg}^{2}},$
(X.14c)
$\displaystyle\frac{\mathcal{A}^{4}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}}{\delta_{eg}^{2}}\frac{\gamma_{fe}}{\delta_{fe}^{2}},$
(X.14d)
$\displaystyle\frac{\mathcal{A}^{5}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\delta_{eg}(-\delta_{fe})}\frac{\gamma_{fg}}{(\omega_{fg}-2\bar{\omega})^{2}+\gamma_{fg}^{2}}.$
(X.14e)

Here the first two terms describe one-photon absorption and its saturation,
and the three remaining terms describe two-photon absorption. To maximize the
two-photon absorption we take the incident field to be resonant with the two-
photon transition, $2\bar{\omega}=\omega_{fg}$. In this limit the two
detunings $\delta_{eg}$ and $\delta_{fe}$ simplify to a single detuning
parameter $\delta$ which we define by $\delta\equiv\delta_{eg}$ and
$\delta_{fe}=-\delta$ (see Fig. 4.b for special case molecule diagram). Then
each absorption term is further simplified to

$\displaystyle\frac{\mathcal{A}^{1}_{\text{coh}}}{T_{p}}=I_{\text{coh}}\frac{|d_{eg}|^{2}\hbar\bar{\omega}}{\epsilon_{0}c}\frac{\gamma_{eg}}{\delta^{2}},$
(X.15a)
$\displaystyle\frac{\mathcal{A}^{2}_{\text{coh}}}{T_{p}}=-I_{\text{coh}}^{2}\frac{|d_{eg}|^{4}}{(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}^{2}}{\delta^{4}},$
(X.15b)
$\displaystyle\frac{\mathcal{A}^{3}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\delta^{2}}\frac{1}{\gamma_{fg}},$
(X.15c)
$\displaystyle\frac{\mathcal{A}^{4}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\bar{\Gamma}_{eg}}\frac{\gamma_{eg}\gamma_{fe}}{\delta^{4}},$
(X.15d)
$\displaystyle\frac{\mathcal{A}^{5}_{\text{coh}}}{T_{p}}=I_{\text{coh}}^{2}\frac{|d_{eg}d_{fe}|^{2}}{2(\epsilon_{0}c)^{2}}\frac{\hbar\bar{\omega}}{\delta^{2}}\frac{1}{\gamma_{fg}}.$
(X.15e)

Each absorption term now has a simple scaling with the detuning $\delta$.
Higher order population terms $\mathcal{A}^{2}_{\text{coh}}$ and
$\mathcal{A}^{4}_{\text{coh}}$ scale as $\delta^{-4}$ and the virtual two-
photon transitions terms $\mathcal{A}^{3}_{\text{coh}}$ and
$\mathcal{A}^{5}_{\text{coh}}$ scale as $\delta^{-2}$.

Figure 4: Three level molecule diagram for a) the general case, and b) when
resonant with the two-photon transition.

In this limit we identify the rate of one-photon absorption

$\frac{\mathcal{O}_{\text{coh}}}{T_{p}}=\frac{\mathcal{A}^{1}_{\text{coh}}+\mathcal{A}^{2}_{\text{coh}}}{T_{p}},$
(X.16)

where $\mathcal{A}^{1}_{\text{coh}}$ is the lowest order term and
$\mathcal{A}^{2}_{\text{coh}}$ is the higher order correction due to
saturation, and the rate of two-photon absorption by

$\frac{\mathcal{T}_{\text{coh}}}{T_{p}}=\frac{\mathcal{A}^{3}_{\text{coh}}+\mathcal{A}^{4}_{\text{coh}}+\mathcal{A}^{5}_{\text{coh}}}{T_{p}}.$
(X.17)

The second contribution to the two-photon absorption is
$\mathcal{A}^{4}_{\text{coh}}$ which is generated from population in the
intermediate state. The first and third contribution is
$\mathcal{A}^{3}_{\text{coh}}+\mathcal{A}^{5}_{\text{coh}}$, which in this
limit are identical, and are generated by coherence between $g$ and $f$
through a virtual two-photon transition.

For our perturbative expansion to be reasonable we must have
$|\mathcal{A}^{1}_{\text{coh}}|\gg|\mathcal{A}^{2}_{\text{coh}}|$, so that the
first saturation correction of the one-photon absorption is small and higher
order corrections can be justifiably neglected. So we will restrict the
possible detunings and intensities we consider so that the magnitude of
$\mathcal{A}^{2}_{\text{coh}}$ satisfies

$|\mathcal{A}^{2}_{\text{coh}}|\leq 0.1|\mathcal{A}^{1}_{\text{coh}}|,$ (X.18)

If this restriction holds, then the rate of one-photon absorption is to good
approximation given by $\mathcal{O}_{\text{coh}}/T_{p}\simeq
I_{\text{coh}}\sigma^{1}_{\text{coh}}$, where

$\sigma^{1}_{\text{coh}}=\hbar\bar{\omega}\frac{|d_{eg}|^{2}}{\epsilon_{0}c}\frac{\gamma_{eg}}{\delta^{2}},$
(X.19)

identifying the one-photon absorption cross-section
$\sigma^{\text{1PA}}_{\text{coh}}\equiv\sigma^{1}_{\text{coh}}$, in agreement
with past results [49, 32].

Moving to the two-photon absorption in the far detuned limit where
$|\delta|\gg\gamma_{fe},\gamma_{ge}$, we find
$\mathcal{A}^{3}_{\text{coh}}+\mathcal{A}^{5}_{\text{coh}}\gg\mathcal{A}^{4}_{\text{coh}}$
(assuming $\gamma_{fe}\simeq\gamma_{fg}$). So to good approximation
$\mathcal{T}_{\text{coh}}/T_{p}\simeq
I_{\text{coh}}^{2}\sigma_{\text{coh}}^{\text{2PA}}$, where
$\sigma_{\text{coh}}^{\text{2PA}}$ is the two-photon absorption cross-section
given by

$\sigma_{\text{coh}}^{\text{2PA}}\simeq\sigma_{\text{coh}}^{3}+\sigma_{\text{coh}}^{5},$
(X.20)

which has units of area/intensity. Equivalently, we can define each cross-
section by

$\acute{\sigma}_{\text{coh}}^{i}=\frac{\hbar\bar{\omega}}{2}\sigma_{\text{coh}}^{i},$
(X.21)

and the two-photon absorption cross-section by
$\acute{\sigma}_{\text{coh}}^{\text{2PA}}\simeq\acute{\sigma}_{\text{coh}}^{3}+\acute{\sigma}_{\text{coh}}^{5}$,
which has units of $\text{area}^{4}\text{s}/\text{photon}^{2}$, and is
typically written using Goeppert Mayer units where
$1\text{GM}=10^{-58}\text{m}^{4}\text{s}/\text{photon}^{2}$. Then the rate of
two-photon absorption is given by

$\begin{split}\frac{\mathcal{T}_{\text{coh}}}{T_{p}}&\simeq
I_{\text{coh}}^{2}\sigma_{\text{coh}}^{\text{2PA}}\\\
&=2\hbar\bar{\omega}F_{\text{coh}}^{2}\acute{\sigma}_{\text{coh}}^{\text{2PA}}\\\
\end{split}$ (X.22)

where $2\hbar\omega$ is the energy of the two-photon transition and

$\acute{\sigma}_{\text{coh}}^{\text{2PA}}\equiv\frac{1}{2}\frac{|d_{fe}d_{eg}|^{2}}{(\epsilon_{0}c)^{2}}\frac{(\hbar\bar{\omega})^{2}}{\delta^{2}}\frac{1}{\gamma_{fg}},$
(X.23)

in agreement with past calculations [49, 32].

###  X.2 Choice of parameters

To estimate these cross-sections, and more generally the behavior of the
absorption terms (X.13) for different detunings, we must choose parameter
values in a large parameter space. Our goal is to choose values that can be
used to compare the results for the absorption of coherent light with the
absorption of squeezed light, to be considered in the next section. Since much
of the recent interest in two-photon absorption of squeezed light has been
generated by the possibility of increasing the efficiency of various
applications that use large fluorescent molecules, we focus on parameters
relevant there. However, the appropriate parameters for atoms or small
molecules may be very different.

Consider first $\gamma_{fg}$, where formally in our model the value is set by
the contributions of dephasing and non-radiative population decay. However, as
discussed earlier, we also want to use $\gamma_{fg}$ to phenomenologically
model the decay from the final state to the ground state due to fluorescence.
So we set $\gamma_{fg}/2\pi=1\text{THz}$, which is on the order of previously
discussed values for typical fluorescent dyes [29]. For an order of magnitude
estimate we set $d_{eg}=d_{fe}=|e|a_{0}/\hbar$, where $e$ is the electronic
charge and $a_{0}$ is the Bohr radius. Then using
$\bar{\omega}/2\pi=357\text{THz}$ ($\bar{\lambda}=840\text{nm}$), the adopted
value for $\gamma_{fg}$ and setting $\delta=\bar{\omega}/4$ in Eq. (X.23), we
calculate a two-photon absorption cross-section of
$\acute{\sigma}^{\text{2PA}}_{\text{coh}}=846\text{GM}$. This is larger, but
on the order of the experimentally determined value for Rhodamine B given by
$\acute{\sigma}^{\text{2PA,exp}}_{\text{coh}}=210\pm 55\text{GM}$, which was
excited with a single-mode CW Ti:sapphire laser with an excitation power
$P_{\text{coh}}=1\text{W}$ at $\bar{\lambda}=840\text{nm}$ focused to a beam
area $A\simeq 2.24\times 10^{-13}\text{m}^{2}$, corresponding to the intensity
$I_{\text{coh}}\simeq 4.4\times 10^{12}\text{W}/\text{m}^{2}$ [55].

To set an order of magnitude estimate of $\gamma_{eg}$ we note that for a
clear signal of two-photon absorption we require
$\mathcal{T}_{\text{coh}}\gg\mathcal{O}_{\text{coh}}$, in the far detuned
limit. There the absorption terms in Eq. (X.15) are valid, and this
restriction simplifies to

$\gamma_{eg}\ll
I_{\text{coh}}\frac{|d_{fe}|^{2}}{\epsilon_{0}c}\frac{1}{\gamma_{fg}}.$ (X.24)

For the values given above this requires $\gamma_{eg}/2\pi\ll 0.17\text{THz}$,
which we use to set the order of magnitude estimate to
$\gamma_{eg}/2\pi=0.01\text{THz}$. As a simple example we take
$\gamma_{eg}\simeq\bar{\Gamma}_{eg}$, that is, the leading contribution to
$\gamma_{eg}$ is due to $\bar{\Gamma}_{eg}$ and not the dephasing decay rate
$\Lambda_{eg}$. Within this limit the decay of state $e$ is dominated by
population decay and $R_{ee}(\omega)\simeq 1$.

Finally, we must set an order of magnitude for $\gamma_{fe}$. We begin by
taking each dephasing decay rate to be approximately equal
$\Lambda_{fe}\simeq\Lambda_{fg}\simeq\Lambda_{eg}\ll\bar{\Gamma}_{eg}$, then
each $\gamma_{ij}\simeq\bar{\Gamma}_{ij}$ given by Eq. (VI.8). Under this
assumption $\gamma_{fg}\simeq\bar{\Gamma}_{fg}$ and
$\gamma_{fe}\simeq\bar{\Gamma}_{eg}+\bar{\Gamma}_{fg}$. Since
$\bar{\Gamma}_{fg}\gg\bar{\Gamma}_{eg}$ which were set above, it necessarily
follows from the definitions in Eq. (VI.8) that
$\gamma_{fe}\simeq\bar{\Gamma}_{fg}$ and we set it to be
$\gamma_{fe}/2\pi=1\text{THz}$.

###  X.3 Absorption of coherent light

With these parameters set we can explore how the absorption varies with the
intensity and detunings using the full expressions for each absorption term in
Eq. (X.13). As above we set the incident field to be resonant with the two-
photon transition, and so there is only one detuning parameter $\delta$, but
we allow it to range from below $\gamma_{eg}$ to above the value of
$\gamma_{fe}$. As mentioned earlier, the signs of the absorption terms
$\mathcal{A}^{3}_{\text{coh}}$ and $\mathcal{A}^{5}_{\text{coh}}$ depend on
the linewidth parameters and detunings. For the incident field resonant with
the two-photon transition ($2\bar{\omega}=\omega_{fg}$), as we assume here,
$\mathcal{A}^{5}_{\text{coh}}$ is strictly positive, but the sign of
$\mathcal{A}^{3}_{\text{coh}}$ depends on the detuning $\delta$; when
$\delta<\gamma_{eg}$ it is negative, and when $\delta>\gamma_{eg}$ it is
positive. For small detuning from intermediate state resonances a simple
separation of the absorption into one- and two-photon contributions breaks
down, and it appears to only make sense physically to talk about the total
absorption.

Nonetheless, for our model parameters and intensities where the perturbation
approach is valid (Eq.X.18) we find
$|\mathcal{A}^{3}_{\text{coh}}|\ll|\mathcal{A}^{1}_{\text{coh}}|$ and
$|\mathcal{A}^{3}_{\text{coh}}|\ll|\mathcal{A}^{2}_{\text{coh}}|$ in the
region where $\mathcal{A}^{3}_{\text{coh}}$ is negative, and thus there it is
negligible compared to other contributions to the absorption. So in the region
where $\mathcal{A}^{3}_{\text{coh}}$ is negative the absorption is dominated
by one-photon absorption whether we group $\mathcal{A}^{3}_{\text{coh}}$ as a
one- or two-photon effect. In the far detuned limit, again given by Eq.
(X.15), we find $\mathcal{A}^{3}_{\text{coh}}=\mathcal{A}^{5}_{\text{coh}}$,
and so in that region it clearly contributes to the two-photon absorption.
Therefore for the parameters we are using, it is physically reasonable to
define the general rate of one- and two-photon absorption by

$\displaystyle\frac{\mathcal{O}_{\text{coh}}}{T_{p}}\equiv\frac{\mathcal{A}^{1}_{\text{coh}}+\mathcal{A}^{2}_{\text{coh}}}{T_{p}},$
(X.25a)
$\displaystyle\frac{\mathcal{T}_{\text{coh}}}{T_{p}}\equiv\frac{\mathcal{A}^{3}_{\text{coh}}+\mathcal{A}^{4}_{\text{coh}}+\mathcal{A}^{5}_{\text{coh}}}{T_{p}},$
(X.25b)

which are plotted in Fig. 5 and 6 respectively, where the white areas denote
the parameter region where perturbation theory is no longer valid due to Eq.
(X.18) not being satisfied. In the far detuned and high intensity limit we
find the absorption is dominated by $\mathcal{T}_{\text{coh}}$, as expected.

Figure 5: Plot of the rate of one-photon absorption for coherent light
($\mathcal{O}_{\text{coh}}/T_{p}$) against the detuning from resonance
($\delta$) and intensity $I_{\text{coh}}$. The “white” area represents the
parameter space where the perturbation theory is not valid. Figure 6: Plot of
the rate of two-photon absorption for coherent light
($\mathcal{T}_{\text{coh}}/T_{p}$) against the detuning from resonance
($\delta$) and intensity $I_{\text{coh}}$. The “white” area represents the
parameter space where the perturbation theory is not valid.

##  XI Squeezed Light

Next we consider degenerate squeezed light, which we model with the state

$\ket{\beta}=S(\beta)\ket{\text{vac}},$ (XI.1)

where the unitary squeezing operator is given by

$S(\beta)\equiv e^{\frac{\beta}{2}\int
d\omega_{1}d\omega_{2}\gamma(\omega_{1},\omega_{2})a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})-\text{H.c.}},$
(XI.2)

with $\beta=|\beta|e^{i\theta}$ the squeezing parameter. The function
$\gamma(\omega_{1},\omega_{2})$ is the joint spectral amplitude (JSA), which
satisfies

$\int d\omega_{1}d\omega_{2}|\gamma(\omega_{1},\omega_{2})|^{2}=1.$ (XI.3)

For a typical spontaneous parametric downconversion process, neglecting
dispersion, the form of the JSA in time and frequency is

$\displaystyle\gamma(t_{1},t_{2})=\alpha\left(\frac{t_{1}+t_{2}}{2}\right)\phi(t_{1}-t_{2}),$
(XI.4a)
$\displaystyle\gamma(\omega_{1},\omega_{2})=\alpha(\omega_{1}+\omega_{2})\phi\left(\frac{\omega_{1}-\omega_{2}}{2}\right).$
(XI.4b)

The function $\alpha(\omega)$ is related to the spectral distribution of the
pump used to generate the squeezed light, and $\phi(\omega)$ is the phase-
matching function describing the dispersion properties of the medium mediating
the generation [56].

For the simplest model of squeezed light, where the squeezing only involves a
single (“super-”) mode of the electromagnetic field and the JSA is
unentangled, the one- and two-photon correlation functions are straightforward
to evaluate using textbook results. It is well known that in this limit the
two-photon correlation functions scale linearly and quadratically with photon
number. For calculations where the spectral and temporal degrees of freedom
are taken into account the same scaling is found [1, 2, 3, 13, 2, 14, 15, 16,
29, 30, 31, 32, 33, 15]. The resulting linear dependence in the two-photon
absorption is understood as “photon pair absorption,” since a squeezed state
is a superposition of states involving only pairs of photons. This has been
the center of experimental and theoretical discussions.

For a general squeezed state characterized by a correlated JSA, one cannot
directly calculate the two-photon correlation functions because a
transformation property of the form $S^{\dagger}(\beta)a(\omega)S(\beta)$
cannot be manipulated into a useful form. There are then three approaches that
are immediately suggested.

The first is to restrict oneself to the low flux or “isolated pair” limit. In
this limit we can approximate the state in Eq. (XI.1) as a state that is
mostly the vacuum state but includes a two-photon state with a small
amplitude. This method was considered in great detail by Raymer et al. [29,
30, 31, 32] with the conclusion that, at low photon fluxes and current
technologies, a linear scaling of two-photon absorption should be
undetectable.

A second approach involves performing a Schmidt decomposition of the JSA,
where many Schmidt modes are present; indeed, if one approximates the JSA as a
double Gaussian, the Schmidt decomposition has an analytic result [15]. This
allows one in principle to evaluate the correlation functions, but in practice
one is left with an infinite sum over Schmidt modes that are all delocalized
in frequency.

A third approach is to directly work with the field operators derived from the
nonlinear light generation; this was employed by Dayan [13], and more recently
by Raymer and Landes [33] with similar results where they treated the narrow
pump limit and included dispersion due to the nonlinear light generation, but
with no resonant excitation of intermediate levels.

We employ a different approach here that is most similar to earlier work by
Dayan and Raymer and Landes [13, 33], where we take the squeezed light to be
generated by a CW source. In this limit the spectral correlations are enhanced
to such a strong degree the complexity of the JSA is _reduced_ , and one can
carefully derive a transformation property for
$S^{\dagger}(\beta)a(\omega)S(\beta)$, allowing for the evaluation of the
correlation functions. We can then provide a calculation of the one- and two-
photon absorption for a highly correlated squeezed state in the CW limit, and
in both low and high intensity regimes.

###  XI.1 CW squeezed light

To model the JSA in the CW limit we take

$\displaystyle\alpha(t)=\frac{e^{-i2\bar{\omega}t}}{\sqrt{T_{p}}},\hskip
14.22636pt-\frac{T_{p}}{2}\leq t\leq\frac{T_{p}}{2},$ (XI.5a)
$\displaystyle\alpha(\omega)=\frac{1}{\sqrt{\Omega_{p}}}\operatorname{sinc}\left(\frac{(\omega-2\bar{\omega})\pi}{\Omega_{p}}\right),$
(XI.5b)

which satisfies

$\int dt|\alpha(t)|^{2}=\int d\omega|\alpha(\omega)|^{2}=1,$ (XI.6)

where $T_{p}$ is the length of the pump pulse and $\Omega_{p}=2\pi/T_{p}$
identifies its bandwidth. The pump frequency is centered at $2\bar{\omega}$ so
that the photon pairs generated will be centered at $\bar{\omega}$. To take
the CW limit we let $T_{p}\to\infty$, $\Omega_{p}\to 0$, so that

$\lim_{\Omega_{p}\to 0}|\alpha(\omega)|^{2}=\delta(\omega-2\bar{\omega}).$
(XI.7)

Although general phase-matching functions could be considered, to allow us to
extract analytic results we neglect dispersion and adopt a simple model where
$\phi(\omega)$ takes on a fixed value at the frequencies where it is nonzero:

$\displaystyle\phi(t)=\frac{1}{\sqrt{T_{e}}}\operatorname{sinc}\left(\frac{\pi
t}{T_{e}}\right),$ (XI.8a)
$\displaystyle\phi(\omega)=\frac{1}{\sqrt{\Omega_{e}}},\hskip
14.22636pt-\frac{\Omega_{e}}{2}\leq\omega\leq\frac{\Omega_{e}}{2},$ (XI.8b)

which are both real functions and satisfy

$\int d\omega|\phi(\omega)|^{2}=\int dt|\phi(t)|^{2}=1.$ (XI.9)

The phase-matching function is characterized by the bandwidth $\Omega_{e}$,
which is the range of frequencies over which the generation is effective: The
frequency components of the generated photons are centered at $\bar{\omega}$,
but range between
$\bar{\omega}-\Omega_{e}/2\leq\omega\leq\bar{\omega}+\Omega_{e}/2$. The
related time parameter is $T_{e}=2\pi/\Omega_{e}$ and is called the
“entanglement time”[52, 2, 4, 19, 24, 26, 29, 30, 1, 14]. The quantity $T_{e}$
is a measure of the “coherence time” of pairs of photons, the time over which
a pair of photons can be thought of as being correlated. Indeed, photons are
generated at times $t_{1}$ and $t_{2}$ according to the JSA (Eq. (XI.4a))
which is proportional to $\phi(t_{1}-t_{2})$, which from Eq. (XI.8a) is peaked
for times $t_{1}-t_{2}<T_{e}$. The form of the JSA in time and frequency
variables is schematically shown in Fig. 7. For a highly correlated state the
parameters satisfy $T_{p}\gg T_{e}$ or equivalently $\Omega_{e}\gg\Omega_{p}$;
we consider this limit in this paper.

Figure 7: Schematic of JSA for the pump and phase matching function given in
Eq. (XI.5) and (XI.8). Our diagram is a schematic representation because we
have neglected to include the sinc oscillations present, however, in the CW
limit when the state is strongly correlated the schematic representation is a
good approximation of the JSA being considered.

In Appendix E we work out the operator $S^{\dagger}(\beta)a(\omega)S(\beta)$
for squeezed light from a medium pumped by a CW source with a general phase-
matching function $\phi(\omega)$ (Eq.(E.15)). Specifying to our simple model
(XI.8) of the phase-matching function, that transformation simplifies to

$\begin{split}&S^{\dagger}(\beta)a(\omega)S(\beta)\\\ &\to
a(\omega)+\frac{[c(\omega)-1]}{\sqrt{\Omega_{p}}}\int
d\omega^{\prime}\alpha\left(2\bar{\omega}-\omega+\omega^{\prime}\right)a(\omega_{2})\\\
&+\frac{s(\omega)e^{i\theta}}{\sqrt{\Omega_{p}}}\int
d\omega^{\prime}\alpha\left(\omega+\omega_{1}\right)a^{\dagger}(\omega^{\prime}),\end{split}$
(XI.10)

where we have set

$\displaystyle s(\omega)=s=\text{sinh}(|\beta_{0}|)\hskip
5.69054pt\text{if}\hskip
5.69054pt|\omega-\bar{\omega}|\leq\frac{\Omega_{e}}{2},$ (XI.11a)
$\displaystyle c(\omega)=c=\text{cosh}(|\beta_{0}|)\hskip
5.69054pt\text{if}\hskip
5.69054pt|\omega-\bar{\omega}|\leq\frac{\Omega_{e}}{2},$ (XI.11b)
$\displaystyle s(\omega)=0\hskip 5.69054pt\text{if}\hskip
5.69054pt|\omega-\bar{\omega}|>\frac{\Omega_{e}}{2},$ (XI.11c) $\displaystyle
c(\omega)=1\hskip 5.69054pt\text{if}\hskip
5.69054pt|\omega-\bar{\omega}|>\frac{\Omega_{e}}{2},$ (XI.11d)

which satisfy $s(2\bar{\omega}-\omega)=s(\omega)$ and
$c(2\bar{\omega}-\omega)=c(\omega)$ and define the squeezing parameter
$\beta_{0}$ by setting $\beta_{0}\equiv\beta\sqrt{T_{e}/T_{p}}$. As we will
see, the introduction of $\beta_{0}$ essentially “normalizes” the squeezing
parameter in a way analogous to the introduction of $\alpha_{0}$ from $\alpha$
in Sec. X. If we keep the intensity in the center of the pulse the same, the
squeezing parameter $\beta$ will be proportional to $\sqrt{T_{p}/T_{e}}$,
because $\beta$ is related to the average photon number. From the above
definition of $\beta_{0}$ one can see that $\beta_{0}$ will remain fixed and
independent of $T_{p}$ and $T_{e}$; although $\beta$ diverges in the CW limit,
$\beta_{0}$ is fixed.

With the transformation property of Eq. (LABEL:eq:squeezed_state_transform_v2)
we can calculate the one- and two-photon correlation functions (Appendix F).
The one-photon correlation function evaluated at two different frequencies
(Eq. (F.6)) is given by

$\begin{split}\bra{\beta}a^{\dagger}(\omega_{2})a(\omega_{1})\ket{\beta}=\frac{T_{p}}{2\pi}s(\omega_{2})s(\omega_{1})\text{sinc}\left(\frac{(\omega_{1}-\omega_{2})\pi}{\Omega_{p}}\right)\end{split}$
(XI.12)

which is proportional to $T_{p}$ as would be expected. Note that the factor
$s(\omega_{2})s(\omega_{1})$ includes the restriction of frequencies to be
within the bandwidth of the squeezed light. Using Eq. (XI.12) we can
immediately calculate the pulse energy, which is

$E_{\text{sq}}=\hbar\bar{\omega}\frac{T_{p}}{T_{e}}\text{sinh}^{2}\left(|\beta_{0}|\right).$
(XI.13)

We can identify the flux of photons to be

$F_{\text{sq}}\equiv\frac{1}{AT_{e}}\text{sinh}^{2}\left(|\beta_{0}|\right)$
(XI.14)

(and the total number of photons in the pulse to be
$N_{\text{sq}}=T_{p}AF_{\text{sq}}$), for then we obtain the correct
expression for the intensity of light in the first line below,

$\begin{split}I_{\text{sq}}&\equiv\frac{E_{\text{sq}}}{AT_{p}}=\hbar\bar{\omega}F_{\text{sq}}\\\
&=\frac{\hbar\bar{\omega}}{AT_{e}}\text{sinh}^{2}\left(|\beta_{0}|\right).\end{split}$
(XI.15)

while from the second line we can identify the number of photons within the
entanglement time $T_{e}$ to be
$N_{T_{e}}\equiv\text{sinh}^{2}\left(|\beta_{0}|\right)$.

The squeezed state can then be understood as follows: Photons pairs are
created throughout the length of the CW beam ($-T_{p}/2\leq t\leq T_{p}/2$)
with pairs being localized within a time $T_{e}$;
$N_{T_{e}}=\text{sinh}^{2}\left(|\beta_{0}|\right)$ is the expected number of
photons within a time $T_{e}$; and since $\beta_{0}$ is fixed, $N_{T_{e}}$
does not depend on $T_{e}$ or $T_{p}$.

To recover the weakly squeezed (isolated pair) limit [31], we take
$|\beta_{0}|\ll 1$ so that we can approximate $N_{T_{e}}\to|\beta_{0}|^{2}\ll
1$. Then we identify $|\beta_{0}|^{2}$ as the number of photons within a time
$T_{e}$; and if $\beta_{0}$ is small enough, we expect there to be at most one
pair of photons within a length of time $T_{e}$. Then the flux is
$F_{\text{sq}}\to|\beta_{0}|^{2}/(AT_{e})$ and the number of photons is given
by $N_{\text{sq}}\to T_{p}|\beta_{0}|^{2}/T_{e}=|\beta|^{2}$, acknowledging
$|\beta|^{2}$ to be the total number of photons in this limit.

Moving to the two-photon correlation function (Eq.
(LABEL:eq:squeezed_state_two-photon_correlation_function)), we have

$\begin{split}&\bra{\beta}a^{\dagger}(\omega_{4})a^{\dagger}(\omega_{3})a(\omega_{2})a(\omega_{1})\ket{\beta}\\\
&=\frac{T_{p}}{2\pi}s(\omega_{3})c(\omega_{3})c(\omega_{1})s(\omega_{1})\alpha(\omega_{1}+\omega_{2})\alpha(\omega_{3}+\omega_{4})\\\
&+\frac{T_{p}}{2\pi}s^{2}(\omega_{2})s^{2}(\omega_{1})\alpha(\omega_{1}-\omega_{4}+2\bar{\omega})\alpha(\omega_{2}-\omega_{3}+2\bar{\omega})\\\
&+\frac{T_{p}}{2\pi}s^{2}(\omega_{2})s^{2}(\omega_{1})\alpha(\omega_{1}-\omega_{3}+2\bar{\omega})\alpha(\omega_{2}-\omega_{4}+2\bar{\omega}),\\\
\end{split}$ (XI.16)

again proportional to $T_{p}$ as expected. The two-photon correlation function
has three contributions, and in each the factors $s(\omega)$ ensure that the
photons are within the bandwidth of the squeezed light. The first term,
previously referred to as the “coherent” contribution [13, 33] is due to
photons that are anti-correlated such that
$\omega_{1}+\omega_{2}=2\bar{\omega}$. Due to the anti-correlations this
contribution is dominant for $\mathcal{A}^{3}_{\text{sq}}$ and
$\mathcal{A}^{5}_{\text{sq}}$, which are derived from absorption pathways
which include a two-photon transition with only a virtual excitation of the
intermediate level.

The second and third contributions to Eq.
(LABEL:eq:squeezed_correlation_fn_2_v1), previously referred to as the
“incoherent contribution [13, 33] are correlated such that
$\omega_{1}=\omega_{4}$ and $\omega_{2}=\omega_{3}$ for the second term and
$\omega_{1}=\omega_{3}$ and $\omega_{2}=\omega_{4}$ for the third term.
Although the third contribution is very similar to the second, the different
frequency correlations can lead to large differences in how they contribute to
the absorption because the broadening function
$\mathcal{R}^{i}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ (see Eq.
(IX.3)) is in general not symmetric under exchange of any of its frequency
arguments. Specifically, we expect the second and third term to dominate for
$\mathcal{A}^{2}_{\text{sq}}$ and $\mathcal{A}^{4}_{\text{sq}}$. But because
the “rephasing” and “nonrephasing” pathway is doubly resonant when
$\omega_{2}=\omega_{3}$, which can be seen by examining
$\mathcal{R}^{2}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ and
$\mathcal{R}^{4}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ given in Eq.
(IX.3), we expect the second contribution to be larger than the third.

Using the model we have developed and the correlation functions we derived in
Eq. (XI.12) and (LABEL:eq:squeezed_correlation_fn_2_v1) we can immediately
calculate the normalized second order correlation function in frequency,
defined by

$g^{(2)}(\omega_{1},\omega_{2})\equiv\frac{\bra{\beta}a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})a(\omega_{2})a(\omega_{1})\ket{\beta}}{\bra{\beta}a^{\dagger}(\omega_{2})a(\omega_{2})\ket{\beta}\bra{\beta}a^{\dagger}(\omega_{1})a(\omega_{1})\ket{\beta}}.$
(XI.17)

We find

$\begin{split}g^{(2)}(\omega_{1},\omega_{2})&=\frac{c^{2}(\omega_{1})}{s^{2}(\omega_{1})}\text{sinc}^{2}\left(\frac{(\omega_{1}+\omega_{2}-2\bar{\omega})\pi}{\Omega_{p}}\right)\\\
&+1+\text{sinc}^{2}\left(\frac{(\omega_{1}-\omega_{2})\pi}{\Omega_{p}}\right)\end{split}$
(XI.18)

which is in agreement with experiment when dispersion is compensated [57]. The
anti-correlated term in which $\omega_{1}+\omega_{2}=2\bar{\omega}$ is the
“cross-correlation” contribution and the correlation along the diagonal where
$\omega_{1}=\omega_{2}$ is the “auto-correlation” of each frequency mode with
itself for indistinguishable photons. Here each correlation is sensitive to
the pump bandwidth set by $\Omega_{p}$.

Alternatively we calculate the normalized second order correlation function in
time, defined by

$g^{(2)}(t_{1},t_{2})\equiv\frac{\langle\beta|a^{\dagger}(t_{1})a^{\dagger}(t_{2})a(t_{2})a(t_{1})|\beta\rangle}{\langle\beta|a^{\dagger}(t_{2})a(t_{2})|\beta\rangle\langle\beta|a^{\dagger}(t_{1})a(t_{1})|\beta\rangle},$
(XI.19)

which we calculate by taking the Fourier transforms of Eq. (XI.12) and
(LABEL:eq:squeezed_correlation_fn_2_v1) to give

$\begin{split}g^{(2)}(t_{1},t_{2})=1+\left(2+\frac{1}{N_{T_{e}}}\right)\text{sinc}^{2}\left(\frac{(t_{1}-t_{2})\pi}{T_{e}}\right)\end{split}$
(XI.20)

which yields the usual result for $t_{1}=t_{2}$ [58, 59, 33, 60, 37] and
satisfies the inequality that $g^{(2)}(0,\tau)\leq g^{(2)}(0,0)$, the
condition for “bunched light”. However, we note that while in frequency space
the auto-correlation term (third contribution of Eq. (XI.18)) is peaked when
$\omega_{1}=\omega_{2}$ and is characterized by the pump width $\Omega_{p}$,
in time this contribution lead to a correlation that $t_{1}=t_{2}$ on time
scales _not_ given by the inverse bandwidth of the pump but by the
entanglement time $T_{e}$.

Using the correlation functions (XI.12), the lowest order rate of absorption
is given by

$\frac{\mathcal{A}^{1}_{\text{sq}}}{T_{p}}=I_{\text{sq}}\sigma^{1}_{\text{sq}},\\\
$ (XI.21)

which scales linearly with intensity as expected, but with a modified cross-
section due to the large bandwidth of the squeezed light given by

$\sigma^{\text{1}}_{\text{sq}}=\int\limits_{\bar{\omega}-\frac{\Omega_{e}}{2}}^{\bar{\omega}+\frac{\Omega_{e}}{2}}\frac{d\omega}{\Omega_{e}}\text{Im}\left[\mathcal{R}^{1}(\omega)\right].$
(XI.22)

For each higher order absorption term we use the correlation function
(LABEL:eq:squeezed_correlation_fn_2_v1), together with the property that in
the CW limit $\alpha(\omega)$ is strongly peaked at $\omega=2\bar{\omega}$,
then the rate of each higher order absorption term is given by

$\frac{\mathcal{A}^{i}_{\text{sq}}}{T_{p}}=I_{\text{sq}}(I_{\text{vac}}+I_{\text{sq}})\sigma^{i}_{\text{sq,{I}}}+2I_{\text{sq}}^{2}\sigma^{i}_{\text{sq,{II}}},$
(XI.23)

where we find it useful to define a quantity called the “vacuum intensity” by

$I_{\text{vac}}=\frac{\hbar\bar{\omega}}{AT_{e}},$ (XI.24)

following the work of Klyshko [9] where a similar quantity was defined, which
we understand as the intensity of one-photon passing through an area $A$
within a time $T_{e}$. We find each higher order absorption term is dependent
on two cross-sections $\sigma^{i}_{\text{sq,{I}}}$ and
$\sigma^{i}_{\text{sq,{II}}}$ for $i=2,3,4,5$. The first cross-section
$\sigma^{i}_{\text{sq,{I}}}$ is given by

$\sigma^{i}_{\text{sq,{I}}}=\int\limits_{\bar{\omega}-\frac{\Omega_{e}}{2}}^{\bar{\omega}+\frac{\Omega_{e}}{2}}\frac{d\omega_{1}d\omega_{2}}{\Omega_{e}^{2}}\text{Im}\left[\mathcal{R}^{i}(\omega_{1},2\bar{\omega}-\omega_{1},\omega_{2},2\bar{\omega}-\omega_{2})\right],$
(XI.25)

where we use the subscript “I” for $\sigma^{i}_{\text{sq,{I}}}$ to denote that
this cross-section corresponds to the anti-correlated term of the correlation
function (first term in Eq. (LABEL:eq:squeezed_correlation_fn_2_v1)). The
second cross-section $\sigma^{i}_{\text{sq,{II}}}$ is given by

$\begin{split}\sigma^{i}_{\text{sq,{II}}}&=\frac{1}{2}\int\limits_{\bar{\omega}-\frac{\Omega_{e}}{2}}^{\bar{\omega}+\frac{\Omega_{e}}{2}}\frac{d\omega_{1}d\omega_{2}}{\Omega_{e}^{2}}\text{Im}\left[\mathcal{R}^{i}(\omega_{1},\omega_{2},\omega_{2},\omega_{1})\right.\\\
&\hskip
113.81102pt+\left.\mathcal{R}^{i}(\omega_{1},\omega_{2},\omega_{1},\omega_{2})\right],\end{split}$
(XI.26)

where we use the subscript “II” for $\sigma^{i}_{\text{sq,{II}}}$ to denote
that this cross-section corresponds to the second and third term from the
correlation function given by Eq. (LABEL:eq:squeezed_correlation_fn_2_v1).
Here we combined the second and third contribution into a single term and
pulled out a factor of $2$, so this can be understood as an average cross-
section.

Considering the form of each higher order absorption term
$\mathcal{A}^{i}_{\text{sq}}$ in Eq. (XI.23) we find two contributions. The
first, which is proportional to $I_{\text{sq}}(I_{\text{vac}}+I_{\text{sq}})$,
is multiplied by the cross-section $\sigma^{i}_{\text{sq,{I}}}$ corresponding
to photons that are anti-correlated in frequency. For weak squeezing such that
$\beta_{0}\ll 1$, then $N_{T_{e}}\ll 1$ and we can think of only a single pair
of photons within a time $T_{e}$ leading to a linear scaling of absorption
with intensity and $I_{\text{vac}}$ playing the role of an enhancement factor.
Then as $\beta_{0}$ increases, $N_{T_{e}}\gg 1$ where there are many photons
within a time $T_{e}$ so the absorption scales quadratically with intensity
but is still multiplied by $\sigma^{i}_{\text{sq,{I}}}$.

While one view is that squeezed light modifies the cross-section, and so an
“entangled cross-section” is sometimes defined [52, 4, 18, 19, 24, 26, 27, 14,
29, 30, 31], and in our notation is given by
$I_{\text{vac}}\sigma^{i}_{\text{sq,{I}}}$, here we take a different approach
because in the current form it will be more straight forward to compare to
coherent light. Further, the contribution from $\sigma^{i}_{\text{sq,{I}}}$
involves not only a linear scaling with intensity, but also a quadratic one.

The second contribution on the right-hand side of Eq. (XI.23), which is due to
$\sigma^{i}_{\text{sq,{II}}}$, scales as intensity squared. The result of
averaging the contributions to the cross-section $\sigma^{i}_{\text{sq,{II}}}$
in Eq. (XI.26) is a factor of 2, for the two ways of choosing two
indistinguishable photons.

To evaluate the squeezed state cross-sections and absorption terms in general
we need to resort to numerical integration. However, we first consider the
limiting case where the squeezed state has a narrow bandwidth $\Omega_{e}$.

####  XI.1.1 Squeezed state absorption: narrow bandwidth limit

From Eq. (XI.22) and (XI.23) we see that to evaluate the cross-sections for
absorption we must integrate over the bandwidth of squeezed light. The
integrals can be simplified if we take the limit when the bandwidth
$\Omega_{e}$ of the squeezed light is small compared to the decay widths
$\gamma_{ij}$. In this limit, the coherence time of the photon pairs is long
compared to the decay times of the molecule but each photon has a well defined
frequency. Since each photon has a well defined frequency, the spectral
correlations within $\sigma^{i}_{\text{sq,{I}}}$ and
$\sigma^{i}_{\text{sq,{II}}}$ ((XI.25) and (XI.26)) are no longer relevant. To
good approximation we can evaluate each cross-section in Eq. (XI.22), (XI.25)
and (XI.26) at the centre frequency $\bar{\omega}$ and we find that

$\displaystyle\left.\sigma_{\text{sq}}^{1}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=\sigma_{\text{coh}}^{1}$ (XI.27a)
$\displaystyle\left.\sigma^{i}_{\text{sq,{I}}},\sigma^{i}_{\text{sq,{II}}}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=\sigma_{\text{coh}}^{i}.$ (XI.27b)

Here each squeezed light cross-section reduces to the coherent light cross-
sections we defined in equation (X.12).

The one- and two-photon absorption rates then simplify to

$\displaystyle\left.\frac{\mathcal{A}^{1}_{\text{sq}}}{T_{p}}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=I_{\text{sq}}\sigma^{1}_{\text{coh}},$ (XI.28a)
$\displaystyle\left.\frac{\mathcal{A}^{i}_{\text{sq}}}{T_{p}}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=(3I_{\text{sq}}^{2}+I_{\text{sq}}I_{\text{vac}})\sigma^{i}_{\text{coh}},$
(XI.28b)

where we expect the second contribution in Eq. (XI.28b) given by
$I_{\text{sq}}I_{\text{vac}}\sigma^{i}_{\text{coh}}$ to be small. In the limit
we are considering, $T_{e}$ is much longer than decay times of the molecule so
we expect the contribution from photon pair absorption to be negligible.

In this limit we have

$\displaystyle\left.\frac{\mathcal{A}^{1}_{\text{sq}}}{\mathcal{A}^{1}_{\text{coh}}}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=1,$ (XI.29a)
$\displaystyle\left.\frac{\mathcal{A}^{i}_{\text{sq}}}{\mathcal{A}^{i}_{\text{coh}}}\right|_{\frac{\Omega_{e}}{\gamma_{ij}}\ll
1}=3+\frac{I_{\text{vac}}}{I_{\text{sq}}}=g^{(2)}(0,0),$ (XI.29b)

where to make the comparison we took $I_{\text{coh}}=I_{\text{sq}}$; these
results are in agreement with Raymer and Landes [33]. We find the lowest order
contribution is identical in this limit, which is expected since each photon
from either coherent or squeezed light is at $\bar{\omega}$ and pair
correlations are irrelevant. For each higher order term we find that for
squeezed light there is a strong enhancement when $I_{\text{vac}}\gg
I_{\text{sq}}$. Taking the model parameters from section X.2 and setting
$\Omega_{e}/2\pi=10^{-4}\text{THz}\ll\gamma_{ij}$ ($T_{e}=10^{7}\text{fs}$),
we calculate $I_{\text{vac}}=105\text{W}/\text{m}^{2}$. Thus _very_ low
intensities would be required to achieve a significant enhancement by using
squeezed light rather than coherent light.

Therefore, in this limit for reasonable intensity regimes where two-photon
absorption experiments are done, $I_{\text{vac}}/I_{\text{sq}}\ll 1$ so
squeezed states provide at most a factor of $3$ enhancement for each
absorption term due to photon bunching.

For a complete characterization of the absorption beyond the narrow bandwidth
limit, which is valid in both the low and high intensity regimes, we resort to
a numerical analysis.

####  XI.1.2 Squeezed state absorption: numerical analysis

We now calculate the $\sigma_{\text{sq}}^{i}$ cross-sections in Eq. (XI.22),
(XI.25) and (XI.26) numerically using the decay parameters $\gamma_{ij}$
discussed in section X.2 and set $\Omega_{e}/2\pi=10\text{THz}$, which
corresponds to the limit

$\gamma_{eg}\ll\gamma_{fe}=\gamma_{fg}<\Omega_{e}.$ (XI.30)

We set $2\bar{\omega}=\omega_{fg}$, and look at the cross-sections as a
function of detuning $\delta$ from resonance with the intermediate state (See
Fig. 4).

We begin with the $\sigma^{1}$ cross-sections for coherent and squeezed light,
plotted in Fig. 8. The coherent light cross-section is fixed near resonance
($\delta<\gamma_{eg}$) and decreases for large detuning
($\delta>\gamma_{eg}$). Due to the large squeezed light bandwidth, the
squeezed light cross-section is smaller but remains fixed for larger values of
the detuning until the detuning is larger than half the squeezed state
bandwidth ($\delta>\Omega_{e}/2$), at which point there is no resonant
overlap. Thus for detunings $\gamma_{eg}\ll\delta<\Omega_{e}/2$, there is
approximately three orders of magnitude difference between the cross-section
for squeezed light and that for coherent light, due simply to the large
bandwidth of squeezed light. In the very far detuned limit
($\delta>\Omega_{e}/2\gg\gamma_{eg}$), $\delta$ is the largest quantity in the
broadening function, so integrating over the squeezed light bandwidth is
negligible compared to $\delta$; therefore the squeezed light cross-section
approaches the coherent limit.

Figure 8: Plotting the $\sigma^{1}$ cross-sections for coherent and squeezed
light vs detuning from resonance $\delta$.

In Fig. 9 we plot the negative of the coherent and squeezed light $\sigma^{2}$
cross-sections. The coherent light cross-section has the same behaviour as
that of $\sigma^{1}_{\text{coh}}$, fixed near resonance and decreasing far
from resonance. For the squeezed light cross-section we plot both
$\sigma_{\text{sq,{I}}}$ and $\sigma_{\text{sq,{II}}}$; note that
$\sigma^{2}_{\text{sq,{II}}}\geq\sigma^{2}_{\text{sq,{I}}}$. To see why,
consider the broadening function involved in this cross-section,

$\mathcal{R}^{2}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})\propto\frac{R_{ee}(\omega_{2}-\omega_{3})}{Q_{eg}(\omega_{4})Q_{eg}^{*}(\omega_{3})Q_{eg}(\omega_{2})}.$
(XI.31)

For anti-correlated photons such that $\omega_{3}+\omega_{4}=2\bar{\omega}$,
this term is small because the two factors of the product
$Q_{eg}(\omega_{4})Q_{eg}^{*}(\omega_{3})$ cannot be resonant at the same time
for most frequencies within the bandwidth; the peak at
$\delta\simeq\Omega_{e}/2$ is nontrivial, and will be discussed in greater
detail below. However, if we consider photons such that
$\omega_{2}=\omega_{3}$ and $\omega_{1}=\omega_{4}$ that product is doubly
resonant, the broadening function simplifies to a squared Lorentzian, and so
has the same behaviour as $\sigma^{1}_{\text{sq}}$. In both cases, when very
far detuned ($\delta>\Omega_{e}/2$), $\delta$ is the largest quantity and each
cross-section approaches the coherent light value.

Figure 9: Plotting the $\sigma^{2}$ cross-sections for coherent and squeezed
light vs detuning from resonance $\delta$.

Next we consider the $\sigma^{3}$ cross-sections, which has a unique
behaviour. When $\delta<\gamma_{eg}$ the coherent light cross-section is
negative, identically zero when $\delta=\gamma_{eg}$, and positive but
decreasing when $\delta>\gamma_{eg}$. The behaviour is not obvious from the
diagram III in Fig. 2 and the perturbative expansion in Fig. 3.c). But by
setting $\omega_{fg}=2\bar{\omega}$ in Eq. (X.13c) for
$\mathcal{A}^{3}_{\text{coh}}$ one can analytically see this scaling with
$\delta$. For the squeezed light cross-sections we see a similar behavior;
however, the turning point from negative to positive is shifted to
$\delta\simeq\Omega_{e}/2$ due to the large bandwidth. In Fig. 10 we show the
different $\sigma^{3}$ cross-sections for the values of $\delta$ for which
they are positive and negative. In the large detuned limit
($\delta>\Omega_{e}/2$) we find
$\sigma^{3}_{\text{sq,{I}}}>\sigma^{3}_{\text{sq,{II}}}$ because the bandwidth
is larger than the virtual two-photon absorption pathway linewidth, i.e.
$\Omega_{e}>\gamma_{fg}$. Absorption from anti-correlated photons is therefore
always resonant with $\omega_{fg}$, while for uncorrelated photons with
frequency $\omega_{1}$ and $\omega_{2}$ there is wasted energy along the
endpoints of the squeezed light bandwidth. Since $\sigma^{3}_{\text{sq,{I}}}$
is always resonant with the virtual two-photon absorption it approaches the
coherent light limit for large detuning ($\delta>\Omega_{e}/2$).

Figure 10: Plotting the $\sigma^{3}$ cross-section for coherent and squeezed
light vs detuning from resonance $\delta$.

In Fig. 11 we plot the $\sigma^{4}$ cross-sections. The coherent light cross-
section is fixed near resonance and then decreases far from resonance with two
different scalings. When $\delta<\gamma_{eg}$ it scales as $\delta^{-2}$, but
as $\delta$ increases so that $\delta>\gamma_{fe}$ the cross-section scales as
$\delta^{-4}$. The squeezed light cross-sections satisfy
$\sigma^{4}_{\text{sq,{II}}}\geq\sigma^{4}_{\text{sq,{I}}}$, and again this
can be understood from the behaviour of the relevant broadening function

$\mathcal{R}^{4}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})\propto\frac{R_{ee}(\omega_{2}-\omega_{3})}{Q_{fe}(\omega_{4})Q_{eg}^{*}(\omega_{3})Q_{eg}(\omega_{2})}.$
(XI.32)

When the photons frequencies are given by $\omega_{1}=\omega_{4}$ and
$\omega_{3}=\omega_{2}$, then both factors in the product
$Q_{eg}^{*}(\omega_{3})Q_{eg}(\omega_{2})$ become small at the same time. This
is the leading contribution to the squeezed state cross-section when
$\delta<\Omega_{e}/2$. However, as $\delta$ approaches $\Omega_{e}/2$, half
the squeezed light bandwidth, $\sigma^{4}_{\text{sq,{I}}}$ quickly increases
to its maximum value. This is a nontrivial scaling with the detuning which we
discuss in detail below. In the very far detuned limit
($\delta>\Omega_{e}/2$), $\delta$ is the leading quantity and each cross-
section approaches the coherent light cross-section. Similar to the
$\sigma^{1}$ and $\sigma^{2}$ cross-sections, $\sigma^{4}_{\text{sq,{II}}}$ is
approximately three orders of magnitude larger than the coherent light cross-
section due to the large squeezed light bandwidth.

Figure 11: Plotting the $\sigma^{4}$ cross-sections for coherent and squeezed
light vs detuning from resonance $\delta$.

Finally we consider the $\sigma^{5}$ cross-sections, which is plotted in Fig.
12. The coherent light cross-section also has two behaviours between
$\delta<\gamma_{eg}$ and $\delta>\gamma_{fe}$. However unlike
$\sigma^{4}_{\text{coh}}$, they both scale as $\delta^{-2}$ because
$\sigma^{5}_{\text{coh}}$ has a $\delta^{2}$ dependence in the numerator,
which can be seen from Eq. (X.13e). We find
$\sigma^{5}_{\text{sq,{I}}}>\sigma^{5}_{\text{sq,{II}}}$ since the
$\sigma^{5}$ cross-sections corresponds to a virtual two-photon transition and
$\Omega_{e}>\gamma_{fg}$, similar to the $\sigma^{3}$ cross-sections. In the
large detuned limit the anti-correlated cross-section approaches the coherent
limit and we find a similar cusp behaviour as in
$\sigma^{2}_{\text{sq,{II}}},\sigma^{3}_{\text{sq,{I}}},\sigma^{3}_{\text{sq,{II}}}$
and $\sigma^{4}_{\text{sq,{I}}}$ which we now discuss.

Figure 12: Plotting the $\sigma^{5}$ cross-sections for coherent and squeezed
light vs detuning from resonance $\delta$.

In each squeezed light cross-section there is a “cusp” behaviour located at
$\delta\simeq\Omega_{e}/2$. This is due to an asymmetric scaling with
frequency, and is $\delta$ dependent. Although we could choose any broadening
function to illustrate this behavior we consider the broadening function
$\mathcal{R}^{5}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$, which is the
simplest mathematically. Consider
$\mathcal{R}^{5}(\omega_{1},\omega_{2},\omega_{3},\omega_{4})$ evaluated at
the anti-correlated frequencies, given by

$\begin{split}&\mathcal{R}^{5}(\omega_{1},2\bar{\omega}-\omega_{1},\omega_{2},2\bar{\omega}-\omega_{2})\propto\\\
&\hskip
85.35826pt\frac{-1}{Q_{fe}(2\bar{\omega}-\omega_{2})Q_{fg}(2\bar{\omega})Q_{eg}(2\bar{\omega}-\omega_{1})},\end{split}$
(XI.33)

which is resonant with the virtual two-photon transition as expected. Taking
the imaginary part and changing variables $\sigma^{5}_{\text{sq,{I}}}$ is
proportional to

$\sigma^{5}_{\text{sq,{I}}}\propto\frac{1}{\gamma_{fg}}\int\limits_{\delta-\frac{\Omega_{e}}{2}}^{\delta+\frac{\Omega_{e}}{2}}\frac{d\omega_{1}d\omega_{2}}{\Omega_{e}^{2}}\frac{\omega_{1}\omega_{2}+\gamma_{fe}\gamma_{eg}}{[\omega_{1}^{2}+\gamma_{fe}^{2}][\omega_{2}^{2}+\gamma_{eg}^{2}]}.$
(XI.34)

Due to the behaviour of the first term there is complete cancellation when
$\delta=0$ and the contribution increases to its maximum when
$\delta=\Omega_{e}/2$, where the range of integration is over a strictly
positive function. In this case the integral is simple enough to be worked out
analytically and the asymmetric dependence on frequencies in the numerator is
responsible for the quick increase in cross-section values as $\delta$
approaches $\Omega_{e}/2$. For each resonant denominator we must evaluate it
at the specified frequencies and then take the imaginary part. Taking the
imaginary part leads to frequency dependent terms in the numerator which are
asymmetric and lead to the resulting “cusp” behavior.

Following the discussion in Sec. (X.3) we define the total squeezed state two-
photon absorption cross-sections by

$\displaystyle\sigma^{\text{2PA}}_{\text{sq,{I}}}=\sigma^{3}_{\text{sq,{I}}}+\sigma^{4}_{\text{sq,{I}}}+\sigma^{5}_{\text{sq,{I}}},$
(XI.35a)
$\displaystyle\sigma^{\text{2PA}}_{\text{sq,{II}}}=\sigma^{3}_{\text{sq,{II}}}+\sigma^{4}_{\text{sq,{II}}}+\sigma^{5}_{\text{sq,{II}}},$
(XI.35b)

which are plotted in Fig. 13 and compared to the total coherent light two-
photon absorption cross-section, given by

$\sigma^{\text{2PA}}_{\text{coh}}=\sigma^{3}_{\text{coh}}+\sigma^{4}_{\text{coh}}+\sigma^{5}_{\text{coh}},$
(XI.36)

where the $\sigma^{3}$ cross-sections include both positive and negative
contributions (see discussion at the end of section X.3).

The cross-sections in Fig. 13 are the sums of the contributions already
discussed above. The coherent light two-photon absorption cross-section is
fixed for $\delta\ll\gamma_{eg}$ and then increases to its maximum value at
$\delta=\gamma_{eg}$; then it begins to decrease with a scaling of
$\delta^{-2}$. The increase in the total coherent cross-section until its
maximum at $\delta=\gamma_{eg}$ is a direct result of the coherent light
cross-section $\sigma^{3}_{\text{coh}}$ being negative when
$\delta<\gamma_{eg}$. When $\delta<\Omega_{e}/2$
$\sigma^{\text{2PA}}_{\text{sq,{II}}}>\sigma^{\text{2PA}}_{\text{sq,{I}}}$;
which is a direct result of the large resonant enhancement due to the large
bandwidth of the squeezed light. However, as soon as $\delta$ becomes larger
than $\Omega_{e}/2$ all resonant enhancement is lost and
$\sigma^{\text{2PA}}_{\text{sq,{I}}}>\sigma^{\text{2PA}}_{\text{sq,{II}}}$
because $\gamma_{fg}<\Omega_{e}$, and there is wasted energy at the endpoints
of the bandwidth. In the limit where $\gamma_{fg}\ll\Omega_{e}$ this decrease
would be even larger and we would find $\sigma^{\text{2PA}}_{\text{sq,{II}}}$
to be negligible, in agreement with Raymer and Landes [33].

Figure 13: Plotting the $\sigma^{\text{2PA}}$ cross-sections for coherent and
squeezed light vs detuning from resonance $\delta$.

The largest enhancement to the squeezed state two-photon absorption cross-