Two-photon absorption (TPA) and other nonlinear interactions of molecules with time–frequency-entangled photon pairs have been predicted to display a variety of fascinating effects. Therefore, their potential use in practical quantum-enhanced molecular spectroscopy requires close examination. This Tutorial presents a detailed theoretical study of one- and two-photon absorption by molecules, focusing on how to treat the quantum nature of light. We review some basic quantum optics theory and then we review the density-matrix (Liouville) derivation of molecular optical response, emphasizing how to incorporate quantum states of light into the treatment. For illustration, we treat in detail the TPA of photon pairs created by spontaneous parametric down conversion, with an emphasis on how quantum light TPA differs from that with classical light. In particular, we treat the question of how much enhancement of the TPA rate can be achieved using entangled states. This Tutorial includes a review of known theoretical methods and results as well as some extensions, especially the comparison of TPA processes that occur via far-off-resonant intermediate states only and those that involve off-resonant intermediate states by virtue of dephasing processes. A brief discussion of the main challenges facing experimental studies of entangled two-photon absorption is also given.

## I. INTRODUCTION

In the past decade, many theoretical studies have proposed that “quantum advantages” in spectroscopy can be achieved by the use of time–frequency-entangled photon pairs (EPPs). For review, see Refs. 1–4. These include, for example, proposals for virtual-state spectroscopy,^{5,6} Raman spectroscopy,^{7} and multi-dimensional optical spectroscopy.^{8,9} Two-photon absorption (TPA) with entangled light can provide a testbed for many of these ideas, as indicated by a flurry of recent papers addressing the experimental feasibility of using entangled two-photon absorption (ETPA) to increase the measured signal in spectroscopic or imaging applications. This Tutorial presents not a complete review, but a tutorial treatment of the underlying quantum optics needed for understanding how quantum light interacts with molecules, especially one- and two-photon absorption. We review some basic quantum optics theory, and then we review the density-matrix (Liouville) derivation of molecular optical response, emphasizing how to incorporate quantum states of light into the treatment. The treatment focuses on photon pairs at very low flux where different pairs do not overlap and addresses in detail the question of the degree of enhancement of TPA by quantum correlations of photon number and in spectral properties. This Tutorial includes a review of known theoretical methods and results, as well as some extensions, especially the comparison of TPA processes that occur via far-off-resonant intermediate states only and those that involve off-resonant intermediate states by virtue of dephasing processes. A brief discussion of the main challenges facing experimental studies of ETPA is also given.

We consider the three well-known “quantum pathways” in the fourth-order perturbation theory of the density matrix and clarify to what extent the TPA probability can be modeled by the conventional second-order perturbation theory using state amplitudes. The treatment contains results not previously published, to our knowledge, regarding the roles of the coherent excitation of multiple nonresonant intermediate states. These pathways are known as “step-wise” pathways or “non-rephasing” and “rephasing” pathways, whereas the conventional second-order perturbation theory includes only the direct “two-quantum” or “double quantum coherence” pathway that proceeds through a set of nonresonant virtual states. The two-quantum pathway is known to be enhanced by frequency anticorrelations of the exciting photons. A major question that we address is whether the step-wise pathways are similarly enhanced. We find that they are not.

The basis of the current thinking on entangled two-photon absorption goes back to the seminal theoretical papers by Klyshko,^{10} Gea-Banacloche,^{11} and Javanainen and Gould.^{12} These initial studies showed that in the regime of “isolated” EPP—defined as the regime of extremely low flux where not more than two photons on average impinge on the molecule within the field’s coherence time and not more than two photons impinge on the molecule within its electronic dephasing time—the TPA probability for exciting a molecule scales linearly with the photon flux because the photons arrive in pairs. Moreover, for EPP created by spontaneous parametric down conversion (SPDC) pumped by a narrow-band laser, an increased bandwidth of the EPP field does not decrease the absorption probability. This effect occurs because the frequencies of the two EPP photons are not random but rather anticorrelated so that they sum to the fixed frequency of the pump laser. We refer to this effect as enhancement due to spectral correlation. The first of these predictions was verified in an experiment using narrow-band EPP to excite a two-photon transition in a vapor of atomic cesium.^{13} The second prediction found experimental support in atomic systems for a slightly different scenario: TPA in the high flux (squeezing) regime.^{14} Subsequent studies using molecular solutions reported many orders-of-magnitude enhancement of TPA using EPP,^{15,16} but those studies have been called into question by recent experimental studies.^{17–19}

Another issue is whether or not a *single-photon* state excites a molecule any differently than does a weak coherent state. Given that a one-photon transition of a single molecule can absorb only one photon at a time, does the molecule “know” what state of light that photon is provided by? The short answer is “no.” However, when more than one molecule is present, even in the absence of direct Coulomb coupling between molecules, quantum entanglement between molecules can occur as a result of them being coupled to the same optical field.

We first introduce the overall problem of one- and two-photon absorption of quantum light. We then review the formalism of the quantized electromagnetic (EM) field, including the concept of temporal modes. We then find and apply perturbative solutions to the density matrix equations of motion to derive excitation probabilities for one-photon and two-photon absorption, emphasizing differences between classical (coherent state) excitation and quantum (isolated EPP) excitation. We review and extend a recent derivation predicting the amount of quantum enhancement of TPA that can be achieved using isolated EPP compared to the same flux of classical light, from which it is concluded that observation of two-photon excitation of molecules by isolated EPP is challenging in practice.^{20–22}

Dispersion of the EPP field by passage through typical linear-optical elements is known to decrease the efficacy of TPA by EPP. Our formalism allows incorporating rigorously the effects of dispersion in TPA, leading to a simple formula quantifying the expected decrease in the EPP-driven TPA probability.

Several theoretical studies with goals similar to ours have appeared previously but considered only a subset of the issues we treat. Dayan emphasized cases with no resonant intermediate states, which is also our focus, and also treated the case of overlapping EPP at high flux.^{23} Schlawin and Buchleitner emphasized cases in which resonant intermediate states play a strong role in TPA.^{24}

## II. THEORY OF ONE- AND TWO-PHOTON ABSORPTION—GENERAL FORMALISM

The theory of interaction of molecules with quantized light, including one- and two-photon absorption, has been treated numerous times. The traditional approach, first developed by Maria Göppert-Mayer, uses perturbation theory for state amplitudes and posits a density of molecular or photonic states to arrive at the Fermi Golden Rule for the *conventional cross section* for TPA.^{25–27} Accessible textbook treatments are given in the quantum optics text by Loudon^{28} and in the nonlinear optics text by Boyd.^{29} When dealing with short light pulses or light having time–frequency correlations (entanglement), a more detailed treatment is needed, and several such treatments have appeared, a few being Refs. 1, 11, 12, 20, and 21.

Figures 1 and 2 define the molecular states of interest and the perturbative sequences for the molecular density-matrix elements leading to excitation by one- or two-photon transitions. Molecular density-matrix elements are given by the total density operator $\rho \u0302$ as

where $j$ is a molecular energy eigenstate (“ket”) and $j$ is its corresponding “bra.” The arrows in the figures correspond to single interactions with the field. Note that a single interaction creates a “coherence” between two states, while a second interaction is needed to create a population. To clarify, if there is no dipole dephasing present, a “population” *ρ*_{ee} is simply the square of a “coherence” *ρ*_{eg}, and the squaring process corresponds to a second interaction as in Figs. 2(e) and 2(f). We will return to these diagrams repeatedly.

We follow, as far as possible, the standard literature in ultrafast spectroscopy, where it is common to use a molecular density-matrix formalism to model absorption or nonlinear-optical response of multilevel systems.^{1–4,30} In a recent paper,^{20} we showed how to obtain the correct results for the Double quantum coherence (DQC) pathway of Fig. 2(d) using an alternative formalism that is common in the quantum optics literature—the “operator optical Bloch equations.”^{31} These equations respect the quantum nature of the field via commutation relations involving raising and lowering operators for the field and the molecular states. However, for treating the Non-rephasing (NRP) and Rephasing (RP) pathways of Figs. 2(e) and 2(f), the operator optical Bloch equations become cumbersome, and the density matrix method appears to be more straightforward.

Consider a single molecule interacting with a quantum light field. The molecular energy eigenstates satisfy the eigenequation $H\u0302Mi=\u210f\omega ii$. The density matrix equations track the time evolution of molecular-state density-matrix elements *ρ*_{ji}, which equate to populations (that is probabilities) for *i* = *j* and “coherences” for *i* ≠ *j*. In the Schrödinger picture, the density matrix equation of motion for the combined molecule-field system (without damping or dephasing interactions) is $\u2202t\rho \u0302=(1/i\u210f)[H\u0302,\rho \u0302]$, where the Hamiltonian is $H\u0302=H\u03020\u2212d\u20d7\u22c5\Sigma \sigma e\sigma E\u0302\sigma $, and the unperturbed Hamiltonian is $H\u03020=H\u0302M+H\u0302F$, with $H\u0302F$ being the energy of the field; $d\u20d7$ is the electric-dipole-vector operator, and $e\sigma E\u0302\sigma $ is the field operator with **e**^{σ} being its linear polarization vector, where the sum over *σ* takes on two values (*σ* = 1, 2).

The electric field operator, at the location of the molecule, is written as $E\u0302\sigma (t)=E\u0302(+)(t)+E\u0302(\u2212)(t)$, where $E\u0302(\u2212)(t)$ is a wide-band (all frequencies) photon creation operator and $E\u0302(+)(t)$ is a wide-band annihilation operator. In the semiclassical approximation, one assumes that the field, with state $\rho \u0302F$, is not quantum correlated with the molecular response, described by the state $\rho \u0302M$, so the combined state is the tensor product $\rho \u0302(t)=\rho \u0302M(t)\u2297\rho \u0302F(t)$. Then, the mean-field approximation corresponds to the replacement $Tr\rho \u0302jiE\u0302\sigma \u2192Tr\rho \u0302MjiTr\rho \u0302FE\u0302\sigma =\rho ijE\u0302\sigma $, which recovers the standard molecular density matrix treatment (standard optical Bloch equations) with the field operator replaced by the mean field. This approximation precludes treatment of quantum states of light, which are characterized by the correlation functions $E\u0302(t1)E\u0302(t2)E\u0302(t3)\u2026$. Therefore, in order to allow for the possibility that the field is quantum correlated with the molecular response, we continue to treat the field as an operator.

It is easiest to express the solutions for the density operator in the interaction picture (indicated by *I* subscripts), in which observables are expressed as $O\u0302I(t)=U\u03020\u2020(t,t0)O\u0302U\u03020(t,t0)$, where $U\u03020(t,t0)=exp\u2212iH\u03020(t\u2212t0)$, with *t*_{0} being an arbitrary time at which the interaction picture and Schrödinger picture are equivalent. We set *t*_{0} = 0 for convenience. The density operator in the new picture is

and system operators in this picture are

where we suppressed (*t*, *t*_{0}). The expectation value of an observable is expressed as (using $U\u03020U\u03020\u2020=1\u0302$)

where we used cyclic permutation under the trace and suppressed the time arguments of $U\u03020(t,t0)$ for clarity.

The interaction picture density operator satisfies

in which the interaction Hamiltonian is

The perturbative solution for the density operator is given by the series^{30,32}

where $\rho \u03020=\rho \u0302M(0)\u2297\rho \u0302F(0)$ is the initial state of the system, assumed to be uncorrelated long before the interaction begins (*t* = −*∞*), and the *n*th iterate is

It is worth noting that we have not set *t*_{0} to be the same as the time at which the initial state is specified, thereby letting us simplify the equations.

The ket–bra operators $ij$ (like Pauli matrices) can be used to represent any molecular operator, such as the electric dipole operator, $d\u0302I(t)=d\u0302I(\u2212)(t)+d\u0302I(+)(t)$, where (setting *t*_{0} = 0)

and the electric-dipole matrix elements are $dij=id\u20d7\u22c5e\sigma j$. These operators act as molecular raising and lowering operators. Note that interaction-picture operators referring to different degrees of freedom (such as molecule and field) commute, even at different times, because they equal Schrödinger-picture operators multiplied by complex-valued (non-operator) functions of time.

In order to identify clearly the various quantum pathways leading to the molecular populations and coherences of interest, we will analyze the perturbative solutions of the equations of motion. Before we do that, we review the quantum theory of light.

## III. QUANTIZATION OF TRAVELING-WAVE OPTICAL FIELDS

How does light travel through vacuum? At every point, even in vacuum, there exist an infinite number of electromagnetic harmonic oscillators, one for each of a continuum of frequencies, which are coupled via the Maxwell equations to their same-frequency neighbors, creating propagation. The optical “modes” of the system are the collective degrees of freedom of all the oscillators that form natural solutions of the Maxwell equations.

The electric field **E**(**r**) represents the summed amplitude of all the oscillators at point **r**. In this view, a “photon” is not a tangible object. The word photon is used merely to name various states of the field: a one-photon state of the field, a two-photon state of the field, etc. This viewpoint, developed by Dirac, is the most common one used in quantum optics, although you will often hear the word photon bandied about in not-so-cautious ways.

In the quantum theory of a harmonic oscillator (e.g., a mass on a spring), the raising operator $a\u0302\u2020$ increases the energy by one unit, given by *ℏω*_{0}, where *ω*_{0} is the angular frequency of the oscillator. In the quantum theory of light, a given raising operator $a\u0302j\u2020$ increases the excitation number in the corresponding mode by one: *n* → *n* + 1. If a raising operator acts on the lowest ground state of the system (the “vacuum”), a one-photon state of the field is created. A raising operator $a\u0302j\u2020$ is also called a “creation operator” not because it creates an object called a photon but because when acting on the vacuum, it creates an excitation of the field, described by a state given the name “photon.” That is, a photon is best thought of as a state of the field, not as a particle.^{33,34}

Many molecular spectroscopic studies are carried out in a liquid or solid host medium. Strictly speaking, excitations of the electromagnetic (EM) field inside a medium should be called polaritons, not photons, because they are collective excitations of the field and electronic states of the medium together.^{35,36} For ease of language, we will continue to use the word photon.

Note that a one-photon state of the field need not be monochromatic because it may be spread over many spectral modes in a coherent quantum superposition.^{37} For example, if the state is generated by spontaneous emission by an excited-state atom, the emitted light is a wave packet and its spectrum is that of the natural linewidth associated with the lifetime of the excited state. Thus, we will discuss in this Tutorial how to quantize the EM field in terms of nonmonochromatic modes, also called temporal modes or pulse modes.^{38,39}

In quantum electromagnetism, the infinite set of raising and lowering operators together are used to construct a global operator $E\u0302(r)$. In the Heisenberg picture, where operators evolve in time and states do not, this operator can be expressed in terms of monochromatic modes $u\u0303j(r)$, each having a particular frequency *ω*_{j}. Here, we focus on cases where propagation occurs approximately in one dimension only (*z*), such as light in a well-collimated beam or in an optical wave guide such as a fiber. We can choose a set of monochromatic beam-like mode functions as **e**_{j}w_{j}(**x**, *z*, *β*_{j})exp(*iβ*_{j}*z*), where **e**_{j} is one of the two polarization unit vectors perpendicular to the direction of propagation. *w*_{j}(**x**, *z*, *β*_{j}) are transverse modes that depend strongly on **x** = (*x*, *y*) and weakly or not at all on *z*. They are orthogonal and square normalized in the transverse variable **x**. Here, *β*_{j} is the longitudinal propagation constant (wave number); periodic boundary conditions require it to have equally spaced values *β*_{j} = *j*2*π*/*L*, where *j* is any integer, typically of order 10^{6} for optical fields. For each value of the longitudinal propagation constant *β*_{j}, there exists a discrete set of transverse modes, which generally have different frequencies depending on their transverse spatial shape. If the medium is only weakly dispersive, the transverse modes with equal *β* are orthonormal in two-dimensions (for any *z* value) according to

When propagation is quasi-one-dimensional, it is useful to adopt angular frequency *ω* as the continuous variable used to label modes, since we typically measure frequency. Note that for each transverse mode *j*, there is a dispersion relation *ω*_{j}(*β*), determined by the solutions of the Maxwell Equations (for example, the Hermite–Gaussian modes propagation in a uniform medium). Considering the optical spectrum in free space as a continuum, the field operator describing a not-too-broad spectral band of light is $E\u0302=E\u0302(+)+E\u0302(\u2212)$, with^{36,40}

The frequency *ω* has units rad/s, *n* is the medium’s refractive index at the central frequency of the spectral band of the field, *ɛ*_{0} is the vacuum permittivity, and *c* is the vacuum speed of light. The dispersion relation is now written as *β*_{j}(*ω*). The continuum operators satisfy the commutator

The so-called negative-frequency part $E\u0302(\u2212)$ is related to the positive-frequency part $E\u0302(+)$ by an operator conjugate: $E\u0302(\u2212)=E\u0302(+)\u2020$. If fields of very different central frequencies are considered, then Eq. (11) needs to be written for each spectral band separately.

Now, consider the operator for the field of a single spatial-polarization mode at the location $xM,zM$ of a molecule. Dropping the mode index *j* and replacing the frequency by its central (mean) value $\omega \u0304$, Eq. (11) becomes (we suppress the *I* subscript indicating the interaction picture)

where the scalar part of the electric field operator is

and $E\u0302(t)=E\u0302(+)(t)+E\u0302(\u2212)(t)$, with $E\u0302(\u2212)=E\u0302(+)\u2020$ and $L(\omega )=(\u210f\omega /2\epsilon 0nc)1/2w(xM,zM,\omega )$, where *w*(**x**_{M}, *z*_{M}, *ω*) is the mode amplitude at the molecule’s location. The phase factor is given by $\theta =\beta (\omega \u0304)zM$, and for simplicity, we can set this to an arbitrary value if only one molecule is being considered. Here and in the following, we use the shorthand notation *đ**ω* = d*ω*/2*π*.

We typically consider a state of the field that has a bandwidth much narrower than its central frequency *ω*_{0}. Then, we approximate the field operator as

where $L0=(\u210f\omega \u0304/2\epsilon 0ncA0)1/2$ and *A*_{0} is the *effective* beam area at the molecule’s location given by $A0\u22121/2\u2261w(xM,zM,\omega \u0304)$. In this approximation, the adjoint operator $E\u0302(\u2212)(t)$ creates a photon at time *t* at exactly the location (**x**_{M}, *z*_{M}).

Strictly speaking, the formalism just given is valid only if the medium is nearly transparent in the spectral range of interest, so only dispersion affects the light.^{35,36} In the present case, we consider only one molecule as an absorber and consider that the host medium is transparent to the light, so the above formalism holds.

## IV. CLASSICAL AND QUANTUM STATES OF LIGHT

Here, we summarize the basic properties of coherent states, single-photon states, and two-photon states. For convenience, we consider the incident light to be in pulses of finite duration. In the case of continuous-wave (CW) excitation, we imagine the field to be made of a series of rectangular pulses with constant mean power and duration *T*, as in Fig. 3(a). When comparing to TPA with short pulses of SPDC or coherent-state light, it suffices to consider only a single pulse occupying the same interaction time window *T*, as in Figs. 3(b) and 3(c).

### A. Coherent state

A pulsed coherent state $\alpha $ with spectral amplitude *α*(*ω*) satisfies $a\u0302(\omega )\alpha =\alpha (\omega )\alpha $ at each frequency^{40} or equivalently in the time domain,

where the field amplitude is

Here, *A*(*t*) is a slowly varying envelope, and *ω*_{0} is the central (carrier) frequency. Because the coherent state is the quantum state that most resembles a field in classical EM theory, we often call the coherent state a “classical” state.

The mean number of photons in the pulse is the time-integral of the flux $A(t)2$,

It is useful to define a unity-normalized spectral amplitude *ϕ*(*ω*) by *α*(*ω*) = *α*_{0}*ϕ*(*ω*), where $\u222b\varphi (\omega )2\u0111\omega =1$. Then, *N* = |*α*_{0}|^{2}.

In the case of a coherent state in a constant-amplitude (“rectangular”) pulse with duration *T*, the needed Fourier-transform pair is

as illustrated in Fig. 4.

### B. Single-photon state

A single-photon state of a particular temporal mode (coherent wave packet) *φ*(*ω*) is described by a superposition of monochromatic one-photon states,

where the spectral density (probability) is normalized as $\u222b\u0111\omega \phi (\omega )2=1$. This state corresponds to light arriving at the molecule in the form of a time-domain wave packet,

Such a state can be created by, for example, heralded SPDC in which one of the photons is detected, announcing or “heralding” the presence of the other.^{41,42}

It is sometimes useful to quantize the field in terms of wave packets rather than monochromatic waves.^{37–40} This can be done by choosing any complete orthonormal set of spectral amplitude functions {*φ*_{j}(*ω*)}, with $\u2211j\phi j*(\omega \u2032)\phi j(\omega )=2\pi \delta (\omega \u2032\u2212\omega )$, and using them to define a set of “wave packet operators,”

These operators satisfy $[a\u0302i,a\u0302j\u2020]=\delta ij$. The continuous set of operators has been converted into a discrete set. One can create a single-photon state of the field by acting such an operator on the vacuum state: $a\u0302j\u2020vac$. This action corresponds to creating a photon in a particular “temporal mode,” as can be seen by inverting Eq. (22),

Hence, the electric field operator Eq. (15) can be expressed as

where v_{j}(*t*) is called a temporal mode,

A useful example of a square-normalized temporal mode spectrum is given by

where Γ is the spectral half-width. Then, the temporal mode is a single-sided exponential

as illustrated in Fig. 5.

### C. Entangled photon pair state

Collinear type-0 or type-I spontaneous parametric down conversion is characterized by the two generated photons having the same polarization. When pumped by a laser pulse of finite duration, the SPDC can be designed to occur into a single spatial-and-polarization mode;^{43} then, the state is described by

where *ɛ*^{2} ≪ 1 is the probability that a given pulse contains a photon pair. Consistent with the labeling of pulses shown in Fig. 3, the mean photon flux (twice the pairs flux), time-averaged over the interval *T*, equals *F*_{EPP} = 2*ɛ*^{2}/*T* for either CW or pulsed cases. We neglect higher-order terms representing generation of multiple pairs in order to satisfy our assumption of an isolated EPP interacting with the molecule. The joint-spectral amplitude (JSA) $\psi (\omega ,\omega \u0303)$ is determined by energy and momentum conservation, given the spectrum of the pumping laser pulse and the phase-matching properties of the nonlinear crystal used as second-order nonlinear medium.^{44} It is normalized as $\u222b\u0111\omega \u222b\u0111\omega \u0303|\psi (\omega ,\omega \u0303)|2=1$. An important property of the JSA when treating a single mode of the field, for which there are no distinguishing labels on the photon creation operators other than frequency, is the required symmetry $\psi (\omega ,\omega \u0303)=\psi (\omega \u0303,\omega )$. This fact is seen easily by swapping the integration variables in Eq. (28). Such a symmetry is not upheld when treating type-II SPDC, which takes place into two or more distinct modes.

For a long narrowband pump pulse with central frequency *ω*_{P}, the JSA is largest typically along the antidiagonal, that is, $\omega +\omega \u0303=\omega P$, where the frequencies are anti-correlated, as illustrated in Fig. 6(a). The two-photon wave packet in the time domain is given by the double Fourier transform

The photon arrival times are positively correlated to within the so-called *entanglement time**T*_{e}, as shown in Fig. 6(b). The entanglement time is given roughly by the inverse of the bandwidth.

Because Eq. (28) represents the state of a single spatial mode, it does not describe transverse spatial entanglement (correlation). The photons, if detected, are assumed to be distributed across the beam independently, as is the case if their generation takes place in a single-mode wave guide or the beam is spatially filtered by passing through a single-mode optical fiber. Therefore, in this case with EPP generated by SPDC, the “entanglement area” is the transverse area of the EPP beam at the molecule’s location, unlike in treatments where transverse spatial correlations are considered.^{5} For a wide-area beam, such correlations can localize photon pairs to transverse areas much smaller than the overall beam area but cannot localize pairs to an area much smaller than the diffraction-limited focus of a well-designed optical system.

## V. INDUCED DIPOLE, ONE-PHOTON ABSORPTION, AND DEPHASING

### A. First-order solution

Here, we give the details of the lowest-order solution of the density matrix and its application to the creation of optical coherence and one-photon absorption by a single-photon wave packet. For simplicity, again we assume that the field has energy in only one polarization state and drop the *σ* index.

The first-order solution for the density matrix describes the coherent linear response of a molecule,

A useful step in formulating the solutions is to transform to a difference-time variable, *τ* = *t* − *t*_{1}. Then, the solution is

where *hc* stands for the Hermitian conjugate, consistent with the requirement that the density operator be Hermitian.

The macroscopic electric-dipole polarization density is equal to the number density of molecules times the expectation of the dipole operator, given by $d\u0302I(t)=d\u0302I(\u2212)(t)+c.c.$, where, using $\rho \u03020=\rho \u0302M\rho \u0302F$ and permutation inside the trace,

where

Throughout this Tutorial, we assume that the molecule starts in its ground state, so $\rho \u0302M=gg$. Then, as a consequence of $gd\u0302I(\u2212)(t),=0$,

We evaluate the remaining term,

where we introduced the notation for the difference frequencies,

The form of the field expectation value $E\u0302I(t\u2212\tau )$ depends, of course, on the state of the field. For a coherent state with a known phase, it is given by Eq. (15) as

The mean value of the field is zero because there is an unbalanced number of raising and lowering operators in the expectation value, resulting in projecting the vacuum state onto a single-photon state. This result can also be interpreted as saying that a single-photon state has no definite phase. Thus, in this case, $d\u0302I(t)=0$.

However, even though the mean dipole is zero, the single-photon wave packet creates quantum state entanglement between the molecule and field, as can be seen by the following. The density operator to first order is $\rho \u0302I(t)=\rho \u03020+\rho \u0302I(1)(t)$, where the initial state is $\rho \u03020=gg\u2297\phi \phi $, with $\phi $ being the initial field state. Then, from Eq. (31), after inserting the field operator and using the rotating-wave approximation (RWA),

where we assumed near resonance to a particular *g-j* transition. (To ensure convergence, we inserted a small damping constant *ζ* in the exponential and set it to zero at the end). We carry out the integral and write the result as $\rho \u0302I1t=\eta jvacg\phi +hc$, where

We can rewrite the result approximately as

which we see represents a pure state $g\phi +\eta jvac$.

The resulting state is entangled, that is, it cannot be written in product form $\psi M\u2297\phi F$. Such entanglement provides the basic resource for many quantum-information techniques.

Note also that if two molecules are exposed to the same single-photon pulse and the photon is later observed to have been absorbed (it fails to trigger a “perfect” detector), then the two molecules are left in an entangled state: they share the excitation of one photon.

Such entanglement is the basis for atomic-ensemble quantum memories, which store states of light coherently in an extended “phased-array” of atoms.^{45} Such an entangled state of many atoms can be thought of an exciton even if the atoms are distant from one another. If the electron spin is part of the state labels, this state is often called a spin wave.

As a second example, if the molecule is driven by a *coherent state*, the polarization density is

To simplify writing the result, we will make the rotating-wave approximation (RWA), in which only the near-resonant term *A**(*t* − *τ*) is kept. The non-resonant (counter rotating) terms do contribute as a correction, which is often, but not always, small. We can express the result as

where we transformed back to *t*_{1} = *t* − *τ*. We observe several facts from this result. If we evaluate the result at very long times, after the pulse has come and gone, then each term in the sum is proportional to the Fourier transform of *A*(*t*), that is, the spectral amplitude of the pulse, evaluated at the corresponding molecular transition frequency. In this case of coherent “impulsive excitation” (and ignoring damping and dephasing interactions), the polarization density oscillates at all of the excited molecular frequencies *ω*_{kg}.

We leave as an exercise to show that after a coherent-state pulse interacts with the molecule, there is no entanglement between the molecule and the field. That is because when a pure coherent state experiences loss or absorption, it is transformed into a pure coherent state of lesser amplitude. By definition, a pure state is not entangled with any other system.

### B. Kubo theory of molecular dephasing

Now, we consider the effects of dephasing interactions of the molecule with its surrounding environment, which may be a gas, liquid, or solid. To treat environmental perturbations rigorously, one should introduce new terms into the Hamiltonian in Eq. (6). It turns out that a much simpler approach yields a realistic model—the Kubo theory of line shapes, in which the environment is treated as a (semiclassical) random process that causes the molecular energies to fluctuate randomly leading to finite transition linewidths. For review, see Refs. 30, 32, and 46–48. To avoid a lengthy treatment, we illustrate only the simplest case of Kubo theory, that with zero-duration impact-like interactions of the molecule and the environment, leading to Lorentzian homogeneous line shapes.

In this phenomenological treatment, the molecular energies *ω*_{jk} are replaced by *ω*_{jk} + *x*_{jk}(*t*), where *x*_{jk}(*t*) are random fluctuating quantities (frequencies) with zero mean. The factors of the form exp(−*iω t*) (dropping the *jk* subscript) that appear in equations such as Eq. (42) are replaced by factors of the form

where the random accumulated phase is

Here, the brackets indicate an ensemble average (not a time average) over possible realizations of the random process *x*(*t*). It is assumed that the random process *x*(*t*) [and thus also *ϕ*(*t*)] obeys Gaussian statistics, namely, the probability for *ϕ*(*t*) to take on a particular value is given by

where the variance of *ϕ* is given by

Thus, the average in Eq. (44) becomes

In the case of zero-duration impact-like interactions of the molecule and the environment, the correlation function is delta-like, that is, $x(t)x(t\u2032)=2\gamma \delta (t\u2212t\u2032)$, where 2*γ* is the linewidth of the considered transition, as we will see. The integral is easily carried out,

thereby giving

Recall that *t* is positive. We then have, replacing Eq. (42),

The simple version of Kubo dephasing theory does not include the possibility of population damping—decay of the probability for the molecule to be in a particular state resulting from spontaneous emission or other incoherent nonradiative processes that couple two states. While such processes can be treated rigorously (e.g., Refs. 46 and 47), here we treat them phenomenologically by adding additional damping rates into the Kubo dephasing rates. The rules are if the population decay rate out of each state *i* is denoted *A*_{i}, then the population *ρ*_{ii} is damped as exp[−*A*_{i}*t*] and the off-diagonal density matrix element *ρ*_{ij} of the transition *i* ↔ *j* is damped as exp[−(*A*_{i}/2 + *A*_{j}/2)*t*] (e.g., Ref. 46, Sec. 14.5). Population damping causes a broadening of the spectral lines, called lifetime broadening. There can also be population-increasing processes, such as spontaneous emission into a state, but for our purposes, these processes will not be important because we treat the problem perturbatively and typically on short time scales compared to the times scales of such effects.

### C. Steady-state induced dipole

Now that we have dephasing included, we can evaluate the polarization density in steady state, in which the field amplitude *A* is constant,

We see that in the steady-state regime, the polarization oscillates at the driving frequency *ω*_{0}, whereas we saw earlier that in the impulsive regime, it oscillates at the molecular frequencies. These results are, of course, the same as obtained in a semiclassical treatment where the field is treated as a classical one.

### D. One-photon-induced population

#### 1. General case

The population $Pe(2)$ in a given excited state $e$ resulting from one-photon absorption equals the expectation value $Pe(2)=Tr\rho \u0302I(2)(t)ee$, where $\rho \u0302I(2)(t)$ is the density operator solution iterated to second order,

Because the initial state is $\rho \u0302M=gg$, we note that $\rho \u03020ee=ee\rho \u03020=0$, and thus, only two of the four terms contribute to the trace,

For near-resonant excitation, we again apply the RWA, retaining only terms in which the molecule and field exchange excitations,

The neglected terms will contribute small corrections to the final result, which are significant only when the excitation is far from resonance. Using Eq. (9) for the dipole operators, we find for the first term in Eq. (54),

where *ω*_{ij} = *ω*_{i} − *ω*_{j}. The second term in Eq. (54) is the same with *t*_{1} and *t*_{2} swapped. Hence, using permutation inside the trace and $Tr(A\u0302\u2020)=Tr(A\u0302)*$ and the initial state $\rho \u03020=gg\u2297\rho \u0302F$, we have for the population, or probability to find the molecule in state *e*,

where the second-order field correlation function is

In the final line, we took the long-time limit, *t* = *∞*, to find the population immediately after the pulse has passed through the molecule. The final population is equivalent to the pathway represented in Figs. 1(b) and 1(c). The second-order field correlation function, when divided by $E\u0302I(\u2212)(t1)E\u0302I(+)(t1)1/2E\u0302I(\u2212)(t2)E\u0302I(+)(t2)1/2$, is often referred to as the degree of first-order temporal coherence *g*^{(1)}(*t*_{1}, *t*_{2}).^{28}

Now, we insert the dephasing factor according to Kubo theory, transform to a difference-time variable, *τ* = *t*_{2} − *t*_{1}, and rename the remaining integration variable *t*_{2} ≡ *t* Then, the solution is

It is useful to transform to the frequency domain by inserting Eq. (15) twice,

This result has a simple interpretation—the photons at frequency *ω*, whose number is $a\u0302\u2020(\omega )a\u0302(\omega )$, drive the transition with the strength given by the Lorentzian line shape.

#### 2. One-photon absorption: Coherent state

With these results in hand, we can evaluate the one-photon excitation probability for a coherent-state pulse, for which $E\u0302I(\u2212)(t2)E\u0302I(+)(t1)=L02A*(t2)A(t1)ei\omega 0(t2\u2212t1)$, and we find

In the frequency domain, for which $a\u0302\u2020(\omega )a\u0302(\omega )=\alpha *(\omega )\alpha (\omega )$, we have $a\u0302\u2020(\omega )a\u0302(\omega )=\alpha 02\phi *(\omega )\phi (\omega )$ and

For quasi-monochromatic light, we can write the result in terms of the absorption cross section $\sigma 1$. The probability to excite the molecule per photon incident should equal the ratio of *σ*^{(1)} to the optical beam’s cross-sectional area *A*_{0}. Recall that the mean number of photons is $N=\u222b\alpha (\omega )2\u0111\omega $. Then, if the spectrum of light is centered at *ω*_{0} and narrow compared to the Lorentzian absorption line, we have (using the definition of *L*_{0})

where the one-photon cross section (with units *m*^{-2}) is

This result agrees with, for example, Ref. 49.

#### 3. One-photon absorption: Single-photon state

Next, we evaluate the one-photon excitation probability when driven by a single-photon state, Eq. (20). First, we note that operating the annihilation operator on the state gives

Thus, we have $a\u0302\u2020(\omega )a\u0302(\omega )=\phi *(\omega )\phi (\omega )$, and Eq. (59) gives

We see that the probabilities to excite the state *e* by a coherent state or by a single-photon state have precisely the same form, except that for the former we have *∫**đ**ω*|*α*(*ω*)|^{2} = *N* and in the latter we have *∫**đ**ω*|*φ*(*ω*)|^{2} = 1. Therefore, when exciting a single molecule in the linear-response regime, a single-photon pulse has the same effect as a coherent-state pulse with the mean photon number equal to one. In this scenario, there is nothing especially “quantum” about single-photon absorption. As pointed out earlier, the story is different when exciting two or more molecules—the single-photon pulse can create entanglement of excitation in the molecules’ joint state, whereas a coherent-state pulse does not.

## VI. TWO-PHOTON ABSORPTION: GENERAL TREATMENT

Here, we give the results of the fourth-order solution of the density matrix, setting up its application to TPA by coherent states of EPP, and we discuss the “quantum advantages” of EPP relative to coherent states.

### A. Conventional TPA cross section

Before giving general results, we touch base with the well-known lore on the two-photon cross section, first calculated in 1931 by Maria Göppert-Mayer, named a Nobel laureate in Physics in 1963. Of special interest is the role, if any, of the NRP and RP pathways in far-off-resonance TPA. In the conventional theory, these terms do not appear, and the question is why? The conventional theory uses second-order perturbation theory for state amplitudes, rather than density-matrix elements, assuming the intermediate states are far from resonance with the exciting field’s frequency, and so dephasing rates are ignored (as they must be in that treatment). See Ref. 29 (Sec. 12.5.3) for a concise review. The result is expressed in terms of a TPA cross section,

where 1/2*πγ*_{fg} plays the role of the density of states for the transition. We will see how to reproduce this result using the density matrix approach, which provides further insight into the role of the different pathways.

### B. Fourth-order TPA solution

In the density matrix approach, the fourth-order solution, needed to calculate the probability for TPA, is

where the sum is over *p*, *q*, *r*, *s* = ±, and for simplicity, we made the RWA and used the compact notation

where $\mu \u0302(\u2212)(t)=d\u0302(\u2212)(t)/\u210f$ is a scaled dipole operator. $V\u0302(+)$ raises the molecule and lowers the field. $V\u0302(\u2212)$ does the opposite. The solution has 16 terms contributing to the sum, many of which can be neglected in most cases. We focus on those TPA terms that lead to population in the *f* state, $Pf=Tr\rho \u0302(4)ff$. Here and in the following we drop the subscript *I* indicating the interaction picture.

Given that the transition of interest is *g* → *f*, the only terms resulting from the nested commutators that are nonzero are those of the form $V(+)V(+)\rho \u03020V(\u2212)V(\u2212)$, giving

where we defined

To interpret these terms, consider that all operators act toward $\rho \u03020$ in the order *t*_{1}, *t*_{2}, *t*_{3}, *t*_{4}. Referring to Fig. 2, we can see clearly the correspondence of each term with each diagram.

Because we are using the interaction picture, operators for the field commute with those of the molecule. Thus, we can separate the correlation functions as follows, keeping in mind that $\rho \u0302M=gg$:

where the molecule correlation functions are

These are the same forms that appear in the semiclassical treatment, where the field is treated as a classical function.^{30,32,46} They are evaluated in Appendix A, which derives the following results using Kubo dephasing theory:

where we introduced the difference-time variables *r* = *t*_{4} − *t*_{3}, *s* = *t*_{3} − *t*_{2}, *τ* = *t*_{2} − *t*_{1}. These variables track the time increases during disjoint time intervals during which dephasing takes place, allowing the dephasing during each interval to be considered separately as in Eq. (74). For each transition at frequency *ω*_{ij}, the corresponding dephasing rate is *γ*_{ij}. Note the plus sign in the exponential argument (*γ*_{fe} + *iω*_{fe})*r* of the RP correlation function, which is related to the well-known effect of photon echoes in special cases.^{31,30}

The field correlation functions in Eq. (72), which can describe quantum states of light, are given by

We can summarize them using a common form, where *t*_{a}, etc., are arbitrary (“dummy”) variables,

where we used cyclic permutation inside the trace. This fourth-order field correlation function, when divided by $E\u0302I(\u2212)(t1)E\u0302I(+)(t1)E\u0302I(\u2212)(t2)E\u0302I(+)(t2)$, is often referred to as the degree of second-order temporal coherence *g*^{(2)}(*t*_{c}, *t*_{d}, *t*_{a}, *t*_{b}).^{28}

For pure states, it can be written as

In the case of a pure two-photon state $\Psi (2)$, the correlation function can be written as

where we inserted unity, $\u2211\Psi \Psi \Psi $, in the center of Eq. (77) and noted that only the vacuum term contributes. It is notable that in this case, the correlation function equals the product of two functions, known as the two-photon detection amplitude,

This result suggests that we can view TPA in a molecule as a two-photon detector.^{9}

### C. Two-photon amplitude

The form of the two-photon detection amplitude reveals interesting aspects of quantum optics, namely, the role played by the boson nature of photons when viewed (with due circumspect) as particles. We insert the two-photon component of the state in Eq. (28) into the expression [Eq. (78)] and insert the frequency representation of the field operators from Eq. (15),

where the frequency-domain correlation function in the vacuum state is

Using the commutator $[a\u0302(\omega ),a\u0302\u2020(\omega \u2032)]=2\pi \delta (\omega \u2212\omega \u2032)$, we find

Then, Eq. (79) becomes

where in the last line, we swapped the $\omega ,\omega \u0303$ variables. We thus see the interesting result that the two-photon detection amplitude depends only on the symmetrized form of the JSA, which we denote as

However, recall we are treating, for simplicity, only the case that the driving light is in a single spatial-polarization mode of the field. After Eq. (28), we noted that in this case, the JSA must be symmetric, that is, $\psi (\omega \u0303,\omega )=\psi (\omega ,\omega \u0303)$. Thus, in this case $\Psi (\omega ,\omega \u0303)=\psi (\omega ,\omega \u0303)$. We can thus write

The mod-square |Φ(*t*_{a}, *t*_{b})|^{2} is the joint probability to detect a photon at time *t*_{a} and a photon at time *t*_{b}, presuming one has sufficiently fast detectors. If one separates the field into its spectrum, say, using a prism, then $|\Psi (\omega ,\omega \u0303)|2$ is seen to be the joint probability to detect a pair of photons at the two indicated frequencies.

We retain the notation $\Psi (\omega ,\omega \u0303)$ because it can be shown straightforwardly to arise automatically even in the case of Type-II SPDC, where the signal and idler modes as distinct, so the state is written with labels on the creation operators,

The JSA here need not be symmetric, but the symmetrized form still determines the two-photon detection amplitude as in Eq. (83). This fact illustrates an important point in quantum optics: modes of the field are distinct and therefore in the quantum sense, they are distinguishable. Therefore, the joint state, $\psi (\omega ,\omega \u0303)$, of two modes need not be symmetric under label exchange. On the other hand, when viewing light as made of photons, which in the quantum sense are indistinguishable “particles,” their joint state $\Psi (\omega ,\omega \u0303)$ must be symmetric. The needed symmetry is automatically satisfied by the math of boson commutators.

### D. TPA probabilities

Transforming to difference-time variables, *r* = *t*_{4} − *t*_{3}, *s* = *t*_{3} − *t*_{2}, *τ* = *t*_{2} − *t*_{1}, and the inverse, *t*_{1} = (*t*_{4} − *r* − *s* − *τ*), *t*_{2} = (*t*_{4} − *r* − *s*), *t*_{3} = (*t*_{4} − *r*), and taking *t* to infinity to encompass the entire excitation pulse, we find from Eqs. (70) and (74),

where again we denote the product of matrix elements by $Mfe\u2032eg=\mu fe\u2032\mu e\u2032g\mu ge\mu ef$ and

These expressions are equivalent to those in Refs. 30, 32, and 46, generalized here to arbitrary states of the field.

Inserting the frequency-domain form of the field operators, Eq. (15), and after some laborious math, the time integrals can be evaluated to give

### E. Four-frequency correlation functions and quantum advantage

If the exciting pulse is a pure coherent state such that $a\u0302(\omega )\alpha =\alpha 0\varphi (\omega )\alpha $, then the frequency-domain correlation function is

where we used that the mean number of photons in the pulse is *N* = |*α*_{0}|^{2}. The temporal amplitude, from Eq. (17), is

which acts essentially like a “classical” optical pulse envelope.

Turning to “quantum” light, if the exciting pulse is a one-photon pulse, then the correlation function equals zero—there can be no TPA. If the exciting pulse is a two-photon pulse, then the correlation function equals

which is the frequency-domain counterpart of Eq. (78). As we did in Eq. (80), we insert the two-photon component of the state in Eq. (28) into the right-most factor in expression Eq. (92),

where we used the delta-function form of the vacuum correlation function $C\omega ,vac(4)$ from Eq. (82) and the symmetrized state $\Psi (\omega ,\omega \u0303)$ in Eq. (84). Then, the four-frequency correlation function for the two-photon state is

Here, we see two essential differences between the classical-state and quantum-state cases as follows:

- The correlation function for the coherent state is the product of a function of a single variable evaluated at four frequencies, which for the two-photon state it is the product of a function of two variables evaluated at two pairs of frequencies. The two descriptions can be compared by identifyingWhen a function of two variables equals the product of two single-variable functions, we say it is(95)$\Psi (\omega ,\omega \u0303)\u225c\varphi (\omega )\varphi (\omega \u0303).$
*separable*. Thus, in this case, we can identify a coherent state as being identical to a two-photon state in its spectral properties. On the other hand, the quantum state $\Psi (\omega ,\omega \u0303)$ has the possibility to represent correlations of the frequencies of the two photons. These are the entangled photon pairs (EPP), and we will see such a correlation can offer a kind of*quantum advantage*in that it can enhance the probability of TPA. The mean number of photons for the coherent state is

*N*, while for the two-photon state, the mean number of photons is twice the mean number of pairs or 2*ɛ*^{2}. Thus, from Eqs. (90) and (94), we see that the correlation function and thus the*f*-state probability scales*quadratically*with the mean photon number for the coherent state, while for the two-photon state, it scales*linearly*with the mean photon number. This linear scaling may offer a second kind of quantum advantage when using a two-photon state to drive TPA in cases where the photon flux is extremely low.

### F. Far-off-resonance approximation

If the exciting field is far from resonance with any intermediate states, as in Fig. 2(a), we approximate Eq. (89) as

where *ω*_{0} is the central frequency of the driving-field spectrum, and we still abbreviate $\omega \u0303\u2032\u2261\omega +\omega \u0303\u2212\omega \u2032$. We can gain some insight by examining the denominators inside the integrals. For DQC, the denominator is minimized when $\omega +\omega \u0303=\omega fg$, that is, the two photon’s frequencies sum to the two-photon resonance frequency. For NRP and RP, the denominators are minimized when $\omega \u2032\u2212\omega =\omega ee\u2032$, that is, the difference of the two photon’s frequencies equals either zero (for *e*′ = *e*, meaning the pathway goes through a “real” population of an intermediate state) or $\omega e\u2032e\u22600$ (meaning the pathway goes through a “coherence” between two states *e* and *e*′).

In this Tutorial, we focus on the case that the DQC term is two-photon resonant or near-resonant, which means that the NRP and RP terms are far-off resonance in most practical situations. An exception is a molecule such as a dimer consisting of two like monomers in which the two-photon resonance exciting the doubly excited state has frequency very near twice the monomer single-photon absorption transitions. In such cases, one should consider the limitations imposed by the rotating-wave approximation and consider additional terms in the expansion of the nested commutators in Eq. (68).

In all the following, we will assume the condition of far-off-resonance intermediated states. We will show that under this condition, the DQC pathway dominates and that it may be strongly enhanced by time–frequency entanglement.

## VII. TWO-PHOTON EXCITATION BY COHERENT STATES

Here, we discuss two-photon excitation by coherent states in several scenarios. The examples serve as baselines for the comparisons to ETPA in Sec.VIII.

For excitation by a coherent state, we insert Eq. (90) into the far-off-resonance expressions Eq. (89) to give

where again $\omega \u0303\u2032\u2261\omega +\omega \u0303\u2212\omega \u2032$.

### A. Coherent-state DQC pathway

Consider a coherent-state pulse with an arbitrary shape and corresponding spectral amplitude *ϕ*(*ω*). In the first line of Eq. (97), we change variables to $x=\omega +\omega \u0303\u22122\omega 0,z=\omega \u2212\omega 0$ (and *z*′ = *ω*′ − *ω*_{0}) and obtain

where we define

where, for convenience, we denote

*K*_{coh}(*x*) is a projection along anti-diagonal lines, as illustrated in Fig. 3. It represents the different combinations of frequencies that effectively create excitation near $\omega +\omega \u0303=2\omega 0$.

As shown in Fig. 7, we also define a *marginal spectrum* for the coherent state as the *vertical* projection onto the *ω* axis,

which is identical to the standard spectrum as measured by a grating spectrometer with a linear-response detector. In contrast, the anti-diagonally projected spectrum *K*_{coh}(*x*) is that “felt” by the molecule, which acts as a spectrally selective two-photon detector.

### B. Coherent-state exponential pulse: DQC pathway

A nice example that permits an analytical result is TPA excited by a coherent-state pulse in the form of a single-sided exponential, as in Eq. (27),

for which the mean total number of photons is $N=\alpha 02=A(0)2$ and the spectral amplitude is

and

We see that in this model, the coherent-state bandwidth simply adds to the dephasing linewidth of the TPA transition. The population created by the DQC pathway, on two-photon resonance, from Eqs. (87) and (106), is

### C. Rectangular coherent state pulse and TPA cross section

An absorption cross section is an effective area that describes the probability per second that a photon (or photons) will be absorbed from a beam with constant flux *F*_{coh} (photons/s) and given area *A*_{0} (m^{2}). To determine the TPA cross section, consider a “rectangular” coherent-state pulse that suddenly turns on with constant amplitude for a duration *T* and is zero afterward, as in Eq. (19). The square-normalized spectral amplitude is

The anti-diagonal-projected spectrum, from Eq. (99), is then

This leads to the result, if the field is two-photon resonant, *ω*_{fg} = 2*ω*_{0},

The exponential term here is a turn-on transient. If the exciting field has long duration and thus is quasi-monochromatic, that is spectrally narrow compared to the TPA transition linewidth (1/*T* ≪ *γ*_{fg}), then Eq. (110) becomes

To write this probability in terms of a cross section, note that the instantaneous constant photon flux for a “rectangular” pulse is *F*_{coh} = *N*/*T*; therefore, $N2/T=(N2/T2)T=Fcoh2T$. We can thus define a two-photon cross section by writing the probability increase per second in terms of a flux density *F*_{coh}/*A*_{0} (photons per s per m^{2}) as

where

where we used $L0=(\u210f\omega 0/2\epsilon 0ncA0)1/2$ and *μ*_{ij} = *d*_{ij}/*ℏ* = **d**_{ij} · **e**/*ℏ*. This result is precisely the far-off-resonance two-photon cross section derived using second-order time-dependent perturbation theory and averaging over the density of states, as can be seen by comparing the double sum to

and noting that at TPA resonance, $\omega e\u2032g\u2212\omega 0=\u2212\omega fe\u2032+\omega 0$. Note that to achieve agreement between our result and the conventional one, we had to consider that the pulse is long enough (*T* ≫ 1/*γ*_{fg}) to allow neglecting the turn-on transient, as is usual in using perturbation theory for determining rates. However, we cannot take the pulse too long because according to this perturbative theory, the population increases linearly in time; that is, there is no steady-state limit in this treatment.

For molecules in solution, $\sigma 2$ is typically of the order of 1 to 1000 GM (where 1 *GM* = 10^{−58} m^{4}s).^{29} Then, for a steady 1-W laser beam with wavelength 800 nm collimated to an area 5 *μ*m^{2}, the flux is *F*_{coh} = 4.0 × 10^{18} photons/s and the flux density is *F*_{coh}/*A*_{0} = 8 × 10^{29} photons/m^{2}s, and the expected TPA rate per molecule is about 64/s. For a 1-ns pulse, roughly the upper limit for the applicability of the present perturbation theory, the probability to excite a given molecule is thus 6.4 × 10^{−8}. Given that 4.0 × 10^{9} photons pass through the beam area at the molecule’s vicinity, we infer an extremely small TPA efficiency per photon per molecule. This low efficiency is one of the motivations for exploring whether time–frequency entanglement of photons can greatly increase this efficiency.

For completeness, Eq. (112) can be generalized for a long arbitrarily shaped quasi-monochromatic light pulse with a time-dependent flux *F*_{coh}(*t*) and a nonzero detuning between 2*ω*_{0} and *ω*_{fg} to

where we used the Fourier relation

### D. Gaussian coherent-state pulse: DQC pathway

As an another example of TPA, consider a Gaussian-pulse coherent state having a spectrum with bandwidth ∼*σ*, i.e., $\alpha (\omega )=\alpha 0(\sigma 2/2\pi )\u22121/4\u2061exp(\u2212(\omega \u2212\omega 0)2/4\sigma 2)$, and duration ∼1/*σ*. The total number of photons is $N=\u222b\alpha (\omega )2\u0111\omega =\alpha 02$. Then, we find for the anti-diagonally projected spectrum *K*_{coh}(*x*) = exp(−*x*^{2}/8*σ*^{2}), and for two-photon resonance, 2*ω*_{0} = *ω*_{fg}, the excitation probability is, from Eq. (100),

where Σ^{(2)} is defined in Eq. (101) and where *ξ*(*z*) = exp(+*z*^{2})erfc(*z*) [sometimes denoted as erfcx(*z*)], which has a maximum value 1 at *z* = 0 (an ultrashort pulse). For a long quasi-monochromatic pulse (large *z*), *ξ*(*z*) decays to zero as ∼1/(*π*^{1/2}*z*). We summarize these two limits as

We see, again, the expected quadratic dependence on the number of photons. Note that in the impulsive limit *σ* → *∞*, the probability does not depend on *γ*_{fg} or *σ* because for the coherent (“instantaneous-response”) nonlinear TPA process, the effect of spreading the spectrum over a broader range is compensated by the increasing peak intensity in the time domain, as verified in Appendix C also for an arbitrary pulse shape. In the opposite limit *σ* → 0, where the pulse duration tends to infinity, the probability goes to zero because the fixed number of photons are spread over a longer and longer time interval, decreasing the chances for accidental coincidences.

### E. Coherent-state NRP and RP pathways

The two pathways labeled NRP and RP in Fig. 2 may also contribute to the probability $Pf=Tr\rho \u0302(4)ff$ for exciting the *f* state. In these pathways, TPA proceeds through a “step-wise” process through the population of state *e*, *ρ*_{gg} → *ρ*_{eg} → *ρ*_{ee} → *ρ*_{fe} → *ρ*_{ff},^{54} or through a coherent process $\rho gg\u2192\rho eg\u2192\rho ee\u2032\u2192\rho fe\u2032\u2192\rho ff$, in which two intermediate states *e* and *e*′ are excited coherently.

We focus on the case that twice the center frequency of the exciting light is near two-photon resonance with the *f*–*g* transition, that is, 2*ω*_{0} ≃ *ω*_{fg}, where the rotating-wave approximation is most reliable. We address whether or not the NRP and RP pathways contribute significantly to the *f*-state probability in comparison to the DQC pathway in the case of coherent-state excitation.

Recall that the DQC pathway creates a TPA cross section in exact agreement with the Göppert–Mayer perturbation theory after any initial transients have damped out, as shown in Sec. VI. We would like to learn whether the NRP and RP terms together can nevertheless contribute to the cross section under certain conditions.

To evaluate the NRP and RP terms, we write the second and third lines of Eq. (96) as

where for arbitrary states of light,

First, we observe from Eq. (119) that for the terms with *e*′ = *e* so that $\omega fe\u2032=\omega fe$, we have $Re,eNRP+Re,eRP=0$ in this far-off-resonance situation where dephasing is unimportant. Thus, if there is only a single relevant intermediate state, then the NRP and RP pathways cancel and do not contribute. If the dephasing linewidth *γ*_{fe} becomes significant, as in the original Eq. (89), it can “destroy” the near-perfect cancellation of $Re,e\u2032NRP$ and $Re,e\u2032RP$, leading to active pathways for TPA that do not exist in the absence of fast dephasing. Study of these pathways requires a more detailed analysis than presented in this tutorial. This “destruction of destructive interference” is similar to that studied in different contexts in, for example, Refs. 50–52, and reviewed in Refs. 53 and 54.

Hence, in the general case, only pathways with *e*′ ≠ *e*, in which two intermediate states *e* and *e*′ are excited coherently, contribute to the *f*-state probability. To analyze cases with *e*′ ≠ *e*, in Eq. (120), we change variables to *y* = *ω* − *ω*′ and write it in the form

where, for a coherent state,

This expression is valid for excitation by coherent states or two-photon states. For a coherent state, it becomes

If we define an effective “two-photon amplitude” for the coherent state as $\Psi coh\omega ,\omega \u0303=\varphi \omega \varphi *\omega \u0303$, then we have

This integral is an autocorrelation of a two-dimensional function $\Psi coh(\omega ,\omega \u0303)$, with the shift being along the anti-diagonal direction in the *ω*, *ω*′ plane. It is not the same as the anti-diagonal projection *K*_{coh}(*x*) that appears in the analogous formula [Eq. (99)] for the DQC pathway.

### F. Exponential coherent state pulse: NRP and RP pathways

Here, we consider under what conditions the NRP and RP contributions to the *f*-state probability cancel exactly or approximately. We use Eqs. (87) and (119) to write their sum as

where

Thus, terms in which two intermediate states are degenerate ($\omega ee\u2032=0$) do not contribute.

To compare with the DQC term, again consider excitation by a coherent-state pulse in the form of a single-sided exponential, as in Eq. (103). Then, *G*_{coh}(*y*) = 4Γ^{2}/(*y*^{2} + 4Γ^{2}), which gives

Consider the special case in which the product of the four dipole matrix elements $Mfe\u2032eg$ is real, the plausibility of which is discussed in Appendix B. Then, the sum of NRP and RP contributions to the *f*-state probability depend only on the real part of $Re,e\u2032NRP+Re,e\u2032RP$,

Compare this result to the corresponding result for DQC when on TPA resonance, from Eq. (106),

The first factors of these two expressions are identical. The magnitude of the middle term in Eq. (128), $\omega ee\u2032/(\u2212\omega fe+\omega 0)$, is of the order of unity or much smaller. The third term has two limiting behaviors,

Although these expressions are rather complex, we conclude that the NRP + RP term can be comparable to the DQC term under some conditions. On the other hand, it will be much smaller than the DQC term if $\omega ee\u2032$ is much smaller in magnitude than (−*ω*_{fe} + *ω*_{0}) and $\gamma ee\u2032$ is comparable to or greater than *γ*_{fg}.

Going beyond the special case that $Mfe\u2032eg$ is real requires more complicated analysis, which we avoid in this Tutorial.

### G. Coherent-state rectangular pulse: NRP and RP pathways

Consider a rectangular coherent-state pulse exciting the NRP and RP pathways. To evaluate Eq. (121), we repeat the method of calculation in Sec. VII C, first using Eq. (108) in Eq. (123), which gives *G*_{coh}(*y*) = sin^{2}(*yT*/2)/(*yT*/2)^{2}. Then, we find

If the rectangular pulse is much longer than the inverse dephasing rate, $T\u226b1/\gamma ee\u2032$, we can neglect the second term, which is an oscillating turn-on transient.

## VIII. TWO-PHOTON EXCITATION BY ENTANGLED TWO-PHOTON STATES

We now come to a most interesting case, which is the primary motivation for this Tutorial—TPA with entangled photon pairs in the isolated-pair regime. As discussed in the Introduction, an important question is how large can the enhancement of TPA by time–frequency entanglement in the EPP state be? We expect such an enhancement if the photon pairs’ nondeterministic frequencies are correlated and sum to the TPA resonance frequency, even if each has an average spectrum that is much broader than the TPA transition linewidth. Again, we focus on the case that the intermediate molecular states are far from resonance. We find that large enhancement by time–frequency entanglement is possible.^{20,21} We show that the DQC pathway dominates, as we showed earlier for coherent state excitation under certain conditions.

### A. Two-photon-state DQC pathway

The probability to excite the *f* state via the far-off-resonance DQC pathway is given, in general, by Eqs. (87), (94), (96), and (101) as

where again $\omega \u0303\u2032\u2261\omega +\omega \u0303\u2212\omega \u2032$.

Changing variables to $x=(\omega +\omega \u0303\u22122\omega 0),z=\omega \u2212\omega 0$ (and *z*′ = *ω*′ − *ω*_{0}) analogously to the coherent-state case, we find

where

As for the coherent state result in Eq. (99), *K*_{ψ}(*x*) is an anti-diagonal projection of the JSA (or two-photon amplitude) onto the spectral region around the two-photon resonance. *K*_{ψ}(*x*) has the same form as the coherent-state DQC result with two replacements: $\varphi (\omega )\varphi (\omega \u0303)\u2192\Psi (\omega ,\omega \u0303)$ and *N*^{2} → 4*ɛ*^{2}. That is, the result for EPP resembles the coherent-state result but with a generalized spectral dependence and with linear instead of quadratic photon-flux dependence. While the coherent state takes the form of a separable two-photon state without frequency correlations, the EPP state is a pure state that can be nonseparable, that is, $\Psi (\omega ,\omega \u0303)\u2260\varphi (\omega )\varphi (\omega \u0303)$, indicating frequency correlations (entanglement). In particular, the JSA may be much narrower in the direction of the diagonal frequency axis than the anti-diagonal axis, as shown in Fig. 8. In this case, the frequencies of the two photons are anticorrelated and, because the state is a pure state, entangled. Such correlations can enhance the rate. While the coherent state relies on “accidental” coincidences for photons to arrive together at the molecule, the EPP state has photons always arriving in pairs, giving linear scaling with flux (proportional to *ɛ*^{2}), as discussed in Sec. VI E.

We consider two cases.

### B. Separable two-photon state: DQC pathway

An ultrashort pump pulse together with a particular phase-matching condition of type-I SPDC can create a separable (factorable) JSA, as in Ref. 55. In this case, there is no spectral entanglement. Then, because of symmetry, $\psi (\omega ,\omega \u0303)=\varphi 0(\omega )\varphi 0(\omega \u0303)$, which gives

where *K*_{sep}(*x*) is given by

This expression has precisely the same spectral dependence as for the coherent-state pulse in Eq. (100), and thus, in this case, although the state is “non-classical,” there is no enhancement of TPA through spectral correlation. There can still be enhancement by photon number correlation.

### C. Anti-diagonally separable two-photon states: DQC pathway

An important example of time–frequency entanglement is EPP generated by type-I SPDC using a narrow-band pump pulse with duration *T* that is long compared to the inverse of the phase-matching bandwidth. Energy conservation localizes the JSA $\Psi (\omega ,\omega \u0303)$ along the antidiagonal, $\omega \u0303+\omega =\omega P$. Properly designed, such a state maximizes the enhancement of TPA by time–frequency EPP.^{20,21}

For degenerate SPDC, the JSA has a single peak at the frequency where both $\omega ,\omega \u0303$ equal *ω*_{0} = *ω*_{P}/2. In the absence of dispersion, which creates phase correlations, we can model the JSA as the product of narrow and broad functions, *ψ*_{N}(*ω*) and *ψ*_{B}(*ω*) respectively, centered at $\omega =\omega \u0303=\omega P/2$ and oriented along diagonal and antidiagonal axes in the $(\omega ,\omega \u0303)$ plane, respectively. The width of *ψ*_{N}(*ω*) is the linewidth (inverse duration) of the pump pulse, and the linewidth of *ψ*_{B}(*ω*) is roughly $2$ times the spectral width of the EPP, set by the phase-matching conditions. Then,

where by state symmetry it is required that *ψ*_{B}(−*x*) = *ψ*_{B}(*x*), and both functions are square-normalized in *đ**x* = d*x*/2*π*. Note that while we call such a state “anti-diagonally separable,” it is an entangled state with regard to the frequencies of the two photons.

Then, we find

This gives the probability

This equation shows that it is the spectrally narrow function *ψ*_{N}(*x*) that effectively drives the molecular transition as a result of the anticorrelation of EPP frequencies.

### D. Gaussian EPP pulse: DQC pathway and EPP enhancement

The above integral can be evaluated by assuming Gaussian forms (valid for type-0 or type-I SPDC in the case of a long narrow-band pump pulse)

where we impose the condition *σ*_{N} < *σ*_{B}. The “narrow” width *σ*_{N} equals the spectral width of the laser pulse driving the SPDC, while the “broad” width *σ*_{B} is determined by phase matching.^{55} In this case,

which leads to the probability

where *ξ*(*z*) = exp(+*z*^{2})erfc(*z*), which we also encountered for the Gaussian shaped coherent-state pulse in Sec. VII D. A related expression was given in Ref. 6. In two limits, this becomes

The factor *σ*_{B}/*γ*_{fg} in the first of the two expressions in Eq. (143) can be interpreted as the number of temporal modes in the EPP state that impinge on the molecule in the molecular coherence time 1/*γ*_{fg}, which in this case is much smaller than the EPP pulse duration $1/\gamma fg\u226a1/\sigma N$. Comparing this result with the comparable one in the first line of Eq. (118), we see that the EPP result differs from the coherent-state result by a factor roughly equal to (4*ɛ*^{2}/*N*^{2})(*σ*_{B}/*σ*). The first of these factors is the ratio of photon fluxes and the second is the ratio of the bandwidths of the excitation pulses. Consider that the photon fluxes are equal, which requires them to be small enough that in the case of EPP, there are no overlapping pairs on average. Then, we see that the factor (*σ*_{B}/*σ*) can be much greater than one if the EPP is generated in a broad band containing many time–frequency modes. In contrast, the coherent state pulse is a single time–frequency mode, which contains, in this case, a single pair of photons on average.

The second of these two expressions, for *γ*_{fg} ≪ *σ*_{N}, is in the impulsive limit with respect to the correlation duration (∼1/*σ*_{B}) of the EPP wave packet. In this Gaussian model of the JSA, the ratio *σ*_{B}/*σ*_{N} is a measure of the number of temporal (time–frequency) modes in the EPP state that impinge on the molecule in a single pulse with duration ∼1/*σ*_{N} and thus is a measure of entanglement.^{56} Given that the EPP always arrive together within a time 1/*σ*_{B}, regardless of the duration (1/*σ*_{N}) of the pulse, the probability for TPA is enhanced by this factor relative to a narrow-band coherent pulse of the same duration, wherein the photons arrive independently.

We note that the Gaussian approximation can also be applied to type-II SPDC, with similar conclusions for the long-pulse case.

### E. Exponential EPP pulse: DQC pathway

For later comparison with the NRP and RP pathways, consider modeling the DQC term with an EPP pulse spectrum given by a complex Lorentzian for the narrow function and a Gaussian for the broad function. In the spectral domain, we have

This model has a single-sided exponential pulse shape in the time domain, as in Eq. (103). Then,

and

On TPA resonance (2*ω*_{0} = *ω*_{fg}), this gives, upon writing Σ^{(2)} explicitly,

This result has limits similar to those in Eq. (143), with Γ and *σ*_{B} playing similar roles. The explanation of the potentially very large enhancement by EPP spectral correlations are similar in this case, showing that enhancement occurs for a variety of pulse shapes.

### F. NRP and RP pathways excited by two-photon states

Here, we address whether or not the NRP and RP pathways are significantly enhanced by time–frequency entanglement of the exciting photon pairs. One might expect the answer to be no because these pathways are impacted by dephasing processes that could disrupt the delicate anticorrelation of the photon pair’s frequencies.

where as in Eq. (121),

and using Eq. (94),

We wish to compare the NRP and RP contributions Eq. (148) to the DQC contribution Eq. (139) in the case that the JSA equals the product of narrow and broad functions, *ψ*_{N}(*ω*) and *ψ*_{B}(*ω*) respectively, centered at $\omega =\omega \u0303=\omega P/2$ and oriented along diagonal and antidiagonal axes in the $(\omega ,\omega \u0303)$ plane, as in Eq. (137). Then, changing variables to $x=\omega +\omega \u0303\u22122\omega 0\u2261\omega D\u22122\omega 0,z=(\omega \u2212\omega 0)/2\u2261\omega A$ gives for G_{Ψ}(y),

The narrow function has dropped out, meaning the NRP and RP terms are independent of pulse duration in this scenario, which includes a single EPP with the entanglement time determined by the broad function. We see that the autocorrelation of the two-dimensional function $\Psi coh(\omega ,\omega \u0303)$ along the anti-diagonal direction in the *ω*, *ω*′ plane [as in Eq. (124)] has been converted to an autocorrelation of a one-dimensional function *ψ*_{B}(*ω*_{A}). Again, it is not the same as the anti-diagonal projection *K*_{coh}(*x*) that appears in the analogous formula, Eq. (99), for the DQC pathway.

The interpretation of Eq. (149) with Eq. (151) and the accompanying illustration in Fig. 9 is that the population damping rate (*γ*_{ee}) of the intermediate state or the dephasing rate ($\gamma ee\u2032$) of the intermediate-state coherence acts as a filter. This filter limits or “windows” the range of frequencies in the exciting field that contribute to exciting TPA via the NRP and RP pathways through the $\omega ee\u2032$ coherence. In contrast, for the DQC pathway, illustrated in Fig. 8, there is no such windowing behavior because the pathway bypasses the intermediate-state populations and coherences.

In contrast to the DQC contribution, which is dominated by the narrow part of the JSA, the NRP and RP contribution are determined by the broad part of the JSA. For this reason, we do not expect enhancement of the NRP and RP processes by spectral correlation.

To find the combined contribution of the NRP and RP probabilities, we sum Eq. (148) and evaluate the result at the TPA resonance frequency *ω*_{0} = *ω*_{fg}/2 to give

To achieve a direct comparison with the DQC term in a particular scenario, consider again modeling the EPP pulse spectrum by a complex Lorentzian for the narrow function and a Gaussian for the broad function, as in Eq. (144). Then, the narrow function drops out and the broad function determines *G*_{Ψ}(*y*) in Eq. (151) to be $G\Psi (y)=4\epsilon 2\u2061exp\u2212y2/8\sigma B2$, leading to

Rather than attempting this integral in closed form, consider two limits of EPP bandwidth relative to $\gamma ee\u2032$,

We write the sum of probabilities as

If $Mfe\u2032eg$ is real, then

Observing that $\gamma ee\u20322/(\gamma ee\u20322+\omega ee\u20322)\u22641$ and that the factor $\omega ee\u2032/(\u2212\omega fe+\omega 0)$ has a magnitude of order unity or less, we can compare this result in either limit with Eq. (147) for that of the DQC pathway under the same form of EPP excitation field. We see that the DQC-pathway probability is greater than that for the NRP + RP pathways by at least a factor roughly equal to

which can be much greater than unity if the number of temporal modes, *σ*_{B}/Γ, in the EPP pulse is large.

To emphasize this most important conclusion, we have shown that whereas the DQC pathway can be greatly enhanced relative to the coherent-state case, as in Eqs. (143) and (147), the NRP + RP pathways collectively are not enhanced by the spectral anticorrelations (time–frequency entanglement) of the EPP field. Thus, for EPP excitation, the DQC is predicted to be dominant.

### G. Effect of dispersion on TPA by EPP

The frequency-dependent refractive index of optical elements that light passes through before reaching the sample is known to broaden ultrashort pulses temporally and reduce their effectiveness in nonlinear-optical processes. To account for such effects in EPP-driven TPA, we incorporate dispersive propagation into the two-photon JSA by replacing^{20}

where *D* is the second-order (group delay) dispersion of the transmitting optical system.

As an example, we insert this expression into Eq. (132) for the DQC term excited by EPP. Then, using the Gaussian forms in Eq. (140) (valid for long SPDC pump) leads to

Comparing to Eq. (141), we see the sole effect of second-order dispersion is to replace *σ*_{B} by

The *x* dependence of *K*_{ψ}(*x*) is not altered; there is only an overall decrease in the magnitude of $K\psi x$ and thus a decrease in the TPA probability.

## IX. EXPERIMENTAL CHALLENGES

Regarding the implementation of experiments on ETPA, it is important to consider three issues: (1) What are the experimental signatures that can provide indisputable evidence for ETPA? (2) What are potential reasons one might miss observing ETPA accidentally? (3) What are possible reasons classical signals might be misidentified as evidence for ETPA? Here, we follow, in part, the outline given in the Supplementary data of Ref. 17, which is based, in part, on Refs. 18 and 20.

Experimental signatures that can provide indisputable evidence for ETPA include a combination of (a) linear scaling with optical flux incident on the sample (but this alone is not sufficient as several classical processes can mimic this), (b) quadratic dependence on linear loss between the SPDC source and the sample, and (c) experimental verification that the flux being measured at the sample consists of photon pairs (by coincidence counting of photons in the sample volume when the sample is removed).

Potential reasons one might miss observing ETPA accidentally include (a) insufficient spatial overlap of photon pairs in the sample, (b) linear dispersion that broadens the EPP correlation time and reduces the effectiveness of ETPA, (c) detector saturation and dead-time effects, (d) insufficient EPP flux and/or fluorescence collection efficiency, and (e) reabsorption of fluorescence in the molecular sample.

Possible reasons classical signals might be misidentified as evidence for ETPA include (a) observation of linear scaling of signal with optical flux incident on the sample without performing the other checks mentioned above and (b) the presence of low-lying resonant intermediate states not recognized for the sample being used, for example, such states created by molecular aggregation.

## X. SUMMARY AND DISCUSSION

Given the challenges cited in Sec. IX, the scientific community that is working to develop ETPA as a tool for quantum-enhanced spectroscopy and imaging of molecular and biological samples is still struggling to identify the techniques and conditions under which such a “quantum advantage” can be achieved. While several experiments have presented evidence that ETPA in molecules does provide such advantages, other experiments, as well as the theory summarized here, have called those conclusions into question.

We have reviewed the quantum optics theory needed for incorporating entangled quantum states of light into the theory of two-photon absorption by atoms or molecules. The density matrix (or operator) in fourth-order perturbation theory plays a central role because it can describe damping and dephasing of the states and their mutual coherences that contribute to the TPA process. This method is in contrast with the conventional second-order perturbation theory that uses state amplitudes and includes transition linewidths only by averaging over the final density of states. (Note that averaging over a density of states is equivalent to homogeneous broadening in lowest-order perturbation theory where there is no saturation of populations.) The conventional theory corresponds to the double quantum coherence (DQC) pathway (double-side Feynman diagram), while the additional pathways included in the density-matrix approach are the nonrephasing (NRP) and rephasing (RP) pathways.

Our treatment clarifies to what extent the predictions of these two approaches differ. We find that if the exciting field is far from resonance with any intermediate states, the conventional DQC pathway typically dominates the TPA process, although under some conditions, the NRP and RP pathways can make significant contributions.

The treatment we developed confirms that the DQC contribution can be greatly enhanced by the presence of frequency anticorrelations in the exciting field composed of photon pairs created in, for example, spontaneous parametric down conversion. The enhancement occurs because the frequencies of the two photons sum to that of the SPDC pump laser, so if that laser has a narrow bandwidth, the sum-frequency variable is “compressed” into the TPA transition line profile. An equivalent explanation can be given in the time domain: the entanglement time of the photon pairs (inverse bandwidth) can be much shorter than the overall duration of the illumination pulse, meaning that the pair behaves as if confined to an ultrashort pulse whose arrival time is indeterminate. Because TPA is a nonlinear-optical, two-photon process, it is enhanced when the exciting light is confined to shorter time intervals. Detailed discussion and quantitative estimates of such effects are given in Ref. 20, and a general derivation of a hard upper bound of such enhancement is given in Ref. 21.

The enhancement of TPA probability by EPP relative to that using a coherent state with a similar pulse shapes can be quantified by, for example, combining results from Eqs. (107) and (146). We see that assuming an exponential pulse shape, the DQC contribution driven by EPP is enhanced relative to the coherent-state case by the ratio, which we call the *quantum enhancement factor*^{20,21}

where we used *N*_{EPP} = 2*ɛ*^{2}. Recall that *σ*_{B} is roughly the bandwidth of the EPP field, while Γ is roughly the inverse of the pulse duration for both the EPP and coherent-state pulses. Similar results are found for various pulse shapes, and a general result independent of pulse shape is given in Ref. 21. Recall the EPP results are valid under the condition of isolated EPP, that is, not overlapping photon pairs, and are therefore restricted to very low flux.

When the mean number of EPP photons in an excitation pulse is equal to the mean number of photons in a weak coherent-state pulse, *N*_{EPP} = 2*ɛ*^{2} = *N*_{coh} ≡ *N*_{both}, this ratio equals (2/*N*_{both}) × (*σ*_{B}/Γ). Thus, if the mean number of photons in a pulse is much less than one, the first factor, which we call the “photon-number enhancement factor,” can be large. The second factor, which we call the “spectral enhancement factor,” can also be large if the EPP field’s bandwidth is much greater than its inverse duration. This condition corresponds to large time–frequency entanglement.

Unfortunately, in many applications of interest, such as spectroscopy or two-photon microscopy, the predictions here indicate that in most practical cases, the predicted final-state population is too small to be detected for typical molecules using typical technology in current use. As discussed in Sec. VII C, conventional two-photon cross sections are extremely small typically. The amount of enhancement that can be achieved by the number- and frequency correlations calculated here is not likely great enough to overcome the small cross section in typical scenarios.

A major question that is addressed by the theory is to what extent the NRP and RP contributions are similarly enhanced by frequency anticorrelations in the exciting field. We find, not surprisingly, that there is no such enhancement because these processes are step-wise, occurring through populations or mutual coherences among intermediate states. The step-wise nature of these processes disrupts the delicate frequency anticorrelations in the exciting field leading to no enhancement. This Tutorial is the first, to our knowledge, pointing out that the NRP and RP pathways do not provide an explanation for the anomalously large ETPA probabilities reported in some experimental studies.

This conclusion helps address the presently controversial issue discussed above concerning the detectability of entangled two-photon absorption (ETPA). Some experimental studies^{15,16} have reported apparent values of ETPA excitation probabilities that greatly exceed the values predicted here [for example, Eq. (142)]. Other studies have recently found upper bounds on EPTA probabilities that are much smaller and in line with the present predictions,^{17–19} but the final word has yet to be spoken on this question.

In the following, we comment on limitations and possible extensions of the theory.

First, we note that the fourth-order perturbation solution for the density matrix that is used here and in many treatments of ultrafast-laser spectroscopy is suited for excitation by short laser pulses but not applicable to calculate steady-state responses. This theory treats the final state as merely an “integrating receptacle” for population. The challenge is to develop a non-perturbative treatment that allows the final state population to be a dynamical variable. This would allow steady-state solutions, including quantum states of the exciting field, and would enable direct comparisons with semiclassical treatments, such as presented in Refs. 50–54.

Another limitation of the present theory is that it treats only isolated photon pairs incident on the atom or molecule. Some experiments have been carried out with large fluxes such that pairs do overlap,^{14,18} and a couple of theoretical treatments cover such cases.^{6,23} The study of TPA using multi-spectral-mode squeezed states of light would be worthwhile, and careful consideration of the roles of NRP and RP pathways in this context would be of interest.

A challenging question is to what extent the frequency correlations that enhance ETPA could be mimicked by “classical” fluctuating fields. Frequency correlations could be built into such a model using a statistical mixture of coherent states, as in Ref. 57, but Schlawin and Buchleitner argued that such states do not enhance the absorption probability above the pure coherent-state case,^{24} and Lerch and Stefanov showed that a statistical mixture of correlated monochromatic states can mimic the frequency correlations but not the time correlations.^{58} Furthermore, the enhancement by photon number correlations (the correlated arrival of pairs) likely cannot be mimicked perfectly by classical fields.

The treatment of collisional dephasing used here models homogeneous line shapes as Lorentzian at all frequencies. While this approach is standard and common, it is an oversimplification of the physics. There are two well-known ways to improve the treatment—the Brownian-oscillator-bath model, which is appropriate for spectroscopy of solvated molecules,^{59} and the non-impact theory of collisional line broadening, which is well developed for spectroscopy of atomic or molecular vapors.^{60,61} For example, the effectiveness of collisional dephasing can be greatly reduced when light is detuned far from the line center. To invert the argument, observing excitation of states far from resonance can be used to characterize the dephasing bath itself and yield important information.

Finally, we mention the possibility to use optical phase modulation techniques^{62–65} or phase-cycling techniques^{66} similar to those used in multidimensional spectroscopy with ultrafast lasers to dissect the separate contributions of the DQC, NRP, and RP pathways. If the challenge of low signal levels can be overcome, then combining EPP excitation with such phase modulation techniques might provide a new avenue for obtaining hard-to-get information on molecular structure and dynamics.

## ACKNOWLEDGMENTS

We thank Perry Rice, Paul Berman, Jeff Cina, Daniel Steck, Steven van Enk, Markus Allgaier, Brian Smith, and Sofiane Merkouche for helpful discussions. This work was supported by a grant from the National Science Foundation RAISE-TAQS Program (Grant No. PHY-1839216).

The authors declare no conflicts of interest.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: MOLECULAR DEPHASING AND CORRELATION FUNCTIONS

The molecular correlation functions introduced in Eq. (73) are evaluated by first writing them in a common form,

where

which gives in the three cases,

where

To apply Kubo dephasing theory (Sec. V B), we need to group the factors, so we can identify disjoint time intervals in which the dephasing interactions occur. This grouping allows treating the dephasing interactions in one interval as statistically independent of those in other intervals. Thus, we introduce the difference-time variables *r* = *t*_{4} − *t*_{3}, *s* = *t*_{3} − *t*_{2}, *τ* = *t*_{2} − *t*_{1}, which means *t*_{1} = (*t*_{4} − *r* − *s* − *τ*), *t*_{2} = (*t*_{4} − *r* − *s*), *t*_{3} = (*t*_{4} − *r*). The integration ranges are [0, *∞*] for *r*, *s*, and *τ* and [−*∞*, *∞*] for *t*_{4}. Then, denoting *ω*_{ij} = *ω*_{i} − *ω*_{j}, we have

Note *t*_{4} has dropped out in these terms. Since the time intervals are disjoint, we can apply Kubo theory in each separately and replace

Then, from Eq. (A3),

which reproduces Eq. (74).

### APPENDIX B: UNDER WHAT CONDITIONS IS THE DIPOLE MATRIX ELEMENT PRODUCT REAL?

In cases where the dipole matrix element product $Mfe\u2032eg=\mu fe\u2032\mu e\u2032g\mu ge\mu ef$ is real, the analysis of expressions like those in Eq. (87) is simplified. Such a simplification is especially useful for the NRP and RP terms, where $Mfe\u2032eg$ does not separate from the frequency-dependent functions as it does for the DQC case, as shown in Eq. (100). While the question is $Mfe\u2032eg$ real is not easily answered in general, it is useful to consider a few special cases.

In one-dimensional systems, all eigenfunctions are real, and thus, $Mfe\u2032eg$ is real.

For TPA in single-electron atoms, the only complex variations in the eigenfunctions enter in the form exp(−*imϕ*), and these factors when integrated will always produce real matrix elements.

In molecules with a high degree of symmetry, $Mfe\u2032eg$ can be proven to be real, at least in special cases. Consider, for example, a symmetric *N*-mer that is *N* identical atoms or molecules coming together to form a symmetric structure. Denote the collective ground state, with all monomers in their lowest-energy state by $g$. Denote the set of singly excited states, with one monomer in its first excited state, by $un$ (*n* = 1 to *N*), and denote a particular doubly excited state, with two monomers in their first excited state, by $f$. Assume by symmetry that all dipole matrix elements connected to the ground state are equal: $un\mu g=\mu g$, and assume that all dipole matrix elements connected to the doubly excited state are equal: $un\mu f=\mu f$. If the degeneracy of the singly excited states is lifted by a symmetric interaction among the systems, the singly excited states are mixed by a unitary transformation to create new eigenstates in the singly excited manifold,

The dipole matrix elements are transformed to

Then,

Conjugating this expression and swapping indices *n* ↔ *m*, *r* ↔ *s* shows $Mfe\u2032eg*=Mfe\u2032eg$, which verifies it is real.

We are not aware of more general proofs on the reality of $Mfe\u2032eg$, but such would be useful for predicting TPA signals in complex molecules.

### APPENDIX C: TPA BY A COHERENT PULSE IN THE IMPULSIVE LIMIT

In an extreme limit, if the coherent state is a pulse much shorter than the molecular dephasing time—the impulsive limit—then Eq. (100) becomes

where we used *K*(0) = *∫A*^{2}(*t*)d*t*. If *A*(*t*) is real and positive, then *K*(0) = 1, reproducing the result in Eq. (118) for the Gaussian coherent state in the impulsive limit *σ* ≫ *γ*_{fg}. In this limit, the probability does not depend on *γ*_{fg} or *σ* because for the nonlinear TPA process, the effect of spreading the spectrum over a broader range is compensated by the increasing intensity in the time domain.

Note that *A*^{2}(*t*) is proportional to the two-photon Rabi frequency^{67} and can be complex. This means that the TPA probability can go to zero in the impulsive limit if the pulse constitutes a two-photon zero-*π* pulse, defined by *K*(0) = 0.