# How to optimize the absorption of two entangled photons

Edoardo G Carnio1,2, Andreas Buchleitner1,2 and Frank Schlawin3,4,5,6 1
Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, Hermann-Herder-
Straße 3, D-79104, Freiburg, Germany  2 EUCOR Centre for Quantum Science and
Quantum Computing, Albert-Ludwigs-Universität Freiburg, Hermann-Herder-Straße
3, D-79104, Freiburg, Germany  3 Clarendon Laboratory, University of Oxford,
Parks Road, Oxford OX1 3PU, United Kingdom  4 The Hamburg Centre for Ultrafast
Imaging, Luruper Chaussee 149, Hamburg D-22761, Germany  5 Max Planck
Institute for the Structure and Dynamics of Matter, Luruper Chaussee 149,
D-22761 Hamburg, Germany  6 Universität Hamburg, Luruper Chaussee 149, D-22761
Hamburg, Germany

(August 27, 2024)

###### Abstract

We investigate how entanglement can enhance two-photon absorption in a three-
level system. First, we employ the Schmidt decomposition to determine the
entanglement properties of the optimal two-photon state to drive such a
transition, and the maximum enhancement which can be achieved in comparison to
the optimal classical pulse. We then adapt the optimization problem to
realistic experimental constraints, where photon pairs from a down-conversion
source are manipulated by local operations such as spatial light modulators.
We derive optimal pulse shaping functions to enhance the absorption
efficiency, and compare the maximal enhancement achievable by entanglement to
the yield of optimally shaped, separable pulses.

##  1 Introduction

Coherent control generally refers to the manipulation of quantum dynamical
processes by suitably shaped control fields or interactions [1, 2, 3]. It has
found widespread application in the control of chemical reactions [4, 5, 6, 7,
8, 9], in spectroscopy [10, 11], in laser cooling [12, 13, 14] or even in
quantum information and computing [15, 16]. In these well-established
scenarios, the control is built on the manipulation of interfering pathways to
maximize a desired outcome or speed up a process using classical laser pulses
[17]. Optimal control theory then consists in finding the optimal laser pulse
shapes and sequences, and the quantum character of light can be safely
neglected.

In recent years, the possible use of quantum light sources for quantum-
enhanced applications in microscopy or spectroscopy has gathered a lot of
attention [18, 19, 20, 21, 22, 23]. Of particular interest is the use of
entangled photon pairs for applications which involve two-photon transitions.
Such pairs can induce two-photon transitions more efficiently than laser
pulses, and promise the use at low photon flux, thus preventing damage in
photosensitive samples [24, 25, 26, 27]. In addition, quantum correlations may
help to further manipulate optical signals [28, 29, 30]. Yet, a key problem in
the practical application, currently, is the low absorption cross section of
many samples [31, 32, 33, 34]. One possible remedy to address this issue and
further this field is the shaping of entangled photonic wave functions [35,
36, 37], to enhance the absorption probability. This strategy was explored in
a number of recent theoretical papers, where the crucial role of quantum
correlations between different travelling modes [38, 39, 40], or of the
quantum statistics of a cavity mode [41], was highlighted. In contrast to
classical control described by optimal control theory, the light fields have
to be treated quantum mechanically, and, due to the small photon number per
mode, perturbation theory can be employed.

Here we present a detailed study of how entanglement shared between two
photons can enhance the probability to induce a two-photon excitation, in a
sample which we describe by a simple three-level toy model, with finite
excited state lifetimes, as depicted in fig. 1(a).

Figure 1: (a) A two-photon transition from ground state $\ket{g}$ to target
state $\ket{f}$, via intermediate state $\ket{e}$, is driven by two incoming,
near-resonant single-photon pulses with carrier frequencies $\omega_{1}$ and
$\omega_{2}$, respectively. The excitation efficiency is subject to finite
decay rates of $\ket{e}$ and $\ket{f}$, the characteristic scale $\omega_{e}$
of the level spacing, the detuning $\Delta$ of the individual nearest
neighbour level spacings from degeneracy (in units of the intermediate state’s
lifetime $\gamma_{e}$), the deviation of the excited state decay rates’ ratio
from two (again in units of $\gamma_{e}$), see (17), and, here of central
interest, the frequency correlations (see panel (b)) inscribed into the
incoming pulses. (b) Density
$|T_{t}(\omega_{1},\omega_{2})|^{2}\propto|\Phi(\omega_{1},\omega_{2})|^{2}$
of the optimal two-photon wave function (31) in the space of rescaled
frequencies $(\omega_{j}-\omega_{e})/\gamma_{e}$, and parametrised by detuning
$\Delta=5$ and deviation $\delta=-1.9$ as specified in panel (a) and (17). (c)
Single-photon density (33) derived from (31) upon tracing over the frequency
of the partner photon in $|\Phi(\omega_{1},\omega_{2})|^{2}$, for the same
parameter values as in panel (b). The double-peaked distributions in (b,c)
feature maxima in the vicinity of the first ($\ket{g}\rightarrow\ket{e}$) and
second ($\ket{e}\rightarrow\ket{f}$) bare transition’s (see panel (a))
rescaled frequencies $0$ and $\Delta$, respectively, with distinct widths
determined by $\Delta$ and $\delta$, see (32,33). The widths of the local
maxima give a qualitative impression of the nonclassical frequency correlation
between both incoming photons.

We carefully examine the optimal pulse shapes, as well as the relation between
quantum correlations and the achievable enhancement. In a second part, we then
demonstrate how this formalism can be applied to derive optimal states under
realistic experimental conditions. In particular, we ask how a given entangled
two-photon state can be optimized by only local operations on the individual
photons, in order to induce the desired two-photon transition most
efficiently. We find that, depending on the initial two-photon state and the
sample, substantial enhancements of the absorption probability can be
achieved.

The paper is organized as follows: In section 2, we formulate the coherent
control problem. We derive and solve functionals for the optimal quantum
states of light to excite a simple three-level quantum system. In section 3,
we describe properties of ideal pulses, and analyze the relation between the
entanglement they encode and the possible quantum advantage they offer (over
optimal separable quantum states of light). Subsequently, in section 4, we
turn to the question of how entangled photons can be optimized for two-photon
absorption (TPA) in experimentally realistic scenarios, before we conclude in
section 5.

##  2 Theoretical framework

###  2.1 The Hamiltonian

We consider the interaction between two (propagating) pulses of the quantized
radiation field (the “field” degrees of freedom) and three electronic energy
levels of an atom or molecule (the “matter” degrees of freedom). These systems
are modelled, respectively, by Hamiltonians $H_{\mathrm{f}}$ and
$H_{\mathrm{m}}$, and are coupled by an interaction term $W$. We describe all
of these terms in the following paragraphs.

The incoming light pulses impinge on the three-level target located at the
origin of our reference frame. Each field is quantized within a (cylindrical)
quantization volume, of cross section $A$, along a distinct propagation
direction, giving rise to modes labelled by a one-dimensional continuous
variable, either the wave vector $k$ or the frequency $\omega$ [42]. Choosing
the latter, for each beam $j$ we obtain annihilation $a_{j}(\omega)$ and
creation $a_{j}^{\dagger}(\omega)$ operators satisfying the commutation
relation
$[a_{j}(\omega),a_{l}^{\dagger}(\omega^{\prime})]=\delta_{jl}\delta(\omega-\omega^{\prime})$.
In this framework the photon number operator for beam $j$ reads
$n_{j}=\int_{0}^{\infty}\rmd\omega\,a_{j}^{\dagger}(\omega)a_{j}(\omega)$,
while the total Hamiltonian for both fields is
$H_{\mathrm{f}}=\sum_{j}\int_{0}^{\infty}\rmd\omega\,\hbar\omega
a^{\dagger}_{j}(\omega)a_{j}(\omega)$. At the origin, where it interacts with
the sample, the positive-frequency part of the electric field operator for the
Hilbert space of photon $j$ (in the interaction picture starting at $t_{0}$
with respect to $H_{\mathrm{f}}$) reads

$E^{+}_{j}(t)=\rmi\int_{0}^{\infty}\rmd\omega\left(\frac{\hbar\omega}{4\pi\epsilon_{0}cA}\right)^{1/2}a_{j}(\omega)\rme^{-\rmi\omega(t-t_{0})},$
(1)

where $\epsilon_{0}$ is the dielectric constant and $c$ the vacuum speed of
light. We assume the fields have parallel polarization, but are distinguished
by their propagation directions.

In the following we consider incoming photon pulses of bandwidth
$\Delta\omega$ much smaller than their central frequency $\omega_{0}$, i.e.
$\Delta\omega\ll\omega_{0}$. We can then employ the _narrow bandwidth
approximation_ , where we extend the range of all above frequency integrals to
$(-\infty,\infty)$, substitute $\omega$ with $\omega_{0}$ within the square
root in the integrand of (1), pull the resulting constant factor out of the
integral, and define the Fourier-transformed annihilation operator

$a_{j}(t)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\rmd\omega\,a_{j}(\omega)\rme^{-\rmi\omega(t-t_{0})}.$
(2)

Together with its adjoint $a_{j}^{\dagger}(t)$, it satisfies the commutation
relation
$[a_{j}(t),a_{l}^{\dagger}(t^{\prime})]=\delta_{jl}\delta(t-t^{\prime})$. The
electric field operator (1) for beam $j$ is finally re-expressed as

$E_{j}^{+}(t)=\rmi\mathcal{E}_{0}a_{j}(t)\,,\,{\rm
with}\;\mathcal{E}_{0}=\sqrt{\hbar\omega_{0}/(2\epsilon_{0}cA)}\,,$ (3)

so that the total (positive-frequency) electric field seen by the atomic
target reads (we neglect any geometry dependence in the coupling factors [23])

$E^{+}(t)=E_{1}^{+}(t)+E_{2}^{+}(t)=\rmi\mathcal{E}_{0}\left[a_{1}(t)+a_{2}(t)\right].$
(4)

As shown in fig. 1(a), the three non-degenerate electronic energy eigenstates
of our target are $\ket{g}$, $\ket{e}$, and $\ket{f}$, with increasing energy.
We define the origin of the energy axis to coincide with the energy of
$\ket{g}$, hence the excited state energies are $\hbar\omega_{e}$ and
$\hbar\omega_{f}$, respectively. With these labels, the matter Hamiltonian is
$H_{\mathrm{m}}=\hbar\omega_{e}\ket{e}\bra{e}+\hbar\omega_{f}\ket{f}\bra{f}$.

The light-matter coupling is mediated by an electric dipole Hamiltonian,
switched on at $t_{0}$, which – in the interaction picture with respect to
$H_{\mathrm{f}}+H_{\mathrm{m}}$, and for near-resonant perturbative driving of
the atomic transitions – reads

$W_{I}(t)=-V(t)E^{-}(t)-V^{\dagger}(t)E^{+}(t),$ (5)

with $E^{-}$ the adjoint of $E^{+}$. For the specific level structure here
considered (including the assumption that the atomic eigenstates have well-
defined parity), the dipole operator, along the fields’ polarization, has the
explicit form

$V(t)=\mu_{ge}\rme^{-\rmi\omega_{e}(t-t_{0})}\ket{g}\bra{e}+\mu_{ef}\rme^{-\rmi(\omega_{f}-\omega_{e})(t-t_{0})}\ket{e}\bra{f}.$
(6)

where the dipole matrix transition elements $\mu_{ge}$ and $\mu_{ef}$ between
$\ket{g}$ and $\ket{e}$, and between $\ket{e}$ and $\ket{f}$, respectively,
can be chosen real-valued.

###  2.2 Two-photon absorption amplitude

Our objective is to identify the optimal two-photon field state $\ket{\Phi}$
that drives the matter degrees of freedom from its initial (at time $t_{0}$)
state $\ket{g}$ into the target state $\ket{f}$ (at time $t$), by TPA via
$\ket{e}$. This is tantamount of maximizing the transition probability

$p_{f}(t)=|\braket{f(t),0}{U_{I}(t)}{g,\Phi}|^{2},$ (7)

where $\ket{0}$ indicates the vacuum state of both injected fields,
$\ket{f(t)}=\rme^{\rmi\omega_{f}(t-t_{0})}\ket{f}$ exhibits the explicit time
dependence of the interaction picture with respect to $H_{\mathrm{m}}$
(whereas $\ket{g}$, at energy zero, remains unaffected), and $U_{I}(t)$ is the
time evolution operator in the interaction picture of
$H_{\mathrm{f}}+H_{\mathrm{m}}$, given by the Dyson series

$U_{I}(t)=\mathds{1}+\sum_{n=1}^{\infty}\left(\frac{1}{\rmi\hbar}\right)^{n}\int^{t}_{t_{0}}\rmd\tau_{n}\int^{\tau_{n}}_{t_{0}}\rmd\tau_{n-1}\ldots\int^{\tau_{2}}_{t_{0}}\rmd\tau_{1}W_{I}(\tau_{n})\ldots
W_{I}(\tau_{1}).$ (8)

Since $\ket{\Phi}$ is to be optimized in (7), while the initial atomic and
field, as well as the final atomic state are fixed, we can extract the
transition amplitude operator

$T_{fg}(t)=\braket{f(t)}{U_{I}(t)}{g}.$ (9)

acting solely on the field degrees of freedom.111Note that (7,9) imply a
strictly unitary dynamics, and that rigorous account of incoherent processes
beyond the purely phenomenological level adopted below would require a
treatment in terms of the matter-field density matrix [43].

###  2.3 Matter response function

Given the faint incoming field implied by (7) upon fixing the impinging two-
photon state $\ket{\Phi}$, we can expand $U_{I}(t)$ in powers of $W_{I}$, such
that the leading perturbative contribution to the desired two-photon
transition reads, with (5,6,9),

$T_{fg}(t)=\left(\frac{1}{\rmi\hbar}\right)^{2}\int^{t}_{t_{0}}\rmd{\tau_{2}}\int^{\tau_{2}}_{t_{0}}\rmd\tau_{1}\,\braket{f(t)}{W_{I}(\tau_{2})W_{I}(\tau_{1})}{g}.$
(10)

Using (4) and (5), the integrand becomes

$\displaystyle\braket{f(t)}{W_{I}(\tau_{2})W_{I}(\tau_{1})}{g}=$
$\displaystyle-\mathcal{E}_{0}^{2}\mu_{ge}\mu_{ef}\rme^{-\rmi\omega_{f}t}\rme^{\rmi(\omega_{f}-\omega_{e})\tau_{2}}\rme^{\rmi\omega_{e}\tau_{1}}\left[a_{1}(\tau_{1})a_{2}(\tau_{2})+a_{1}(\tau_{2})a_{2}(\tau_{1})\right],$
(11)

where we dropped the terms $a_{i}(\tau_{1})a_{i}(\tau_{2})$ which result from
$E^{+}(\tau_{2})E^{+}(\tau_{1})$ by (4), since the sought-after pulse
$\ket{\Phi}$ has only one photon in each mode. Equation (2) and an exchange of
the frequency and time integrations gives

$T_{fg}(t)=\int_{-\infty}^{\infty}\rmd\omega_{1}\int_{-\infty}^{\infty}\rmd\omega_{2}\,T_{t;t_{0}}(\omega_{1},\omega_{2})a_{2}(\omega_{2})a_{1}(\omega_{1}),$
(12)

where $T_{t;t_{0}}(\omega_{1},\omega_{2})$ is the _matter response function_
under weak field driving and for a finite interaction time starting at
$t_{0}$, as made explicit by the index. Note that
$T_{t;t_{0}}(\omega_{1},\omega_{2})$ must also encode the symmetry of (2.3)
under exchange of the photon from the first or the second pulse being absorbed
first. This will become explicit in (2.3).

To compactify the explicit expressions for
$T_{t;t_{0}}(\omega_{1},\omega_{2})$, we introduce the line shape functions

$\mathcal{L}_{s}(\omega)=\frac{\rmi\mathcal{E}_{0}}{\sqrt{2\pi}\hbar}\frac{\mu_{s-1\,s}}{\omega-\omega_{s}+\rmi\gamma_{s}},$
(13)

where the dipole matrix element $\mu_{s-1\,s}$ connects the state $s$ to the
next energetically lower lying state $s-1$ (e.g., $\ket{e}$ to $\ket{g}$), and
$\gamma_{s}$ phenomenologically accounts for finite decay rates of $\ket{f}$
and $\ket{e}$ ($\ket{g}$, being the ground state, cannot decay). With
(10,2.3,13) we can thus extract from (12)

$\displaystyle T_{t;t_{0}}(\omega_{1},\omega_{2})$
$\displaystyle=\frac{\mathcal{E}_{0}^{2}}{2\pi\hbar^{2}}\mu_{ge}\mu_{ef}\rme^{-\rmi(\omega_{f}-\rmi\gamma_{f})t}\rme^{\rmi(\omega_{1}+\omega_{2})t_{0}}$
$\displaystyle\times\int^{t}_{t_{0}}\rmd{\tau_{2}}\int^{\tau_{2}}_{t_{0}}\rmd{\tau_{1}}\,\rme^{-\rmi[\omega_{2}-(\omega_{f}-\omega_{e})+\rmi(\gamma_{f}-\gamma_{e})]\tau_{2}}\rme^{-\rmi(\omega_{1}-\omega_{e}+\rmi\gamma_{e})\tau_{1}}+\big{(}\omega_{1}\leftrightarrow\omega_{2}\big{)}$
$\displaystyle=\mathcal{L}_{e}(\omega_{1})\left\\{\left[\rme^{-\rmi(\omega_{1}+\omega_{2})(t-t_{0})}-\rme^{-\rmi(\omega_{f}-\rmi\gamma_{f})(t-t_{0})}\right]\mathcal{L}_{f}(\omega_{1}+\omega_{2})\right.$
$\displaystyle\left.-\left[\rme^{-\rmi(\omega_{2}+\omega_{e}-\rmi\gamma_{e})(t-t_{0})}-\rme^{-\rmi(\omega_{f}-\rmi\gamma_{f})(t-t_{0})}\right]\mathcal{L}_{f}(\omega_{2}+\omega_{e})\right\\}+\big{(}\omega_{1}\leftrightarrow\omega_{2}\big{)},$

where $(\omega_{1}\leftrightarrow\omega_{2})$ in the last line indicates
identical terms with the frequencies $\omega_{1}$ and $\omega_{2}$ of the two
photons exchanged.

###  2.4 Infinitely extended pulses

Equation (2.3) describes the matter response at time $t$ to a pulse switched
on at $t_{0}$, and thus vanishes for $t=t_{0}$. In the following, we will
consider a scenario where the light-matter interaction is always switched on,
while the case of finite $t-t_{0}$ will be discussed elsewhere. We therefore
take the limit $t-t_{0}\rightarrow\infty$, and, as a consequence, the real
exponential factors in (2.3) vanish due to the excited state lifetimes, and we
obtain the simpler expression [40]

$T_{t}(\omega_{1},\omega_{2})=\rme^{-\rmi(\omega_{1}+\omega_{2})(t-t_{0})}\left[\mathcal{L}_{e}(\omega_{1})+\mathcal{L}_{e}(\omega_{2})\right]\mathcal{L}_{f}(\omega_{1}+\omega_{2}),$
(15)

which we will focus on hereafter. The remaining, global phase factor expresses
the time dependences of $a_{1,2}(t)$ from (2) in the interaction picture of
$H_{\mathrm{f}}$; it therefore carries no physical significance – it drops out
in the calculation of the probability (7) – and can be ignored in all our
subsequent calculations.

Let us also observe that (15) can be easily adapted to a more intricate level
structure of the target, with several intermediate states, as presented in
[40]. All is needed is a $\sum_{e}$ in (6), which, by linearity, can be
carried directly throughout all derivations above, to obtain:

$T_{t}(\omega_{1},\omega_{2})=\rme^{-\rmi(\omega_{1}+\omega_{2})(t-t_{0})}\sum_{e}\left[\mathcal{L}_{e}(\omega_{1})+\mathcal{L}_{e}(\omega_{2})\right]\mathcal{L}_{f}(\omega_{1}+\omega_{2}).$
(16)

The structure of the matter response function in (15) directly reflects the
absorption process in the matter: while the term
$\mathcal{L}_{e}(\omega_{1})+\mathcal{L}_{e}(\omega_{2})$ describes the
transition $\ket{g}\rightarrow\ket{e}$ induced by either one of the two
photons, the term $\mathcal{L}_{f}(\omega_{1}+\omega_{2})$ correlates the two
photons’ frequencies by requiring that their sum be resonant with the two-
photon transition $\ket{g}\rightarrow\ket{f}$. For this reason, the analysis
presented in the next section rather depends on the detuning between the
constituent one photon transitions of the two photon process under study. To
quantify the departure from the spectrum of two non-interacting two-level
systems (where $\omega_{f}=2\omega_{e}$ and $\gamma_{f}=2\gamma_{e}$), which
can be independently driven by one single-photon pulse each, we introduce the
(dimensionless) detuning $\Delta$ and the deviation $\delta$ as

$\eqalign{\Delta=\left(\omega_{f}-2\omega_{e}\right)/\gamma_{e},\\\
\delta=\gamma_{f}/\gamma_{e}-2,}$ (17)

as indicated in figure 1(a).

###  2.5 On notation

Let us conclude this theoretical introduction by clarifying, in one single
place, the notation that we use throughout the next sections.

We deal with one-photon and two-photon states: the former are indicated with
small, the latter with capital Greek letters. Consider a one-photon state, for
a beam whose index we momentarily ignore. To refer to its Hilbert space
element we use the common ket notation: $\ket{\psi}$. The _frequency
representation_ $\psi(\omega)$, or _wave function_ , of the one-photon state
$\ket{\psi}$ is its amplitude in the mode $\omega$:

$\psi(\omega)=\braket{\omega}{\psi}\iff\ket{\psi}=\int\rmd\omega\,\psi(\omega)a^{\dagger}(\omega)\ket{0},$
(18)

where we introduced the continuous-mode single-photon state [42]

$\ket{\omega}=\ket{1_{\omega}}=a^{\dagger}(\omega)\ket{0},$ (19)

and the continuous-mode creation operator $a^{\dagger}(\omega)$ as presented
in section 2.1. For a two-photon state the notation and its interpretation are
completely analogous, except for the indices of the two fields:

$\Psi(\omega_{1},\omega_{2})=\braket{\omega_{1},\omega_{2}}{\Psi}\iff\ket{\Psi}=\iint\rmd\omega_{1}\rmd\omega_{2}\,\Psi(\omega_{1},\omega_{2})a^{\dagger}(\omega_{1})a^{\dagger}(\omega_{2})\ket{0}.$
(20)

Across section 3 and section 4 we switch between these two notations. The
frequency representation eases the physical interpretation and makes some
derivations mathematically more transparent. However, we switch to the Hilbert
space when the derivations benefit from a more compact notation. We spell out,
to exemplify the change of notation, the functional derivative [44] of the
overlap of two one-photon wave functions. In the frequency representation we
write

$\mathcal{F}[\phi,\psi^{*}]=\int\rmd\omega\,\phi(\omega)\psi^{*}(\omega)\iff\frac{\delta\mathcal{F}}{\delta\psi^{*}}=\phi(\omega),$
(21)

which is equivalent to writing, in the Hilbert space,

$\mathcal{F}[\ket{\phi},\bra{\psi}]=\braket{\psi}{\phi}\iff\frac{\delta\mathcal{F}}{\delta\bra{\psi}}=\ket{\phi}.$
(22)

##  3 Ideal pulses

###  3.1 Fluctuation-constrained optimization

Our goal is to maximize the population $p_{f}(t)$ in state $\ket{f}$ at the
time $t$, by appropriate choice of the two-photon state $\ket{\Phi}$. The
latter is determined [40] by the extrema of the functional

$J[\ket{\Phi}]=p_{f}(t)-\lambda\left(\braket{\Phi}{n_{1}n_{2}}{\Phi}-1\right),$
(23)

where the dependence of $p_{f}(t)$ on the state $\ket{\Phi}$ can be unveiled
by inserting (9) into (7), to obtain

$p_{f}(t)=|\braket{0}{T_{fg}(t)}{\Phi}|^{2}=\braket{\Phi}{T^{\dagger}_{fg}(t)}{0}\braket{0}{T_{fg}(t)}{\Phi}.$
(24)

The second term in (23) constrains the distribution of the injected two
photons over the incoming fields via the Lagrange multiplier $\lambda$: the
expectation value $\langle n_{1}n_{2}\rangle$ limits the search to two-photon
states where each beam is populated by one photon222The expectation value is
non-vanishing only for this configuration. Furthermore, (20) implies
$\braket{\Phi}{n_{1}n_{2}}{\Phi}=\braket{\Phi}{\Phi}$, such that the Lagrange
multiplier also enforces the normalization of the state.. Rather than by
saturating the number of photons in either beam [19], then, we allow the
maximization of $p_{f}(t)$ via quantum correlations [45], the potential of
which we want to scrutinize here.

At an extremum of $J[\ket{\Phi}]$, its functional derivative, (21,22), must
vanish,

$\frac{\delta J}{\delta\bra{\Phi}}=0.$ (25)

With (24) in (23), this requirement results in the eigenvalue problem

$T^{\dagger}_{fg}(t)\ket{0}\braket{0}{T_{fg}(t)}{\Phi}=\lambda\ket{\Phi},$
(26)

which, with the definition

$\ket{T}=T^{\dagger}_{fg}(t)\ket{0},$ (27)

turns into

$\ket{T}\braket{T}{\Phi}=\lambda\ket{\Phi}.$ (28)

In addition, requiring the variation of the functional (23) with respect to
the Lagrange multiplier to vanish, i.e. $\delta J/\delta\lambda=0$, enforces
the normalization of the two-photon state, such that we arrive at

$\ket{\Phi}=\mathcal{N}^{-1/2}T^{\dagger}_{fg}(t)\ket{0}$ (29)

as the only solution333Since $\ket{T}\bra{T}$ is a one-dimensional projector,
it has only one non-trivial eigenstate.. The maximal population can thus be
expressed directly in terms of $\mathcal{N}$, by (29) in (24), together with
$\braket{\Phi}{\Phi}=1$ as imposed by the constraint in (23), or in frequency
space, with (12,13,15) in (29) in (24), to obtain an explicit expression for
$\mathcal{N}$ in terms of the system parameters [40]:

$p_{f}(t)=\mathcal{N}=\iint|T_{t}(\omega_{1},\omega_{2})|^{2}\,\rmd\omega_{1}\,\rmd\omega_{2}=\frac{2\pi^{2}(\mu_{ge}\mathcal{E}_{0})^{2}(\mu_{ef}\mathcal{E}_{0})^{2}}{\hbar^{4}\gamma_{e}\gamma_{f}}.$
(30)

We remark here that the maximal population does _not_ depend on time, yet it
is attained at $t$ given the initial time $t_{0}$. This is due to the time
dependence carried by the matter response function (15) in the interaction
picture. As long as $t-t_{0}\gg\gamma_{e}^{-1}$, the optimization problem will
yield the optimal population (30) for arbitrary choice of $t_{0}$ and $t$.

###  3.2 Optimal two-photon states

To discuss the properties of the optimal two-photon state, we move to the
frequency representation, following the prescriptions illustrated in section
2.5. With (12) in (29), the optimal two-photon wave function (20) is given by

$\Phi(\omega_{1},\omega_{2})=\braket{\omega_{1},\omega_{2}}{\Phi}=\frac{T^{*}_{t}(\omega_{1},\omega_{2})}{\sqrt{\mathcal{N}}},$
(31)

which determines all the statistical properties of the two injected photons.
Up to normalization, the optimal two-photon wave function in frequency space
is the complex conjugate of the matter response function (15). This means that
the state which maximizes the population in $\ket{f}$ at time $t$ is given by
the time-reversed two-photon state emitted by the three-level system initially
(at $t$) prepared in $\ket{f}$. This is the direct two-photon analogue of the
well-known result of the optimal, “exponentially rising” single-photon state
to excite a two-level atom [46, 47], which is simply the time reversed version
of the single photon wavepacket emitted by an excited two-level atom.

Figure 1(b) shows the properties of the optimal two-photon state as defined by
(31), and parametrised by detuning $\Delta$ and deviation $\delta$. The
concentration of
$|T_{t}(\omega_{1},\omega_{2})|^{2}\propto|\Phi(\omega_{1},\omega_{2})|^{2}$
around the anti-diagonal $\omega_{1}+\omega_{2}=\omega_{f}$ reflects the
constraint, introduced by $\mathcal{L}_{f}(\omega_{1}+\omega_{2})$ in (15),
and discussed in 2.4, that the total frequency of the injected photons be
resonant with the two-photon transition $\ket{g}\rightarrow\ket{f}$. The
frequency sum distribution is readily extracted from
$|\Phi(\omega_{1},\omega_{2})|^{2}$, by changing variables to
$\omega_{\pm}=\omega_{1}\pm\omega_{2}$ and integrating over $\omega_{-}$:

$p_{\rm
sum}(\omega_{+})=\frac{1}{\pi}\frac{\gamma_{f}}{(\omega_{+}-\omega_{f})^{2}+\gamma_{f}^{2}},$
(32)

i.e. a Lorentzian of width $\gamma_{f}=(2+\delta)\gamma_{e}$, centred at
$\omega_{f}=2\omega_{e}+\Delta\gamma_{e}$.

In figure 1(c), instead, we plot the single-photon distribution, i.e. the
marginal distribution of $|\Phi(\omega_{1},\omega_{2})|^{2}$ with respect to
either frequency (here the first), since the optimal two-photon wave function
(31) is symmetric:

$p_{1}(\omega)=\frac{\gamma_{e}(\gamma_{e}+\gamma_{f})(4\gamma_{e}+\gamma_{f})+\gamma_{e}(\omega_{f}-2\omega_{e})^{2}+\gamma_{f}(\omega-\omega_{e})^{2}}{2\pi\left[(\omega-\omega_{e})^{2}+\gamma_{e}^{2}\right]\left[(\omega-\omega_{f}+\omega_{e})^{2}+(\gamma_{e}+\gamma_{f})^{2}\right]},$
(33)

where we can identify two distinct peaks near444The peaks do not exactly match
the transition frequencies because
$|\mathcal{L}_{e}(\omega_{1})+\mathcal{L}_{e}(\omega_{2})|^{2}\neq|\mathcal{L}_{e}(\omega_{1})|^{2}+|\mathcal{L}_{e}(\omega_{2})|^{2}$,
and interference terms contribute to $p_{1}(\omega)$. the transitions
$\ket{g}\rightarrow\ket{e}$ (frequency $\omega_{e}$ and line width
$\gamma_{e}$) and $\ket{e}\rightarrow\ket{f}$ (frequency
$\omega_{f}-\omega_{e}=\omega_{e}+\Delta\gamma_{e}$ and line width
$\gamma_{f}+\gamma_{e}=(3+\delta)\gamma_{e}$). Notice that, since
$\Phi_{t}(\omega_{1},\omega_{2})$ is symmetric under exchange of variables,
the single-photon frequency distribution is the same for both fields: each
photon has the same probability to induce the transition to $\ket{e}$, and the
second completes it to $\ket{f}$.

When $\delta\rightarrow-2$, that is in the limit of a final state with
vanishing linewidth $\gamma_{f}\rightarrow 0$, $p_{\rm sum}(\omega_{+})$ tends
to $\delta(\omega_{+}-\omega_{f})$. This describes the situation where the
frequency of one “emitted” photon strictly determines that of the other, which
means that each photon has the same probability of exciting either the
$\ket{g}\rightarrow\ket{e}$ or the $\ket{e}\rightarrow\ket{f}$ transition. In
this limit, $p_{1}(\omega)$ exhibits two peaks of equal (unit) width at
$\omega_{e}$ and $\omega_{f}-\omega_{e}$, as discussed above.

###  3.3 Entanglement entropy

When $\Delta=\delta=0$, i.e., by (17), $\omega_{f}=2\omega_{e}$ and
$\gamma_{f}=2\gamma_{e}$, the optimal two-photon wave function given by (31)
is separable: either photon independently excites one of two non-interacting
two-level systems, see section 2.4. Otherwise, the frequency degrees of
freedom of the two optimally prepared incoming pulses are, in general,
entangled [40], with non-trivial (i.e., $L>1$) Schmidt decomposition into
Schmidt modes $\varphi^{\ast}_{k}(\omega_{1}),\psi^{\ast}_{k}(\omega_{2})$ and
their non-increasingly ordered, non-negative Schmidt coefficients $r_{k}$,

$\displaystyle\Phi(\omega_{1},\omega_{2})$
$\displaystyle=\sum_{k=1}^{L}r_{k}\varphi^{\ast}_{k}(\omega_{1})\psi^{\ast}_{k}(\omega_{2}),$
(34)

where normalization via $\mathcal{N}$ in (31) implies $\sum_{k}r_{k}^{2}=1$
($0\leq r_{k}\leq 1$), and complex conjugation is inherited from that of
$T_{t}(\omega_{1},\omega_{2})$, also in (31). The entanglement encoded in
$\Phi(\omega_{1},\omega_{2})$ can be quantified by the entanglement entropy
[45],

$S=-\sum_{k=1}^{L}r_{k}^{2}\log_{2}r_{k}^{2}.$ (35)

We now inspect how the degree of frequency entanglement of the incoming pulses
correlates with enhanced absorption probabilities.

###  3.4 Quantum enhancement

By (30), the maximal population of the state $\ket{f}$ is $\mathcal{N}$. To
capture the genuine advantage due to entanglement, an entangled fields’ yield
is to be compared with that of the optimal _separable_ two-photon state

$\Phi_{\mathrm{sep}}(\omega_{1},\omega_{2})=\varphi_{1}^{*}(\omega_{1})\psi_{1}^{*}(\omega_{2}),$
(36)

given by the modes pertaining to the largest Schmidt coefficient $r_{1}$ in
(35) [40]. Such an optimal classical state yields an excited state population
$r_{1}^{2}\mathcal{N}\leq\mathcal{N}$, where equality is achieved only for
$\Delta=\delta=0$. The _quantum enhancement_ of TPA achievable by an entangled
two-photon state, through the quantum correlations between the incoming
fields, is thus given by the ratio

$\displaystyle E_{\mathrm{q}}=r_{1}^{-2}\geq 1$ (37)

of those two excitation probabilities.

According to (35), strong entanglement $S$ requires $L\gg 1$ and thus small
$r_{1}$, which implies strong quantum enhancement $E_{\mathrm{q}}$, by (37).
We test this mutual relationship in figure 2, where we evaluate (36) and (37)
with the optimal two-photon wave function (31), for variable deviation
$\delta$ (panel (a)) and detuning $\Delta$ (panel (b)), respectively, as well
as the maximally achievable values of entanglement and enhancement,
$S_{\infty}$ and $E_{\infty}$, respectively, in the limit of very large
$\Delta$ (panel (c)), which we discuss separately in section 3.5. 555All
results here displayed are based on the Schmidt decomposition (34) of the two-
photon wave function $\Phi(\omega_{1},\omega_{2})$ discretized in frequency
space, i.e., we used a linear algebra package (Wolfram Mathematica) to
calculate the singular values of a matrix. Frequencies were discretized on a
grid of size $\pm 200\gamma_{f}$ and resolution $\gamma_{e}/5$. This choice of
grid size and resolution was validated by a normalization of the discretized
state which was systematically bounded from below by $0.99$.

Figure 2: Entanglement entropy $S$ as defined in (35), and quantum enhancement
$E_{\mathrm{q}}$ given in (37), of the excited state population $p_{f}(t)$
(24) achieved by the optimal two photon state (31) for (a) fixed detuning
$\Delta=5$ with variable deviation $\delta\in(-2;2]$ (note the colour code,
and specific values indicated by labels), and (b) fixed deviation
$\delta=-1.9$ with variable detuning $\Delta\in[0.1;100]$ (colour coded and
labeled). The dashed lines indicate the maximally achievable entanglement and
quantum enhancement, respectively, in the limit $\Delta\rightarrow\infty$. (c)
Maximally achievable enhancement $E_{\infty}$, see (42) and Sec. 3.5, and
associated entanglement entropy $S_{\infty}$, given by (43), as a function of
$\delta$, in the limit $\Delta\rightarrow\infty$ (achieved by
$\gamma_{e}\rightarrow 0$).

In panels (a,b), we observe a monotonic increase of $E_{\mathrm{q}}$ with $S$,
and a plateau of the achievable entanglement and quantum enhancement emerges
for large values of the detuning $\Delta$, which we address in the next
subsection. The increase of $S$ with $\Delta$ and as $\delta\rightarrow-2$ can
be attributed to a narrowing of the frequency-sum distribution (32), in unison
with a broad single frequency distribution (33) (due to its composition by two
distant peaks). This distribution signifies strong frequency anti-
correlations, and indeed, since we are here considering pure quantum states of
light, correspond to a strongly entangled wavefunction [19].

###  3.5 Maximal quantum enhancement

Figure 3: First two pairs $\varphi_{k}(\omega)$, $\psi_{k}(\omega)$ of Schmidt
modes computed for a response function (15) with $\Delta=100$ and
$\delta=-1.5$, discretized on a grid $[-500;500]$ with a resolution
$(2+\delta)/2$. These pairs correspond to the Schmidt coefficients
$r_{1}=0.272557$ and $r_{2}=0.272348$. Because of the large value of $\Delta$,
the Schmidt modes are linear combinations of modes $\alpha_{k}$ and
$\beta_{k}$ centred, respectively, at rescaled frequencies
$(\omega-\omega_{e})/\gamma_{e}$ $0$ and $\Delta$.

Let us now discuss the quantum enhancement $E_{\mathrm{q}}$ in the limit
$\Delta\gg 1$ (or, more physically, $\omega_{f}-2\omega_{e}\gg\gamma_{e}$),
which figure 2(b) suggests to be finite. As discussed in section 2.4, the
matter response function (15) is symmetric in the two photon frequencies,
which implies that its Schmidt modes $\varphi_{k}$, $\psi_{k}$ in (34) are
equal for each $k$, possibly up to a phase. For $\Delta\gg 1$, as in figure 3,
said modes appear as orthogonal superpositions of non-overlapping, complex-
valued line shapes $\alpha_{k}$ (centred at $\omega_{e}$) and $\beta_{k}$
(centred at $\omega_{f}-\omega_{e}=\omega_{e}+\Delta\gamma_{e}$). Moreover,
mutually orthogonal linear combinations of the same modes, e.g.
$\varphi_{1,2}=(\alpha_{1}\pm\beta_{1})/\sqrt{2}$ in figure 3, are associated
with Schmidt coefficients $r_{1}\approx r_{2}$, at least for finite $\Delta$.

Let us now understand the connection between the Schmidt modes $\varphi_{k}$,
$\psi_{k}$ and the line shapes $\alpha_{k}$, $\beta_{k}$. The latter can be
considered the Schmidt modes of a response function where the first photon is
resonant with the $\ket{g}\rightarrow\ket{e}$ transition, while the second
completes the two-photon transition by inducing $\ket{e}\rightarrow\ket{f}$.
The construction of such a response function yields

$Q_{t}(\omega_{1},\omega_{2})=\rme^{-\rmi(\omega_{1}+\omega_{2})(t-t_{0})}\mathcal{L}_{e}(\omega_{1})\mathcal{L}_{f}(\omega_{1}+\omega_{2}).$
(38)

The analogy to (15) becomes clear by noticing that

$T_{t}(\omega_{1},\omega_{2})=Q_{t}(\omega_{1},\omega_{2})+Q_{t}(\omega_{2},\omega_{1}),$
(39)

that is, $T_{t}$ is the symmetrized666It hence allows each photon to begin the
two-photon transition, and the other to complete it. version of $Q_{t}$. For
visualization, $Q_{t}(\omega_{1},\omega_{2})$ and
$Q_{t}(\omega_{2},\omega_{1})$ describe, respectively, the top-left and
bottom-right peaks of $T_{t}(\omega_{1},\omega_{2})$ in figure 1(b).

Let us now write the Schmidt decomposition of the asymmetric response function
(38) as

$Q_{t}(\omega_{1},\omega_{2})=\sqrt{\frac{\mathcal{N}}{2}}\sum_{k=1}^{L^{\prime}}s_{k}\alpha_{k}(\omega_{1})\beta_{k}(\omega_{2}),$
(40)

where the prefactor ensures $\sum_{k}s_{k}^{2}=1$. The Schmidt coefficients
must be independent of $\Delta$, since this parameter translates the top-left
peak in figure 1(b) without changing its anti-diagonal structure. This cannot
be true for $T_{t}$ in (15), where a change in $\Delta$ implies both a
vertical and a horizontal translation of, respectively, the top-left and
bottom-right peaks. In this case, then, the structure on the anti-diagonal
changes and so do the correlations between the two photon frequencies.

Substituting (34) (on the left-hand side, via (31)) and (40) (on the right-
hand side) in the identity (39) we obtain

$\frac{T_{t}(\omega_{1},\omega_{2})}{\sqrt{\mathcal{N}}}=\sum_{k=1}^{L}r_{k}\varphi_{k}(\omega_{1})\psi_{k}(\omega_{1})=\sum_{l=1}^{L^{\prime}}\frac{s_{l}}{\sqrt{2}}\left[\alpha_{l}(\omega_{1})\beta_{l}(\omega_{2})+\beta_{l}(\omega_{1})\alpha_{l}(\omega_{2})\right].$

Going back to the example of figure 3, if we assume $r_{1}=r_{2}$ and multiply
out the combinations proposed for $\varphi_{k}(\omega_{1})$ and
$\psi_{k}(\omega_{2})$, for $k=1,2$, we obtain precisely the products of
$\alpha_{1}$ and $\beta_{1}$ on the right-hand side of (3.5). Notice that
these combinations are not arbitrary: as discussed earlier, because $T_{t}$ is
a symmetric complex function, the modes $\varphi_{k}$ and $\psi_{k}$ might
differ by a phase that ensures the positivity of the singular values. All in
all, we must satisfy
$\varphi_{1}(\omega_{1})\psi_{1}(\omega_{2})+\varphi_{2}(\omega_{1})\psi_{2}(\omega_{2})=\alpha_{1}(\omega_{1})\beta_{1}(\omega_{2})+\beta_{1}(\omega_{1})\alpha_{1}(\omega_{2})$.

The right-hand side of (3.5) is a valid Schmidt decomposition _only_ in the
limit $\Delta\rightarrow\infty$, where the function bases
$\\{\alpha_{k}\\}_{k}$ and $\\{\beta_{k}\\}_{k}$ are also mutually
orthogonal777The overlap of these line shapes goes to zero in the limit
$\Delta\rightarrow\infty$., since – to reconstruct the symmetry of $T_{t}$ –
they both have to appear as Schmidt modes for each photon frequency. Under
this condition, then, the Schmidt coefficients $r_{k}$ of the response
function $T_{t}$ must come in pairs. In summary, in the limit
$\Delta\rightarrow\infty$, eq. (3.5) is a valid Schmidt decomposition and we
can identify

$r_{2k}=r_{2k-1}=\frac{s_{k}}{\sqrt{2}},\quad\varphi_{2k-1}=\psi_{2k}=\alpha_{k},\quad\varphi_{2k}=\psi_{2k-1}=\beta_{k},\qquad
k=1,\ldots,L^{\prime}.$ (41)

The pairwise appearance of Schmidt coefficients $r_{1}=s_{1}/\sqrt{2}$, by
(41), implies that the enhancement $E_{\infty}$, when $\Delta\gg 1$, is twice
the enhancement $E_{\mathrm{a}}$ obtainable with the asymmetric matter
response function (38):

$E_{\infty}=2E_{\mathrm{a}}.$ (42)

For the entropy, instead, we can write

$S_{\infty}=-\sum_{k=1}^{L}r_{k}^{2}\log_{2}r_{k}^{2}=-2\sum_{l=1}^{L^{\prime}}\frac{s_{l}^{2}}{2}\log_{2}\left(\frac{s_{l}^{2}}{2}\right)=1+S_{\mathrm{a}},$
(43)

where $S_{\mathrm{a}}$ is the entanglement entropy of (40). For any value of
$\delta$, then, the quantum enhancement induced by the optimal pulse (31), and
the entanglement between the two photons’ frequencies must be bounded,
respectively, by $E_{\infty}$ and $S_{\infty}$, which are plotted in figure 2
(c). For $\delta\rightarrow-2$ we observe the same steep increase in
enhancement due to the “strict” correlation between the photon frequencies
discussed in sections 3.2 and 3.4.

##  4 Realistic pulses

Figure 4:  Exemplary experimental scenarios discussed in (a) section 4.3 and
(b) section 4.4. In (a) two photons of frequencies $\omega_{1}$ and
$\omega_{2}$ are produced via spontaneous parametric down-conversion (SPDC)
from a pump photon of frequency $\omega_{p}$. Both photons are then frequency-
modulated by the same pulse shaping operator $M$. In (b), instead, it is the
pump photon that is shaped by an optimal operator $M$, before being split into
two photons by SPDC. (c) recalls, for convenience, the parameters describing
the matter degrees of freedom as in figure 1(a).

So far, our investigation targeted the optimal two-photon quantum state which
can be constructed theoretically. However, the experimental manipulation of
the two-photon state is subject to further constraints, and it is therefore
necessary to reformulate the above theory for experimentally implementable
transformations. This is the purpose of the present section. We consider an
experiment in which a two-photon state $\ket{\Sigma}$ can be generated. The
frequency components of this state will then be manipulated, e.g. using a
spatial light modulator [35, 36, 37], to enhance the propensity of the state
to excite a two-photon transition in the three-level sample it impinges on, as
in figure 4(c). First we analyze a general class of transformations, and then
consider two realistic examples (see figures 4(a) and (b)), where the two-
photon state generated e.g. by spontaneous parametric down-conversion is
spectrally shaped (a), or where the (classical) pump pulse that creates the
entangled photon pair by down-conversion is transformed (b). As we will see,
both situations can be analyzed in terms of unitary transformations acting on
the two-photon state $\ket{\Sigma}$.

###  4.1 Optimization procedure

We want to find unitary operations, $M_{1}$ and $M_{2}$, which act on the
Hilbert spaces of the first and of the second photon, respectively, and
maximize the TPA probability induced by the manipulated state
$M_{1}M_{2}\ket{\Sigma}$. If we continue to denote with $\ket{\Phi}$ the
optimal state (29) analyzed in section 3, we intuitively want
$M_{1}M_{2}\ket{\Sigma}$ to most closely resemble $\ket{\Phi}$.

In the frequency representation of section 2.5, we write a pulse shaping
_function_ as

$\braket{\omega_{j}}{{M}_{j}}{\nu_{j}}=M_{j}(\nu_{j},\omega_{j}).$ (44)

In section 4.3, the arguments in the above equation refers to the frequency
components of the photon wavepacket in beam $j$. As we will see in section
4.4, the same formalism can be applied to the shaping of the pump pulse in
type-I parametric down-conversion, where, however, the frequency arguments
represent the components of the sum frequency $\omega_{1}+\omega_{2}$.

While in sections 4.2–4.4 we assume $M$ to be diagonal in frequency space,
here we include operations that change the spectral shape of the photon
wavepacket, as described, e.g., in section V of [48]. We do, however, assume
that the photon number is conserved, i.e. squeezing or displacement operations
are not considered. As in section 3, we want to maximize the final state
population (7), but with the shaped two-photon state $M_{1}M_{2}\ket{\Sigma}$:

$p_{f}(t)=|\braket{0}{{T}_{fg}(t)M_{1}M_{2}}{\Sigma}|^{2}=\mathcal{N}|\braket{\Phi}{M_{1}M_{2}}{\Sigma}|^{2}.$
(45)

The factor $\mathcal{N}$ stems from (29), and clarifies the meaning of this
equation: if $\mathcal{N}$ is the maximal population achievable by using the
optimal state $\ket{\Phi}$, then the population in the case of a shaped
realistic state $M_{1}M_{2}\ket{\Sigma}$ is reduced by the overlap between
these two initial states.

The new functional to optimize – instead of (23) – is then

$J[M_{1},M_{2}]=p_{f}(t)-\sum_{i=1}^{2}\left[\braket{\psi_{j}}{M_{j}^{\dagger}M_{j}}{\psi_{j}}-\braket{\psi_{j}}{\psi_{j}}\right].$
(46)

We remark here that (23) was maximized over the space of two-photon states,
whereas now we are optimizing in the space of _operators_ on the Hilbert
spaces of the individual photons. In (46) we constrain the search to local
operators that do not change the number of photons of undetermined, and
potentially unnormalized, single-photon states $\ket{\psi_{j}}$ as per (18).
The summation term in (46) therefore limits our search to unitary local
operators:

$\displaystyle M_{j}^{\dagger}M_{j}=\mathbb{I}.$ (47)

The single-photon states $\ket{\psi_{j}}$ do not carry further relevance
beyond imposing the constraint (47), and effectively play the same role of the
Lagrange multiplier of (23).

Due to the higher dimensionality of the search space, we cannot find an
explicit _general_ expression of the optimal pulse shaping operators $M_{j}$.
These are solutions of coupled integral equations (86), whose derivation we
defer to A.

###  4.2 Diagonal pulse shaping operators

We now collect the theoretical underpinning of our subsequent discussion of
examples in sections 4.3 and 4.4. We focus on _diagonal_ pulse shaping
operators, i.e. of the form

$\braket{\omega_{j}}{M_{j}}{\nu_{j}}=M_{j}(\omega_{j})\delta(\omega_{j}-\nu_{j}).$
(48)

This means that we restrict ourselves to linear optical elements that do not
change the photon energies, such as the aforementioned manipulation by a
spatial light modulator. Substituting (48) in the derivation of A, we obtain,
in place of (86),

$M_{j}(\omega_{j})=\frac{\sqrt{\mathcal{N}}\mathcal{A}[M_{1},M_{2}]}{|\psi_{j}(\omega_{j})|^{2}}\int\mathcal{W}_{t}^{*}(\omega_{j},\omega_{k})M_{k}^{*}(\omega_{k})\,\rmd\omega_{k}.$
(49)

The functional $\mathcal{A}[M_{1},M_{2}]$ is defined in (80), and we have also
introduced the _effective_ response function $\mathcal{W}_{t}$ as the product
of the matter response function (15, 31) with the (frequency representation of
the) input state $\ket{\Sigma}$:

$\mathcal{W}_{t}(\omega_{1},\omega_{2})=\Sigma(\omega_{1},\omega_{2})T_{t}(\omega_{1},\omega_{2})=\sqrt{\mathcal{N}}\Sigma(\omega_{1},\omega_{2})\Phi^{*}(\omega_{1},\omega_{2}).$
(50)

Using the effective response function, the final state population (45)
achievable via the diagonal pulse shaping functions $M_{1}$ and $M_{2}$ of
(48) reads

$p_{f}(t)=\left|\iint\mathcal{W}_{t}(\omega_{1},\omega_{2})M_{1}(\omega_{1})M_{2}(\omega_{2})\,\rmd\omega_{1}\,\rmd\omega_{2}\right|^{2}.$
(51)

Given the fact that $M_{j}$ is a unitary diagonal operator, see (47), we can
write

$M_{j}(\omega_{j})=\rme^{-\rmi F_{j}(\omega_{j})},$ (52)

with $F_{j}(\omega_{j})$ a real function. Substitution of this ansatz in (49)
yields

$\displaystyle|\psi_{j}(\omega_{j})|^{2}\rme^{-\rmi F_{j}(\omega_{j})}=$
$\displaystyle\sqrt{\mathcal{N}}\iint\mathcal{W}_{t}(\omega_{1},\omega_{2})\rme^{-\rmi\left[F_{1}(\omega_{1})+F_{2}(\omega_{2})\right]}\rmd\omega_{1}\rmd\omega_{2}$
(53)
$\displaystyle\times\int\mathcal{W}_{t}^{*}(\omega_{j},\omega_{k})\rme^{\rmi
F_{k}(\omega_{k})}\rmd\omega_{k}.$

If we also cast the effective response function in polar form,

$\mathcal{W}_{t}(\omega_{1},\omega_{2})=|\mathcal{W}_{t}(\omega_{1},\omega_{2})|\rme^{\rmi\mathcal{S}(\omega_{1},\omega_{2})},$
(54)

the identity (53) is satisfied by functions $F_{1}$ and $F_{2}$ such that

$F_{1}(\omega_{1})+F_{2}(\omega_{2})=\mathcal{S}(\omega_{1},\omega_{2}),$ (55)

and by a Lagrange multiplier $\psi_{j}$ with

$|\psi_{j}(\omega_{j})|^{2}=\sqrt{\mathcal{N}}\iint|\mathcal{W}_{t}(\omega_{1},\omega_{2})|\rmd\omega_{1}\rmd\omega_{2}\int|\mathcal{W}_{t}(\omega_{j},\omega_{k})|\rmd\omega_{k}.$
(56)

Substituting (52) and (54) in (51), and assuming (55) is satisfied, we find
the final state population achievable upon manipulating (“s” for “shaped”) the
input state $\ket{\Sigma}$:

$p_{f}^{\mathrm{(s)}}(t)=\left(\iint|\mathcal{W}_{t}(\omega_{1},\omega_{2})|\,\rmd\omega_{1}\,\rmd\omega_{2}\right)^{2}.$
(57)

On the other hand, in the absence of shaping operators, that is with
$M_{1}=M_{2}\equiv 1$, (51) yields

$p_{f}^{\mathrm{(u)}}(t)=\left|\iint\mathcal{W}_{t}(\omega_{1},\omega_{2})\rmd\omega_{1}\rmd\omega_{2}\right|^{2}.$
(58)

This quantity is the “unshaped” or “unoptimized” population achievable when we
do not manipulate the incoming realistic pulse. Its value is determined by the
sign taken by the effective response function $\mathcal{W}_{t}$, i.e. by the
phase difference between the optimal $\ket{\Phi}$ and the realistic state
$\ket{\Sigma}$. These phase differences are precisely what the
optimal888Applying the Hölder inequality [49] to (51), and using the unitarity
of the shaper functions (47), we can write the inequality
$p_{f}(t)=\left|\iint\mathcal{W}_{t}(\omega_{1},\omega_{2})M_{1}(\omega_{1})M_{2}(\omega_{2})\,\rmd\omega_{1}\,\rmd\omega_{2}\right|^{2}\leq
p_{f}^{\mathrm{(s)}}(t)$ (59) Because equality holds for the shaper functions
satisfying (55), these must define the optimal solution. pulse shaping
functions compensate, by transforming $\mathcal{W}_{t}(\omega_{1},\omega_{2})$
into $|\mathcal{W}_{t}(\omega_{1},\omega_{2})|$ via (52) with (55). To capture
the enhancement obtained by pulse shaping we therefore introduce the
_optimization ratio_ as the ratio between optimized and unoptimized
populations:

$E_{\mathrm{opt}}=\frac{p_{f}^{\mathrm{(s)}}(t)}{p_{f}^{\mathrm{(u)}}(t)}.$
(60)

To conclude, in the following examples we will identify the effective response
function $\mathcal{W}_{t}$. If we can determine, via its phase, the argument
$F_{j}(\omega_{j})$ of the pulse shaping operators, such that (55) is
satisfied, then the optimal final state population is directly given by (57).

###  4.3 Example I: shaping down-converted photons

The first example we consider is depicted in figure 4(a): a pump photon of
frequency $\omega_{p}$ is split via SPDC into two photons of frequencies
$\omega_{1}$ and $\omega_{2}$, which are modulated by the same pulse shaping
operator $M$. This setup describes the experiment carried out by the
Silberberg group in 2005 [25]. We consider a frequency representation (20) of
$\ket{\Sigma}$ given by

$\Sigma(\omega_{1},\omega_{2})=\delta(\omega_{1}+\omega_{2}-\omega_{p})G(\omega_{1}-\omega_{p}/2),$
(61)

where $\omega_{p}$ is the frequency of the pump photon creating the entangled
pair, which we assume to be on resonance with the two-photon transition
($\omega_{p}=\omega_{f}$) in the matter, see figure 4(c). The delta function
stems from a narrowband continuous-wave (cw) pump laser and replaces the
Lorentzian distribution of the sum frequencies we encountered in (32). The cw
nature of the pump pulse implies that the arrival time of the entangled photon
pair is completely undetermined. Consequently, a targeted excitation at a
particular time $t$, as considered in the previous section, is impossible.
Rather, within this approximation we describe a steady state experiment where
a constant stream of entangled photon pairs gives rise to a finite $f$-state
population [50]. Our goal is to optimise this steady state population by
shaping the two-photon state (61).

We consider frequency-degenerate down-converted photons, i.e. the individual
photon wave packets are described by a real function $G$ symmetric around
$\omega_{f}/2$ [51]. The two photons are therefore equally detuned from either
transitions $\ket{g}\rightarrow\ket{e}$ and $\ket{e}\rightarrow\ket{f}$, since
$\omega_{f}/2=\omega_{e}+(\Delta/2)\gamma_{e}=(\omega_{f}-\omega_{e})-(\Delta/2)\gamma_{e}$.
For compactness of notation, then, we shift the frequencies as
$\Omega_{j}=\omega_{j}-\omega_{f}/2$. The effective response function (50) for
this example then reads

$\displaystyle\mathcal{W}_{t}(\Omega_{1},\Omega_{2})$
$\displaystyle=\delta(\Omega_{1}+\Omega_{2})G(\Omega_{1})T_{t}(\Omega_{1}+\omega_{f}/2,\Omega_{2}+\omega_{f}/2)$
(62)
$\displaystyle\propto\delta(\Omega_{1}+\Omega_{2})G(\Omega_{1})\left[\frac{1}{\Omega_{1}+\left(\frac{\Delta}{2}+\rmi\right)\gamma_{e}}+\frac{1}{\Omega_{2}+\left(\frac{\Delta}{2}+\rmi\right)\gamma_{e}}\right].$

In this expression we dropped all prefactors, because the optimal pulse
shaping operator $M$ is determined by the phase of
$\mathcal{W}_{t}(\Omega_{1},\Omega_{2})$, as defined in (54,55), and therefore
only999The deviation $\delta$, which determines the width of the frequency-sum
distribution (32), is implicitly contained in the normalization of
$\ket{\Sigma}$, but for sharply correlated frequencies via the Dirac delta in
(61) it does not affect the optimization. depends on the detuning $\Delta$
(see (17)). It is worth anticipating here that the time dependence of $T_{t}$
from (15) becomes a trivial phase factor $\rme^{-\rmi\omega_{f}(t-t_{0})}$,
once the Dirac delta is integrated over, as we do later; for this reason, we
do not spell it out explicitly in (62).

In the shifted frequencies $\Omega_{j}$ above, and for identical SLMs, i.e.
$M_{1}=M_{2}=M$, Eq. (49) becomes

$|\psi(\Omega)|^{2}M(\Omega)=\sqrt{\mathcal{N}}\mathcal{W}^{*}_{t}(\Omega,-\Omega)M^{*}(-\Omega)\int\mathcal{W}_{t}(\Omega^{\prime},-\Omega^{\prime})M(\Omega^{\prime})M(-\Omega^{\prime})\,\rmd\Omega^{\prime},$
(63)

where we integrated over the variables involved in the Dirac deltas of (62).
Mirroring equations (54)–(56), then, the optimal pulse shaping function reads

$M(\Omega)=\rme^{-\rmi\mathcal{S}(\Omega,-\Omega)/2},$ (64)

with $\mathcal{S}$ given in (54), while

$|\psi(\Omega)|^{2}=\sqrt{\mathcal{N}}|\mathcal{W}_{t}(\Omega,-\Omega)|\int|{\mathcal{W}}_{t}(\Omega,-\Omega)|\,\rmd\Omega.$
(65)

These expressions presuppose that the effective response function (62) is
symmetric in $\Omega$, and that $G(\Omega)$ is real, as assumed above.

In this example, the optimization ratio (60) reads

$E_{\mathrm{opt}}=\frac{\left(\int|\mathcal{W}_{t}(\Omega,-\Omega)|\,\rmd\Omega\right)^{2}}{\left|\int\mathcal{W}_{t}(\Omega,-\Omega)\,\rmd\Omega\right|^{2}},$
(66)

and we plot its values in figure 5(a) for the individual photons described by
the Gaussian profile

$G(\Omega)=\frac{1}{(\pi\sigma^{2})^{1/4}}\rme^{-\Omega^{2}/(2\sigma^{2})},$
(67)

where $\sigma$ denotes the bandwidth of the photon wave packet.

Let us remind the reader of the physical quantities we are comparing here. For
the matter degrees of freedom, we want to drive the two transitions
$\ket{g}\rightarrow\ket{e}$ and $\ket{e}\rightarrow\ket{f}$, of frequencies,
respectively, $\omega_{e}$ and $\omega_{e}+\Delta\gamma_{e}$, see figure 4(c).
We excite the transitions with two photons with frequency profiles centred on
$\omega_{f}/2=\omega_{e}+(\Delta/2)\gamma_{e}$ and width $\sigma$. As
discussed above, the two photons are always detuned by $(\Delta/2)\gamma_{e}$
from the electronic transitions. The probability of exciting the individual
transitions depends, therefore, on the width $\sigma$ of the single-photon
frequency distributions.

If $\Delta=0$, i.e. in the case of no detuning, each photon can resonantly
drive either transition $\ket{g}\rightarrow\ket{e}$ or
$\ket{e}\rightarrow\ket{f}$. In this case, then, there is nothing to optimize.
This can be shown formally with little algebraic manipulation:
$\mathcal{W}_{t}(\Omega,-\Omega)$ from (62) becomes a real function, and hence
carries no phase to compensate via (55). When $\sigma\ll\Delta$, instead, the
photon frequencies are so narrowly distributed around $\omega_{f}/2$ that they
are always off resonant (by $\Delta/2$) with respect to the two electronic
transitions; shaping the incoming entangled wave function barely affects the
excitation probability in this case. Between these two extreme cases, for
fixed detuning $\Delta$, $E_{\mathrm{opt}}$ reaches a maximum and then
saturates for $\sigma\gtrsim\Delta$, when the individual photon pulses are so
broad to be resonant with any transition in the matter. In figure 5(b) we thus
consider the case where we can freely tune the bandwidth $\sigma$ of the
photon pulses to match the detuning $\Delta$, and plot the maximal final state
population $p_{f}$ achievable by both unshaped (61) and shaped (64) Gaussian
wave packets. Remarkably, the latter departs rather slowly from the maximal
value $\mathcal{N}$ (30) achievable by the optimal state (29). Hence, in this
case, shaping can almost saturate the optimal bound $E_{\mathrm{q}}$ allowed
by quantum mechanics.

Figure 5:  (a) Optimization ratio $E_{\mathrm{opt}}$ (60) as a function of
$\Delta$ (50), obtained by optimally shaping a realistic the two-photon wave
function (61), given the single-photon Gaussian profile (67), for three
different values of the latter’s width $\sigma$. (b) Maximum achievable
population, in units of $\mathcal{N}$ (30), for optimally shaped [(64), solid
line] and unshaped [(61), dashed line], realistic two-photon wave functions,
with Gaussian single-photon profiles (67) of width $\sigma=\Delta$, with
$\Delta\geq 0.1$.  Figure 6: Comparison of the quantum enhancement
$E_{\mathrm{q}}$, with respect to the optimal separable pulses (36),
achievable (a) via the optimal entangled two-photon wave function (31), (b) by
optimal shaping (64) of a realistic pulse (61) with Gaussian single-photon
distribution (67), and (c) by an unshaped realistic pulse (61), with Gaussian
profile (67). In all panels we explore the ranges $\Delta\in[0.1;5]$ and
$\delta\in[-1.9;0]$. In panels (b) and (c) the width of the Gaussian
distribution is set to $\sigma=\Delta$ for each pair $(\delta,\Delta)$, to
enable the comparison with the optimal separable pulses. The dashed line in
(b) and (c) traces the pairs $(\delta,\Delta)$ where the realistic two-photon
wave functions [shaped by (64) and unoptimized (61), respectively] yield
$E_{\mathrm{q}}=1$, i.e. it demarcates the region (above/to the left of the
dashed line) where the optimized states can outperform optimal separable
states. Panel (a) reproduces the numerical results of [40].

Finally, the remaining open question is how optimal pulse shaping (64) of the
realistic entangled state $\ket{\Sigma}$ compares to the classical limit set
by optimal separable pulses (36). As in figure 5(b), in figure 6 we match the
width $\sigma$ of the Gaussian pulse (67) to the detuning $\Delta$, to cover
both transition frequencies. As discussed for (62), the deviation $\delta$
(see (17)) does not affect $\mathcal{W}_{t}(\Omega,-\Omega)$; the dependence
on $\delta$ visible in figure 6(b) and (c) stems from the optimal population
achieved by the corresponding separable pulse. Our analysis finds large areas
of parameter space, i.e. the entire upper right part of figure 6(b), for which
only shaped Gaussian pulses can violate the classical limit (compare the areas
in (b) and (c) which lie on the right/below the dashed line – these are the
areas where an enhancement $E_{\mathrm{q}}\leq 1$ with respect to the optimal
separable pulse (36) is achieved).

###  4.4 Example II: shaping the pump

As sketched in figure 4(b), in this example we consider shaping a pump pulse
that is subsequently down-converted. Hence, in contrast to our previous
example, we will not approximate the pump laser by a spectral delta function
and consider down-conversion driven by a finite-bandwidth pulse instead. The
effect of pulse shaping on the entanglement of down-converted photon pairs was
recently studied in [52] for a spectrally chirped pump pulse. Spectral chirp
is a common effect in experiments that manifests as a quadratic phase in the
pulse shape of the pump photon. Here we discuss how this phase can be
counteracted via pulse shaping as described in section 4.2.

We consider a two-photon state of the form

$\Sigma(\omega_{1},\omega_{2})=\alpha(\omega_{1}+\omega_{2})\beta(\omega_{1}-\omega_{2}).$
(68)

This is an appropriate model, for instance, for photon pairs created by type-I
down-conversion [53, 54, 55]. Here, $\alpha(\omega)$ is proportional to the
amplitude of the pump laser driving the down-conversion, and $\beta(\omega)$
denotes the phase-matching function. To apply the optimization procedure of
section 4.2, we change variables to $\omega_{\pm}=\omega_{1}\pm\omega_{2}$ and
write the effective response function (50) as

$\mathcal{W}_{t}(\omega_{+},\omega_{-})=\frac{1}{2}\alpha(\omega_{+})\beta(\omega_{-})T_{t}\left(\frac{\omega_{+}+\omega_{-}}{2},\frac{\omega_{+}-\omega_{-}}{2}\right),$
(69)

where the factor $1/2$ ensures the correct change of variables under
integration. Since the pump frequency determines the sum of the down-converted
photons’ frequencies, we want to shape the distribution of
$\omega_{+}=\omega_{1}+\omega_{2}$, which, in contrast to the previous
example, is not a LOCC (local operation with classical communication, [56])
and can change the amount of entanglement in the final two-photon state. When
we spectrally shape the pump pulse, we carry out a unitary transformation on
$\alpha(\omega_{+})$, i.e. $\alpha(\omega_{+})\rightarrow
M_{1}(\omega_{+})\alpha(\omega_{+})$. Consequently, in this example we
consider $M_{1}(\omega_{+})=M(\omega_{+})$ and $M_{2}(\omega_{-})\equiv 1$, to
write (49) as

$\displaystyle|\psi(\omega_{+})|^{2}M(\omega_{+})=\sqrt{\mathcal{N}}\xi^{*}(\omega_{+})\int\xi(\omega_{+})M(\omega_{+})\rmd\omega_{+},$
(70)

with

$\xi(\omega_{+})=\int\mathcal{W}_{t}(\omega_{+},\omega_{-})\,\rmd\omega_{-}.$
(71)

We can now apply our previous results: Adopting the ansatz (52), and mirroring
equations (54)–(56), we write
$\xi(\omega_{+})=|\xi(\omega_{+})|\rme^{\rmi\vartheta(\omega_{+})}$. The
optimal pulse shaping function is then

$M(\omega_{+})=\rme^{-\rmi\vartheta(\omega_{+})},$ (72)

with the Lagrange multiplier

$|\psi(\omega_{+})|^{2}=\sqrt{\mathcal{N}}|\xi(\omega_{+})|\int|\xi(\omega_{+})|\,\rmd\omega_{+}.$
(73)

Finally, the optimization ratio (60) becomes

$E_{\mathrm{opt}}=\frac{\left(\int|\xi(\omega)|\,\rmd\omega\right)^{2}}{\left|\int\xi(\omega)\,\rmd\omega\right|^{2}}.$
(74)

As in section 4.3, to assess the usefulness of this shaping procedure we have
to make $\alpha(\omega_{+})$ and $\beta(\omega_{-})$ in (68) explicit.
Following [52], we describe the pump with a Gaussian wave packet with spectral
chirp:

$\alpha(\omega_{+})=\frac{\rme^{-(\omega_{+}-\omega_{f})^{2}/(2\sigma^{2})}}{(\pi\sigma^{2})^{1/4}}\rme^{\rmi\frac{\phi}{2}(\omega_{+}-\omega_{f})^{2}},$
(75)

where $\omega_{f}$ is the $\ket{g}\rightarrow\ket{f}$ transition frequency,
$\sigma$ the pulse width, and $\phi$ the quadratic phase determining the
chirp. To compare these parameters to the detuning $\Delta$ and deviation
$\delta$ characterizing the matter spectrum, see (17), we first consider an
“infinitely” broad phase-matching function, to avoid adding further
parameters. Once this context has been analyzed, we introduce a finite width
$\zeta$ also for $\beta(\omega_{-})$.

#### Infinitely broad phase matching

The phase-matching function describes the probability that the frequencies of
the two down-converted photons differ by $\omega_{-}=\omega_{1}-\omega_{2}$.
In the limit of infinite phase matching, then, any detuning between the two
photons is equally probable. Taking the pulse structure of figure 1(b) as an
example, we are here considering a pulse of infinite extension along the anti-
diagonal. If any difference between photon frequencies is equally probable,
$\Delta$ (see (17)), the detuning between the electronic transitions
$\ket{g}\rightarrow\ket{e}$ and $\ket{e}\rightarrow\ket{f}$, only determines
$\omega_{f}$. The incoming two-photon wave function can be approximated by
$\alpha(\omega_{+})$, which we can take out of the integral in (71), to write

$\xi^{\infty}(\omega_{+})=\alpha(\omega_{+})\eta^{\infty}(\omega_{+}),$ (76)

where we have integrated the matter response function (15) as

$\eta^{\infty}(\omega_{+})=\frac{1}{2}\int
T_{t}\left(\frac{\omega_{+}+\omega_{-}}{2},\frac{\omega_{+}-\omega_{-}}{2}\right)\,\rmd\omega_{-}=\frac{-2\rmi\sqrt{\gamma_{e}\gamma_{f}}\rme^{-\rmi\omega_{+}t}}{\omega_{+}-\omega_{f}+\rmi\gamma_{f}}.$
(77)

In figure 7(a) and (b) we plot the enhancement $E_{\mathrm{opt}}$ in this
limit, calculated using $\xi^{\infty}$ in (74). In the absence of the chirp,
panel (a), we are comparing the bandwidth of the pulse $\sigma$ to the
linewidth $\gamma_{f}=(2+\delta)\gamma_{e}$ of the two-photon transition.
Following the remarks closing section 4.2, the pulse shaping imposes the
correct phase structure, entirely encoded in $\eta^{\infty}(\omega_{+})$, to
the state $\ket{\Sigma}$, maximizing its overlap with the optimal state
$\ket{\Phi}$. The enhancement therefore increases when we move to negative
$\delta$, i.e. when the phase of $\eta^{\infty}$ rapidly changes with
$\omega_{+}$ within the bandwidth $\sigma$ of the wave packet $\alpha$.

The same reasoning can be applied when the chirp is present and induces a non-
monotonic increase of $E_{\mathrm{opt}}$ with $\sigma$ and $\delta$, see
figure 7(b). For large $\sigma$, the chirp introduces significant oscillations
in $\xi^{\infty}(\omega_{+})$: compare panels (c) and (d). When $\gamma_{f}$
is also large, for $\delta>0$ in panel (b), these oscillations also affect
$\xi^{\infty}$. In this case, then, shaping via (72) adds exactly the right
frequency dependence to the phase of $\ket{\Sigma}$ to counteract these
oscillations.

Figure 7:  Left: Optimization ratio $E_{\mathrm{opt}}$ (74), as a function of
the deviation $\delta$, obtained by shaping via (72) a pump pulse with a
profile (75) of bandwidth $\sigma$. Panel (a) shows $E_{\mathrm{opt}}$ for a
quadratic phase $\phi=0$, while (b) for $\phi=1$. Right: Comparison of the
one-photon distributions $\alpha$ (75), the integrated matter response
function $\eta^{\infty}$ (77), and their product $\xi^{\infty}$ (76). Here
$\Delta=\delta=0$, $\phi=1$ and, respectively, (c) $\sigma=0.5$ and (d)
$\sigma=5$. Continuous (dashed) lines indicate the real (imaginary) parts of
the functions.

#### Gaussian phase matching

To discuss the case of a phase-matching function with a finite width $\zeta$,
we consider the Gaussian distribution

$\beta(\omega_{-})=\frac{\rme^{-\omega_{-}^{2}/(2\zeta^{2})}}{(\pi\zeta^{2})^{1/4}}.$
(78)

With the definition $\chi=\omega_{+}-2(\omega_{e}-\rmi\gamma_{e})$ we can then
write

$\eta(\omega_{+})=\int\beta(\omega_{-})T_{t}\left(\frac{\omega_{+}+\omega_{-}}{2},\frac{\omega_{+}-\omega_{-}}{2}\right)\,\rmd\omega_{-}=2\eta^{\infty}(\omega_{+})\beta(\chi)\mathcal{G}(\rmi\chi/\zeta),$
(79)

where $\mathcal{G}$ denotes the cumulative function of the standard normal
distribution [49].

In figure 8(b,c) we fix the quadratic phase at $\phi=1$, and tune the
parameters $\sigma$ and $\zeta$ characterizing the incoming pulse (68) via,
respectively, $\alpha(\omega_{+})$ and $\beta(\omega_{-})$. The former
determines the width of the pump pulse, and hence should be compared to the
deviation $\delta$, which sets the width of the optimal two-photon wave
function (31), see figure 1(b) and (32). The latter determines the frequency
difference between the two down-converted photons, and thus should be compared
to the deviation $\Delta$ between the transition frequencies
$\ket{g}\rightarrow\ket{e}$ and $\ket{e}\rightarrow\ket{f}$, see figure 1(c)
and (33).

As in section 4.3, we adjust the pulse parameters ($\sigma,\zeta$) to the
matter characteristics ($\delta,\Delta$) by setting101010 To choose the
specific functional dependence of $\sigma$ and $\zeta$ we initially set
$\sigma=A\gamma_{f}$ and $\zeta=B\gamma_{e}(2+\Delta)$, and computed the
enhancement for several choices of $A$ and $B$. This coarse numerical analysis
showed that $A=3$ and $B=1$ are convenient values, yielding an enhancement
close to the maximal, and perfectly sufficient for the remarks we intend to
make in this example.  $\sigma=3\gamma_{f}=3\gamma_{e}(2+\delta)$ and
$\zeta=\gamma_{e}(2+\Delta)$. This allows us to explore in figure 8 the same
parameter space ($\delta,\Delta$) that was considered in section 3. In panel
(a) we show the enhancement $E_{\mathrm{q}}$ in final-state population
achieved by the optimal state (31) over the classical limit given by (36),
which was already presented in [40] and discussed in section 3.4. In panels
(b) and (c) we analogously compare the enhancement attained by, respectively,
the Gaussian pulse (68) upon shaping the pump with (72), and without shaping.
As shown in section 3.4, the ideal two-photon state (15) maximally populates
the state $\ket{f}$, and hence yields the largest $E_{\mathrm{q}}$ for any
fixed value of $\Delta$ and $\delta$. In panels (b) and (c) we also observe
$E_{\mathrm{q}}\geq 1$ in portions (on the left) of the parameter space that
we demarcate with a dashed lines. There, also Gaussian pulses are able to
enhance the TPA with respect to the yield of the optimal classical pulses. We
recall what we concluded in section 4.2: the shaped pulses can yield larger
enhancement than the unoptimized pulses because the pulse shaper (72)
perfectly counteracts the chirp.

Figure 8: Comparison of the quantum enhancement $E_{\mathrm{q}}$, with respect
to the optimal separable pulses (36), achievable via (a) the optimal entangled
two-photon wave function (15), (b) the unshaped Gaussian pulse of (68), and
(c) the optimized Gaussian via the pulse shaping function of (72). For each
$(\delta,\Delta)$ defining the matter response function (15), we attune the
Gaussian pulse by setting its widths as $\sigma=3\gamma_{e}(2+\delta)$ and
$\zeta=\gamma_{e}(2+\Delta)$ (see text). The quadratic phase in (75) is
$\phi=1$. On the left of the dashed line the various two-photon wave functions
perform better ($E_{\mathrm{q}}\geq 1$) than the optimal separable pulses
(36). Panel (a) reproduces the numerical results of [40].

##  5 Conclusions

To conclude, we analyzed optimal two-photon states to drive a two-photon
transition. We first investigated the optimal state, quantified its quantum
correlations via the entanglement entropy, and discussed the relation of the
latter to the quantum enhancement in TPA such a state can achieve. We then
drew a comparison to the optimal separable pulse, computed from the Schmidt
decomposition of the response function, such that the enhancement really stems
from the entanglement in the state. For the maximally achievable enhancement
we also provided bounds that depend on how strongly anticorrelated the joint
distribution of the frequencies is.

We then considered more realistic scenarios where a given initial two-photon
state is manipulated in order to enhance the two-photon transition. We defined
a new optimization problem, this time for unitary operators representing local
transformations of the individual photons. We derived a self-consistent
equation that can be solved analytically when we consider two photons created
in spontaneous parametric down-conversion. In particular, we inspected the
case of spatial light modulators to shape two photons converted from a
monochromatic pump laser, and then considered shaping the pump pulse directly.
We found that initial two-photon states with sufficiently strong entanglement
can sustain a substantial enhancement of the TPA probability, of the same
order of magnitude as what is ideally achievable.

##  Appendix A Maximization of (46)

Let us write all the terms appearing in (46) in the frequency representation,
by resolving the identity in the appropriate one- or two-photon frequency
space, see section 2.5. If we call
$\mathcal{A}[M_{1},M_{2}]=\braket{\Phi}{M_{1}M_{2}}{\Sigma}$, we then have

$\mathcal{A}[M_{1},M_{2}]=\iint_{\mathbb{R}^{2}}\rmd\omega_{1}\rmd\omega_{2}\iint_{\mathbb{R}^{2}}\rmd\nu_{1}\rmd\nu_{2}\,\Phi^{*}(\omega_{1},\omega_{2})M_{1}(\nu_{1},\omega_{1})M_{2}(\nu_{2},\omega_{2})\Sigma(\nu_{1},\nu_{2}).$
(80)

Similarly, we define

$\mathcal{B}_{j}[M_{j}]=\braket{\psi_{j}}{M_{j}^{\dagger}M_{j}}{\psi_{j}}=\iiint_{\mathbb{R}^{3}}\rmd\omega_{j}\rmd\nu_{j}\rmd\mu_{j}\,\psi_{j}^{*}(\omega_{j})M_{j}^{*}(\nu_{j},\omega_{j})M_{j}(\mu_{j},\nu_{j})\psi_{j}(\mu_{j})$
(81)

and

$\mathcal{C}_{j}=\int_{\mathbb{R}}\rmd\omega_{j}\psi_{j}^{*}(\omega_{j})\psi_{j}(\omega_{j}).$
(82)

We remind the reader that $\mathcal{B}_{j}$ and $\mathcal{C}_{j}$ are also
functionals of $\psi_{j}(\omega_{j})$ and $\psi^{*}_{j}(\omega_{j})$, the wave
function (and corresponding conjugate) representing $\ket{\psi_{j}}$. Since
the latter play the role of Lagrange multipliers, however, we do not write
them as arguments on the left-hand side of (81,82).

We can then rewrite (46) in the frequency representation as

$J[M_{1},M_{2}]=\mathcal{N}\mathcal{A}[M_{1},M_{2}]\mathcal{A}^{*}[M^{*}_{1},M^{*}_{2}]-\sum_{j=1}^{2}\\{\mathcal{B}_{j}[M_{j}]-\mathcal{C}_{j}\\},$
(83)

where the normalization $\mathcal{C}_{j}$ of the $\ket{\psi_{j}}$ states,
being a number, does not carry any dependence on the pulse shaping functions,
over which we are optimizing.

To find the solution to (83), we require its functional derivatives with
respect to the pulse shaping function $M^{*}_{j}$ (to obtain an expression for
$M_{j}$) and to the Lagrange multipliers, respectively, to vanish:

$\displaystyle\frac{\delta J}{\delta M_{j}^{*}}=0,$ (84)

and

$\displaystyle\frac{\delta J}{\delta\psi_{j}^{*}}=0.$ (85)

Equation (85) enforces the unitarity (47), while (84) yields the general
integral equation

$\psi_{j}^{*}(\omega_{j})\int_{\mathbb{R}}\rmd\mu_{j}\,M_{j}(\mu_{j},\nu_{j})\psi_{j}(\mu_{j})=\mathcal{N}\mathcal{A}[M_{1},M_{2}]\frac{\delta\mathcal{A}^{*}}{\delta
M_{j}^{*}}(\nu_{j},\omega_{j}),$ (86)

where $j,k\in\\{1,2\\},j\neq k$, and

$\frac{\delta\mathcal{A}^{*}}{\delta
M_{j}^{*}}(\nu_{j},\omega_{j})=\iint_{\mathbb{R}^{2}}\rmd\omega_{k}\rmd\nu_{k}\,\Sigma^{*}(\omega_{j},\omega_{k})M_{k}^{*}(\nu_{k},\omega_{k})\Phi(\nu_{j},\nu_{k}).$
(87)

This work is supported by the European Research Council under the European
Union’s Seventh Framework Programme (FP7/2007-2013) Grant Agreement No. 319286
Q-MAC. E.G.C. acknowledges support from the Georg H. Endress foundation. F.S.
acknowledges support from the Cluster of Excellence ’Advanced Imaging of
Matter’ of the Deutsche Forschungsgemeinschaft (DFG) - EXC 2056 - project ID
390715994.

## References

## References

  * [1] Rabitz H, de Vivie-Riedle R, Motzkus M and Kompa K 2000 Science 288 824
  * [2] Brif C, Chakrabarti R and Rabitz H 2010 New J. Phys. 12 075008
  * [3] Glaser S J, Boscain U, Calarco T, Koch C P, Köckenberger W, Kosloff R, Kuprov I, Luy B, Schirmer S, Schulte-Herbrüggen T, Sugny D and Wilhelm F K 2015 Eur. Phys. J. D 69 279
  * [4] Tannor D J, Kosloff R and Rice S A 1986 J. Chem. Phys. 85 5805
  * [5] Peirce A P, Dahleh M A and Rabitz H 1988 Phys. Rev. A 37 4950
  * [6] Kosloff R, Rice S A, Gaspard P, Tersigni S and Tannor D J 1989 Chem. Phys. 139 201
  * [7] Judson R S and Rabitz H 1992 Phys. Rev. Lett. 68 1500
  * [8] Assion A 1998 Science 282 919
  * [9] Prokhorenko V I, Nagy A M, Waschuk S A, Brown L S, Birge R R and Miller R J D 2006 Science 313 1257
  * [10] Wohlleben W, Buckup T, Herek J L and Motzkus M 2005 ChemPhysChem 6 850
  * [11] Silberberg Y 2009 Annu. Rev. Phys. Chem. 60 277
  * [12] Bartana A, Kosloff R and Tannor D J 1993 J. Chem. Phys. 99 196
  * [13] Bartana A, Kosloff R and Tannor D J 1997 J. Chem. Phys. 106 1435
  * [14] Frimmer M, Gieseler J and Novotny L 2016 Phys. Rev. Lett. 117 163601
  * [15] Platzer F, Mintert F and Buchleitner A 2010 Phys. Rev. Lett. 105 020501
  * [16] Schäfer V M, Ballance C J, Thirumalai K, Stephenson L J, Ballance T G, Steane A M and Lucas D M 2018 Nature 555 75
  * [17] Shapiro M and Brumer P 2012 Quantum Control of Molecular Processes (Hoboken, USA: Wiley)
  * [18] Dorfman K E, Schlawin F and Mukamel S 2016 Rev. Mod. Phys. 88 045008
  * [19] Schlawin F, Dorfman K E and Mukamel S 2018 Acc. Chem. Res. 51 2207
  * [20] Gilaberte Basset M, Setzpfandt F, Steinlechner F, Beckert E, Pertsch T and Gräfe M 2019 Laser Photonics Rev. 13 1900097
  * [21] Mukamel S, Freyberger M, Schleich W, Bellini M, Zavatta A, Leuchs G, Silberhorn C, Boyd R W, Sánchez-Soto L L, Stefanov A, Barbieri M, Paterova A, Krivitsky L, Shwartz S, Tamasaku K, Dorfman K, Schlawin F, Sandoghdar V, Raymer M, Marcus A, Varnavski O, Goodson T, Zhou Z Y, Shi B S, Asban S, Scully M, Agarwal G, Peng T, Sokolov A V, Zhang Z D, Zubairy M S, Vartanyants I A, del Valle E and Laussy F 2020 J. Phys. B At. Mol. Opt. Phys. 53 072002
  * [22] Szoke S, Liu H, Hickam B P, He M and Cushing S K 2020 J. Mater. Chem. C 8 10732
  * [23] Raymer M G, Landes T and Marcus A H 2021 Entangled two-photon absorption by atoms and molecules: A quantum optics tutorial (Preprint 2103.02551)
  * [24] Georgiades N P, Polzik E S, Edamatsu K, Kimble H J and Parkins A S 1995 Phys. Rev. Lett. 75 3426
  * [25] Dayan B, Pe’er A, Friesem A A and Silberberg Y 2004 Phys. Rev. Lett. 93 023005
  * [26] Lee D I and Goodson T 2006 J. Phys. Chem. B 110 25582
  * [27] Upton L, Harpham M, Suzer O, Richter M, Mukamel S and Goodson T 2013 J. Phys. Chem. Lett. 4 2046
  * [28] Schlawin F, Dorfman K E, Fingerhut B P and Mukamel S 2013 Nat. Commun. 4 1782
  * [29] Raymer M G, Marcus A H, Widom J R and Vitullo D L 2013 J. Phys. Chem. B 117 15559
  * [30] León-Montiel R d J, Svozilík J c v, Torres J P and U’Ren A B 2019 Phys. Rev. Lett. 123 023601
  * [31] Tabakaev D, Montagnese M, Haack G, Bonacina L, Wolf J P, Zbinden H and Thew R T 2021 Phys. Rev. A 103 033701
  * [32] Parzuchowski K M 2021 Phys. Rev. Appl. 15
  * [33] Landes T, Allgaier M, Merkouche S, Smith B J, Marcus A H and Raymer M G 2020 (Preprint 2012.06736)
  * [34] Raymer M G, Landes T, Allgaier M, Merkouche S, Smith B J and Marcus A H 2020 (Preprint 2012.05375)
  * [35] Pe’er A, Dayan B, Friesem A A and Silberberg Y 2005 Phys. Rev. Lett. 94 1
  * [36] Bernhard C, Bessire B, Feurer T and Stefanov A 2013 Phys. Rev. A 88 032322
  * [37] Bessire B, Bernhard C, Feurer T and Stefanov A 2014 New J. Phys. 16 033017
  * [38] Castro A, Appel H and Rubio A 2019 Eur. Phys. J. B 92 223
  * [39] Oka H 2018 Phys. Rev. A 97 1
  * [40] Schlawin F and Buchleitner A 2017 New J. Phys. 19 013009
  * [41] Csehi A, Vibók A, Halász G J and Kowalewski M 2019 Phys. Rev. A 100 053421
  * [42] Loudon R 2000 The quantum theory of light 3rd ed (Oxford: Oxford University Press)
  * [43] Schlawin F 2017 J. Phys. B At. Mol. Opt. Phys. 50 203001
  * [44] Scheck F 2010 Mechanics (Springer Berlin Heidelberg)
  * [45] Law C K, Walmsley I A and Eberly J H 2000 Phys. Rev. Lett. 84 5304
  * [46] Stobińska M, Alber G and Leuchs G 2009 EPL 86 14007
  * [47] Aljunid S A, Maslennikov G, Wang Y, Dao H L, Scarani V and Kurtsiefer C 2013 Phys. Rev. Lett. 111 103001
  * [48] Brecht B, Reddy D V, Silberhorn C and Raymer M G 2015 Phys. Rev. X 5 041017
  * [49] Bronshtein I N, Semendyayev K A, Musiol G and Mühlig H 2007 Handbook of Mathematics 5th ed (Springer Berlin Heidelberg)
  * [50] Zheng Z, Saldanha P L, Rios Leite J R and Fabre C 2013 Phys. Rev. A 88 033822
  * [51] Shih Y 2003 Reports on Progress in Physics 66 1009
  * [52] Arzani F, Fabre C and Treps N 2018 Phys. Rev. A 97 033808
  * [53] Keller T E and Rubin M H 1997 Phys. Rev. A 56 1534
  * [54] Nambu Y, Usami K, Tsuda Y, Matsumoto K and Nakamura K 2002 Phys. Rev. A 66 033816
  * [55] U’Ren A B, Banaszek K and Walmsley I A 2003 Quantum Info. Comput. 3 480
  * [56] Bengtsson I and Zyczkowski K 2006 Geometry of Quantum States (Cambridge University Press)