# Optimization of selective two-photon absorption in cavity polaritons

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

(August 28, 2024)

###### Abstract

We investigate optimal states of photon pairs to excite a target transition in
a multilevel quantum system. With the help of coherent control theory for two-
photon absorption with quantum light, we infer the maximal population
achievable by optimal entangled vs. separable states of light. Interference
between excitation pathways, as well as the presence of nearby states, may
hamper the selective excitation of a particular target state, but we show that
quantum correlations can help to overcome this problem, and enhance the
achievable “selectivity” between two energy levels, i.e. the relative
difference in population transferred into each of them. We find that the added
value of optimal entangled states of light increases with broadening
linewidths of the target states.

##  I Introduction

The theoretical description of (nonlinear) spectroscopy is conventionally
based on a semiclassical approach, where the light fields are treated
classically and only the sample system is treated fully quantum mechanically
Mukamel (1995); Hamm and Zanni (2011). In most situations, this approximation
is extremely well justified, owing to the weak nonlinearity of light-matter
interactions in free space Boyd and Boyd (1992). Nevertheless, recent years
have seen the rapid rise of theoretical investigations, as well as first
proof-of-concept experiments that challenge this convention, and investigate
how quantum properties of light can be applied or exploited beneficially in
spectroscopic applications Dorfman, Schlawin, and Mukamel (2016); Schlawin,
Dorfman, and Mukamel (2018); Gilaberte Basset  _et al._ (2019); Mukamel  _et
al._ (2020); Szoke  _et al._ (2020). This includes the use of photon
correlation measurements to analyse the light fields emitted by a sample del
Valle  _et al._ (2012); Gonzalez-Tudela  _et al._ (2013); Wientjes  _et al._
(2014); Holdaway, Notararigo, and Olaya-Castro (2018); Sánchez Muñoz and
Schlawin (2020) in single-molecule spectroscopy Krüger  _et al._ (2010);
Hildner  _et al._ (2013); Kyeyune  _et al._ (2019), or to exploit coincidence
measurements to detect particular spectral features in the sample Dorfman and
Mukamel (2012); Schlawin, Dorfman, and Mukamel (2016); Dorfman and Mukamel
(2018); Zhang  _et al._ (2018), as well as the generation of photonic
entanglement in fluorescent proteins Shi, Kumar, and Lee (2017).

However, the arguably most active field of research concerns the use of
quantum light, and in particular of entangled photons, to excite the sample.
Squeezed states can improve linear absorption measurements Li  _et al._
(2020a, b). Nonlinear optical signals such as two-photon absorption scale
linearly with the photon flux Georgiades  _et al._ (1995); Dayan  _et al._
(2004); Lee and Goodson (2006), which could enable nonlinear spectroscopy of
photosensitive samples at very low intensities. Following earlier experiments
on two-photon absorption in biomolecules Lee and Goodson (2006); Upton  _et
al._ (2013), a series of recent experiments have scrutinised the situation
Parzuchowski  _et al._ (2020); Landes  _et al._ (2020); Tabakaev  _et al._
(2021) and report widely differing entangled two-photon absorption cross
sections. They inspired new theoretical investigations into the enhancement
that entanglement can provide in two-photon absorption Raymer  _et al._
(2020); Raymer, Landes, and Marcus (2021). Apart from the linear scaling,
theoretical proposals show that spectral quantum correlations of entangled
photon pairs could further help to disentangle complex optical signals and
reveal otherwise hidden features Schlawin  _et al._ (2013); Raymer  _et al._
(2013); León-Montiel  _et al._ (2019), and enable ultrafast spectroscopy in a
cw setup Ishizaki (2020). They could also be used for the generation of
pseudo-sunlight to imitate natural conditions Fujihashi, Shimizu, and Ishizaki
(2020) and provide sensitive probes for dynamical symmetry breaking Li  _et
al._ (2017) and many-body correlations Li  _et al._ (2019). The control of
quantum correlations using temperature León-Montiel  _et al._ (2019), photon
statistics Castro, Appel, and Rubio (2019); Csehi  _et al._ (2019) or spectral
shaping Oka (2018) could further enhance these beneficial properties and
provide experimentalists with new handles to manipulate optical signals in a
way that is not possible in laser-based spectroscopies.

One pertinent question regarding the application of quantum light to
spectroscopy is how much the quantum nature of the former can enhance a given
spectroscopic task. To ultimately decide this in an unbiased way, it is
necessary to compare the performance with optimised (classical) laser pulses,
i.e. with the optimal performance achievable by classical means. The optimal
control with shaped classical laser pulses is a well established field of
research Rabitz  _et al._ (2000); Brif, Chakrabarti, and Rabitz (2010). In
particular, optimal control of two-photon absorption was described and
implemented experimentally in the late 1990’s Dubov and Rabitz (1996); Assion
(1998); Meshulach and Silberberg (1998, 1999); Dudovich  _et al._ (2001). In
these applications, which typically rely on strong laser fields that can
manipulate the interference between excitation pathways, quantum light is seen
as rather detrimental, as pointed out in Ref. Shapiro and Brumer, 2012.
However, we recently showed that this is not true for weak broadband fields,
where quantum correlations of light can enhance excitation probabilities
Schlawin and Buchleitner (2017); Carnio, Buchleitner, and Schlawin (2021).

To find these optimal quantum states of light to drive a two-photon
transition, we used a coherent control theory for continuous-mode quantum
light. In particular, this enabled us to quantify the possible enhancement of
the two-photon absorption probability due to quantum correlations between
frequency components of the injected pulses, which go hand in hand with strong
correlations of the arrival time of the photons. We have so far restricted
this analysis to a simple three-level system. In this paper, we extend this
theory to two-photon excitations in multilevel systems.

We connect to our previous work by briefly recalling in Sec. II the
theoretical framework of light-matter interactions. In Sec. III, we recollect
the theory of optimal states driving a three-level system, and then generalize
this to multilevel targets, where the excitation of undesired, nearby states
can prevent that of the particular target states. Finally, in Sec. IV we apply
the formalism to a cavity polariton system, where an interesting interference
effect involving entangled photons was recently reported in Ref. Gu and
Mukamel, 2020. We conclude with Sec. V.

##  II Theoretical framework

We consider two pulses – each carrying a single photon – impinging on the
atomic target. These fields and matter degrees of freedom are described,
respectively, by Hamiltonians $H_{\text{f}}$ and $H_{\text{m}}$, and are
coupled by a light-matter interaction Hamiltonian $W$, which we present in the
following paragraphs. The total Hamiltonian thus reads
$H=H_{\text{f}}+H_{\text{m}}+W$.

Both quantized pulses are injected along fixed and distinct spatial
directions, and are described by continuous-mode operators, $a_{1}(\omega)$
and $a_{2}(\omega)$, respectivelyLoudon (2000). The fields’ Hamiltonian,
ignoring the vacuum energy, is then
$H_{\text{f}}=\sum_{j}\int_{0}^{\infty}\mathrm{d}\omega\,\hbar\omega
a_{j}^{\dagger}(\omega)a_{j}(\omega)$. The positive-frequency electric field
operator acting on the Hilbert space of photon $j$ (in the interaction picture
with respect to $H_{\text{f}}$) reads

$E_{j}^{+}(z,t)=\mathrm{i}\int_{0}^{\infty}\\!\\!\mathrm{d}\omega\;\sqrt{\frac{\hbar\omega}{4\pi\epsilon_{0}cA}}a_{j}(\omega)\mathrm{e}^{\mathrm{i}[k(\omega)z-\omega
t]}.$ (1)

Here, $z$ is the position along the propagation direction, $A$ the
quantization area perpendicular to it, and $c$ the speed of light, and we
assume a parallel polarization of the pulses. We consider the target sample
placed at $z=0$, and much smaller than the wavelength of the light field, such
that we can drop the spatial modulation of the field operator (1).
Furthermore, we only consider field states characterized by narrow pulse
shapes, of bandwidth $\Delta\omega$, distributed around a central frequency
$\omega_{0}\gg\Delta\omega$. Since all expectation values are calculated with
respect to these states, we can safely extend the range of frequency
integration and write the electric field operator (1) in the _narrow bandwidth
approximation_ Loudon (2000) as

$E_{j}^{+}(t)=\mathrm{i}\mathcal{E}_{0}\int_{-\infty}^{\infty}\\!\\!\mathrm{d}\omega\;a_{j}(\omega)\mathrm{e}^{-\mathrm{i}\omega
t},$ (2)

where $\mathcal{E}_{0}=(\hbar\omega_{0}/4\pi\epsilon_{0}cA)^{1/2}$
approximates the field normalization of (1). The range of integration
$(-\infty,\infty)$ will be assumed from now on.

As for the matter degrees of freedom, we consider a system with a ground state
$\ket{g}$, multiple intermediate states $\ket{e_{1}},\ket{e_{2}},\ldots$ and
multiple final states $\ket{f_{1}},\ket{f_{2}},\ldots$. With each state
$\ket{s}$ we associate an energy $\hbar\omega_{s}$ (we set
$\hbar\omega_{g}=0$), an inverse lifetime $\gamma_{s}$ ($\gamma_{g}=0$ since
the ground state cannot decay), and a Lorentzian line shape
$\mathcal{L}_{s}(\omega)=(\omega-z_{s})^{-1}$, where
$z_{s}=\omega_{s}-\mathrm{i}\gamma_{s}$. For eigenstates with a finite
lifetime, we consider an effective, non-HermitianBreuer and Petruccione (2007)
matter Hamiltonian $H_{\text{m}}=\sum_{j}\hbar
z_{e_{j}}\ket{e_{j}}\bra{e_{j}}+\sum_{k}\hbar
z_{f_{k}}\ket{f_{k}}\bra{f_{k}}$. Adjacent manifolds are dipole-coupled with
dipole matrix elements (along the fields’ polarization) $\mu_{ge_{j}}$ and
$\mu_{e_{j}f_{k}}$, respectively.

The light-matter interaction Hamiltonian in the rotating wave approximation –
which is certainly valid in the present near-resonant perturbative regime
where $\mu\mathcal{E}_{0}\ll\hbar\omega_{0}$ (for any of the dipole matrix
elements $\mu$ above) – and, as the subscript indicates, in the interaction
picture with respect to $H_{\text{f}}+H_{\text{m}}$, is therefore given by
Loudon (2000)

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

where the components of the dipole operator which annihilate an electronic
excitation may be written as

$\displaystyle V(t)=$
$\displaystyle\sum_{j,k}\mu_{ge_{j}}\mathrm{e}^{-\mathrm{i}z_{e_{j}}t}\ket{g}\bra{e_{j}}$
$\displaystyle+\mu_{e_{j}f_{k}}\mathrm{e}^{-\mathrm{i}(z_{f_{k}}-z_{e_{j}})t}\ket{e_{j}}\bra{f_{k}},$
(4)

and $E^{+(-)}=E^{+(-)}_{1}+E^{+(-)}_{2}$, since we assume that each pulse
couples identically to the matter degrees of freedom.

##  III Optimization procedure

###  III.1 Single final state

We seek a two-photon state $\ket{\Phi_{f}}$ of the incoming fields which
optimizes the two-photon transition from the ground state $\ket{g}$ to a
single final state $\ket{f}$, via a manifold of $n_{e}$ intermediate states
$\ket{e_{j}}$.

Mathematically, we consider an initial state $\ket{\Phi_{f}}\otimes\ket{g}$ in
the combined matter-field Hilbert space
$\mathcal{H}=\mathcal{H}_{\text{f}}\otimes\mathcal{H}_{\text{m}}$ (mirroring
the subscripts of the Hamiltonians in the previous section), and determine
$\ket{\Phi_{f}}$ such that the time evolution operator $U_{I}(t,t_{0})$ in the
interaction picture of $H_{\text{f}}+H_{\text{m}}$ maximises the population of
the target state $\ket{0}\otimes\ket{f}$, where both photons have been
absorbed in order to drive the matter degrees of freedom to the final state
$\ket{f}$. As described in detail in Ref. Carnio, Buchleitner, and Schlawin,
2021, we can find this transition amplitude perturbatively. Tracing out the
matter degrees of freedom, the desired transition is mediated by an operator
acting on $\mathcal{H}_{\text{f}}$ alone,

$\displaystyle\hat{T}_{f}$ $\displaystyle=$ $\displaystyle\langle
f(t)|U_{I}(t,t_{0})|g\rangle$ (5) $\displaystyle=$ $\displaystyle\iint
T_{f}(\omega_{1},\omega_{2})a_{1}(\omega_{1})a_{2}(\omega_{2})\mathrm{d}\omega_{1}\mathrm{d}\omega_{2},$

with the explicit expression of the matter response function

$\displaystyle T_{f}(\omega_{1},\omega_{2})=$
$\displaystyle\left(\frac{\mathrm{i}\mathcal{E}_{0}}{\hbar}\right)^{2}\sum_{j}\mu_{ge_{j}}\mu_{e_{j}f}\left[\mathcal{L}_{e_{j}}(\omega_{1})+\mathcal{L}_{e_{j}}(\omega_{2})\right]$
$\displaystyle\times\mathcal{L}_{f}(\omega_{1}+\omega_{2})\mathrm{e}^{-\mathrm{i}(\omega_{1}+\omega_{2})t}.$
(6)

The matter response function (III.1) is to be derived for the excitation of
the population of $\ket{f}$ at a fixed time $t$, with the interaction
$W_{I}(t)$ turned on at $t_{0}\rightarrow-\infty$. Without loss of generality,
we can therefore set $t=0$. The state $\ket{\Phi_{f}}$ can now be found by
variation of the functional Schlawin and Buchleitner (2017); Carnio,
Buchleitner, and Schlawin (2021)

$J_{\mathrm{single}}[\ket{\Phi}]=p_{f}-\lambda\left(\braket{\Phi}{\Phi}-1\right),$
(7)

where $p_{f}=|\braket{0}{\hat{T}_{f}}{\Phi}|^{2}$ is the population in the
target state $\ket{f}$ created by the absorption of $\ket{\Phi}$, leaving the
field in the vacuum state $\ket{0}$. The Lagrange multiplier $\lambda$
constrains the optimization to normalized states111We remark here that, if we
restrict our optimization problem to two-photon states, with one photon per
mode, as we do here, the normalization constraint
$\braket{\Phi}{n_{1}n_{2}}{\Phi}$ used in Schlawin and Buchleitner, 2017;
Carnio, Buchleitner, and Schlawin, 2021, where $n_{j}$ is the photon number
operator for beam $j$, corresponds to $\braket{\Phi}{\Phi}$ used here..
Defining the (unnormalized) two-photon state
$\ket{T_{f}}=\hat{T}^{\dagger}_{f}\ket{0}$ via (5), the functional (7) is
maximizedSchlawin and Buchleitner (2017); Carnio, Buchleitner, and Schlawin
(2021) by the state $\ket{\Phi_{f}}=\mathcal{N}_{f}^{-1/2}\ket{T_{f}}$, whose
associated _wave function_ reads

$\Phi_{f}(\omega_{1},\omega_{2})=\braket{0}{a_{1}(\omega_{1})a_{2}(\omega_{2})}{\Phi_{f}}=\frac{T_{f}^{*}(\omega_{1},\omega_{2})}{\sqrt{\mathcal{N}_{f}}},$
(8)

given the appropriate normalization $\mathcal{N}_{f}=\braket{T_{f}}{T_{f}}$,
which we determine analytically in the next subsection.

###  III.2 Multiple final states

We now generalize the formalism to the case of a manifold with $n_{f}$ target
states. Our aim is to find the optimal two-photon state that maximally
populates a given final state $\ket{f_{1}}$, while minimizing the population
of all other energetically near-degenerate states
$\ket{f_{2}},\ket{f_{3}},\ldots$, which are equally reachable in terms of the
energy of the incoming radiation (and assuming that no selection rule prevents
this). To this end, we need to generalize the functional (7) such that the
target state’s population $p_{f_{1}}$ is maximized, while the excitation of
any other states $\ket{f_{j\neq 1}}$ is penalized:

$J[\ket{\Phi}]=p_{f_{1}}-\sum_{j\neq
1}p_{f_{j}}-\lambda\left(\braket{\Phi}{\Phi}-1\right).$ (9)

Notice that enforcing strictly vanishing populations in the states
$\ket{f_{j}}\neq\ket{f_{1}}$ is prevented by the fact that the matter response
functions (III.1) for different target states $\ket{f_{j}}$ in general exhibit
finite overlap.

In analogy to the previous section, to maximize (9) we first define two-photon
states $\ket{T_{f_{j}}}=\hat{T}_{f_{j}}^{\dagger}\ket{0}$. We remind here that
these states are not normalized, but their normalization will be taken care of
in (12). By writing $p_{f_{j}}=|\braket{T_{f_{j}}}{\Phi}|^{2}$ the functional
(9) transforms into

$J[\ket{\Phi}]=|\braket{T_{f_{1}}}{\Phi}|^{2}-\sum_{j\neq
1}|\braket{T_{f_{j}}}{\Phi}|^{2}-\lambda\left(\braket{\Phi}{\Phi}-1\right).$
(10)

The state $\ket{\tilde{\Phi}}$ that maximizes this functional is found by
requiring the variational derivative with respect to the dual stateCarnio,
Buchleitner, and Schlawin (2021) $\bra{\Phi}$ to vanish:

$\frac{\delta J}{\delta\bra{\Phi}}=0\iff\left(K_{f_{1}}-\sum_{j\neq
1}K_{f_{j}}\right)\ket{\Phi}=\lambda\ket{\Phi},$ (11)

with $K_{f_{j}}=\ket{T_{f_{j}}}\bra{T_{f_{j}}}$. A direct and robust way to
solve this eigenvalue problem – compared, e.g., to the introduction of an
orthonormal basis via an orthogonalization procedure – is to formulate it on
the subspace spanned by the non-orthogonal states $\\{\ket{T_{f_{j}}}\\}_{j}$.
We then obtain the generalized eigenvalue problemBronshtein  _et al._ (2007)
(GEP)

$\operatorname{diag}(1,-1,-1,\ldots)\cdot\boldsymbol{v}=\lambda\,M\cdot\boldsymbol{v}\quad
M_{jk}=\frac{\Sigma_{{f_{j}}{f_{k}}}}{\sqrt{\mathcal{N}_{f_{j}}\mathcal{N}_{f_{k}}}},$
(12)

where the matrix $M$ is given by the overlaps

$\displaystyle\Sigma_{f_{j}f_{k}}$
$\displaystyle=\braket{T_{f_{j}}}{T_{f_{k}}}$
$\displaystyle=-8\pi^{2}\left(\frac{\mathcal{E}_{0}}{\hbar}\right)^{4}\sum_{m,n}\frac{\mu_{ge_{m}}\mu_{e_{m}f_{j}}\mu_{ge_{n}}\mu_{e_{n}f_{k}}}{(z_{e_{m}}-z^{*}_{e_{n}})(z_{f}-z^{*}_{f^{\prime}})},$
(13)

and the normalization constants

$\mathcal{N}_{f}=\Sigma_{ff}=\frac{8\pi^{2}\mathcal{E}_{0}^{4}}{\hbar^{4}\gamma_{f}}\sum_{m<n}\frac{\mu_{ge_{m}}\mu_{e_{m}f}\mu_{ge_{n}}\mu_{e_{n}f}(\gamma_{e_{m}}+\gamma_{e_{n}})}{(\omega_{e_{m}}-\omega_{e_{n}})^{2}+(\gamma_{e_{m}}+\gamma_{e_{n}})^{2}}.$
(14)

We remark here that the original optimization problem (9) on the space of
square-integrable functions on $\mathbb{R}^{2}$ has been reduced to an
eigenvalue problem of dimension $n_{f}$ on the target state manifold. All
matrices which enter (12) are known analytically, and exclusively depend on
the spectral properties (eigenenergies, dipole matrix elements, lifetimes) of
the matter degrees of freedom (the factors $\mathcal{E}_{0}/\hbar$ drop out).

The _selective_ optimal two-photon wave function
$\tilde{\Phi}(\omega_{1},\omega_{2})$ is determined by the eigenvector
$\tilde{\boldsymbol{v}}$ associated with the largest eigenvalue of (12), since
we seek to maximize (9). Its components provide the optimal linear combination
of the _indistinctive_ optimal two-photon wave functions $\Phi_{f_{j}}$ from
(8):

$\tilde{\Phi}(\omega_{1},\omega_{2})=\sum_{j}\tilde{v}_{j}\frac{T^{*}_{f_{j}}(\omega_{1},\omega_{2})}{\sqrt{\mathcal{N}_{f_{j}}}}.$
(15)

From the eigenvectors we can also immediately compute the maximal final-state
populations excited by the selecive optimal state:

$p_{f_{k}}=|(M\cdot\tilde{\boldsymbol{v}})_{k}|^{2},$ (16)

where we remind the reader that the eigenvectors of a GEP are orthonormal with
respect to the scalar product induced by the $M$ matrix:
$(\boldsymbol{v}_{j},\boldsymbol{v}_{k})_{M}=\boldsymbol{v}_{j}^{*}\cdot
M\cdot\boldsymbol{v}_{k}=\delta_{jk}$.

##  IV Application to a near-degenerate manifold

The coherent superposition (15) of indistinctive optimal states (8) suggests
that interference effects between the latter might play an important role in
maximizing the target state’s population $p_{f_{1}}$. However, this
information must be encoded in the coefficients $\tilde{v}_{j}$, whose
specific value, even if analytically computable, might not prove insightful.
Therefore, to assess the effectiveness of our optimization method in
selectively driving a pre-defined two-photon transition, we want to apply it
to a manifold with both near-degenerate and non-degenerate states, as a
benchmark. Such a structure is given, e.g., by the atom-field Hamiltonian
describing $N$ atoms “dressed” by a quantized cavity modeCohen-Tannoudji,
Dupont-Roc, and Grynberg (1998). After defining the model, we study how the
structure of the optimal two-photon states reflects the different _excitation
pathways_ $\ket{g}\rightarrow\ket{e_{j}}\rightarrow\ket{f_{k}}$, i.e. the
possible transitions, through the intermediate manifold of $\ket{e_{j}}$
states, that maximize the yield of the target two-photon transition
$\ket{g}\rightarrow\ket{f_{k}}$. In particular, we discuss the _selectivity_
with which a specific state, from a near-degenerate pair, can be excited by
entangled or separable states of light.

$\omega_{r}/\omega_{0}$ | $-2.04\text{\times}{10}^{-02}$ | $0.698$ | $0.979$ | $1.26\text{\times}{10}^{00}$ | $1.58\text{\times}{10}^{00}$ | $1.97\text{\times}{10}^{00}$ | $2.00\text{\times}{10}^{00}$ | $2.37\text{\times}{10}^{00}$
---|---|---|---|---|---|---|---|---
$\mu_{rs}$ | $g$ | $e_{1}$ | $e_{2}$ | $e_{3}$ | $f_{1}$ | $f_{2}$ | $f_{3}$ | $f_{4}$
$g$ | 0 | $1.09\text{\times}{10}^{00}$ | $1.61\text{\times}{10}^{-03}$ | $0.939$ | 0 | 0 | 0 | 0
$e_{1}$ | $1.09\text{\times}{10}^{00}$ | 0 | 0 | 0 | $0.891$ | $0.475$ | $0.705$ | $0.125$
$e_{2}$ | $1.61\text{\times}{10}^{-03}$ | 0 | 0 | 0 | $0.757$ | $-7.95\text{\times}{10}^{-02}$ | $-7.34\text{\times}{10}^{-02}$ | $0.688$
$e_{3}$ | $0.939$ | 0 | 0 | 0 | $-0.205$ | $0.515$ | $0.706$ | $-0.779$
$f_{1}$ | 0 | 0.891 | 0.757 | $-0.205$ | 0 | 0 | 0 | 0
$f_{2}$ | 0 | 0.475 | $-7.95\text{\times}{10}^{-02}$ | 0.515 | 0 | 0 | 0 | 0
$f_{3}$ | 0 | 0.705 | $-7.34\text{\times}{10}^{-02}$ | 0.706 | 0 | 0 | 0 | 0
$f_{4}$ | 0 | 0.125 | 0.688 | $-0.779$ | 0 | 0 | 0 | 0
Table 1: Eigenenergies $\omega_{r}$ (in units of the cavity mode frequency
$\omega_{0}$, with $\hbar\equiv 1$) of the dressed states $\ket{g}$,
$\ket{e_{j}}$, $\ket{f_{k}}$, $j=1,2,3$, $k=1,\ldots,4$, of the Hamiltonian
(17) for $N=2$ two-level systems with eigenfrequencies
$\omega_{1}=0.8\omega_{0}$ and $\omega_{2}=1.2\omega_{0}$, respectively, and
their mutual dipole matrix elements $\mu_{rs}$. Due to the weak coupling
$g_{1},g_{2}=$0.14\text{\,}\mathrm{eV}$\ll\omega_{0}=$1\text{\,}\mathrm{eV}$\simeq\omega_{1}\simeq\omega_{2}$,
the spectral structure of the uncoupled ($g_{1}=g_{2}=0$) dressed system – one
non-degenerate ground state, three degenerate single-excitation states, four
degenerate double-excitation states – remain essentially intact, apart from
the degeneracies being lifted.

###  IV.1 Dressed-state Hamiltonian

The $N$-atom-field Hamiltonian we work with reads

$H=\sum_{n=1}^{N}\omega_{n}\sigma_{n}^{\dagger}\sigma_{n}+\omega_{0}b^{\dagger}b+\sum_{n=1}^{N}g_{n}(\sigma_{n}^{\dagger}+\sigma_{n})(b+b^{\dagger}),$
(17)

with $\hbar\equiv 1$. We consider $N=2$ two-level atoms, each with its
individual frequency $\omega_{n}$ and raising (lowering) Pauli operators
$\sigma_{n}$ ($\sigma^{\dagger}_{n}$). Each atom couples with the strength
$g_{n}$ to a monochromatic field, i.e. to a quantized harmonic oscillator with
frequency $\omega_{0}$ and annihilation (creation) operator $b$
($b^{\dagger}$). When the rotating wave approximation is applied, i.e. when
$g_{n}\ll|\omega_{0}-\omega_{n}|\ll\omega_{n}$, $\forall n$, the counter-
rotating terms $\sigma_{n}b$ and $\sigma_{n}^{\dagger}b^{\dagger}$ can be
ignored to obtain the Jaynes-CummingsJaynes and Cummings (1963) (for $N=1$) or
Tavis-CummingsTavis and Cummings (1968) ($N>1$) model. In this work, since
$N=2$, we target the two-excitation manifold of the Tavis-Cummings
Hamiltonian, which is spanned by the basis states222The notation indicates the
states ($\ket{g}$ or $\ket{e}$) of the two two-level atoms, followed by the
number of photons left in the cavity mode. $\ket{gg;2}$, $\ket{ge;1}$,
$\ket{eg;1}$, and $\ket{ee;0}$, that contain, in total, two shared excitations
between field and atomic degrees of freedom. Diagonalization of the Tavis-
Cummings Hamiltonian in this manifold yields the four dressed states
$\ket{f_{1}},\ldots,\ket{f_{4}}$, where the two central states $\ket{f_{2}}$
and $\ket{f_{3}}$ are degenerate with energy $2\omega_{0}$; however, when
$g_{n}\sim|\omega_{0}-\omega_{n}|$, the counter-rotating terms of the full
Hamiltonian (17) lift this degeneracy. The discrimination of either state
against the other, then, depends on the competition between their energy
separation and their linewidths due to lifetime broadening. When these are
similar, we have the ideal scenario to test the ability of our method to
selectively excite just one of them.

The dressed electronic states of the atom-field Hamiltonian are also called
(cavity) _polaritons_ , since (17) constitutes a minimal model of the coupling
between photons and the oscillating electric dipoles in a loosely bound
crystal Mahan (2000). In this context, and in the text that follows, states
from the single- and double-excitation manifolds are addressed, respectively,
as “polaritons” and “bipolaritons”. Two-photon absorption to the manifold of
bipolariton states was recently considered in Gu and Mukamel, 2020. The
transition was excited by an entangled biphoton state, created by a cw pump
laser, with a frequency sum matching the excitation energy
$\omega_{1}+\omega_{2}$ of the targeted two-excitation manifold. It was shown
that when the bandwidth of the individual photons, which is determined by the
so-called entanglement time, becomes very narrow, such that the biphoton state
becomes effectively separable, certain bipolariton states become unaccessible
(“dark”Gu and Mukamel (2020)) due to destructive interference between
excitation pathways.

Here, instead, we are interested in the interference between different
excitation pathways that manifests in the general spectral structure of the
optimal states (15), with no further constraints beyond those introduced by
the functional (9). We therefore gather all spectral information on the
(dressed) matter degrees of freedom by diagonalizing the Hamiltonian (17) with
the maximal number of excitations in the cavity mode fixed at $15$ photons
(which is sufficient for numerical convergence). To ease comparison, we use
the parameter values extracted from Ref. Gu and Mukamel, 2020:
$\omega_{0}=$1\text{\,}\mathrm{e}\mathrm{V}$$, $\omega_{1}=0.8\omega_{0}$,
$\omega_{2}=1.2\omega_{0}$, $g_{1}=g_{2}=$0.14\text{\,}\mathrm{e}\mathrm{V}$$.
Of the polaritonic spectrum we consider the non-degenerate ground state
$\ket{g}$, the $n_{e}=3$ polaritons in the single-excitation manifold, and the
$n_{f}=4$ bipolaritons in the double-excitation manifold (note that, under the
above assumption of sufficiently weak coupling $g_{n}$, the uncoupled
polaritonic and bipolaritonic manifolds remain spectrally well-separated, such
that the attributes “(bi)polaritonic” still remain meaningful). As in Sec. II,
we set the origin of our energy scale at $\omega_{g}$, although, for
completeness, we report the unshifted energy levels of the dressed states in
Table 1. To characterize the available excitation pathways from the ground
state to the two-excitation manifold, we compute the dipole matrix elements
entering (II). The element
$\mu_{rs}=\sum_{n}\braket{r}{\sigma_{n}^{\dagger}+\sigma_{n}}{s}$ is taken
along the polarization direction of the electric field, and connects the
eigenstates $\ket{r}$ and $\ket{s}$ of (17) via the excitation of either two-
level system, which we take as aligned dipoles. The values of the dipole
matrix elements are also reported in Table 1.333Note that, as a peculiar
feature of the dressed state scenario, the matrix elements also include
transitions via sidebands, since
$\mu_{rs}=\braket{r}{D}{s}=\sum_{\alpha\alpha^{\prime}m}c^{(s)}_{\alpha
m}[c^{(r)}_{\alpha^{\prime}m}]^{*}\braket{\alpha^{\prime};m}{D}{\alpha;m},$
(18) where $D=\sum_{n}\sigma^{\dagger}_{n}+\sigma_{n}$ and we expand the
dressed states as $\ket{r}=\sum_{\alpha m}c^{(r)}_{\alpha m}\ket{\alpha;m}$,
with $\alpha,\alpha^{\prime}\in\\{gg,ge,eg,ee\\}$. For the case considered
here of small coupling $g_{n}$, however, these sideband transitions are orders
of magnitude less likely than the resonant transitions we address.

The last quantity we need to calculate the overlaps in (III.2) are the
linewidths of the (bi)polariton states. As we mentioned above, we wish to
assess how well our method can resolve two near-degenerate states in the two-
excitation manifold. For this reason, we set, for now,
$\gamma_{f}=0.01\omega_{0}\approx(\omega_{f_{3}}-\omega_{f_{2}})/3$. This
implies a significant overlap between the indistinctive optimal pulses (8)
that excite each bipolariton individually, as one finds directly with (III.2).
For the single-excitation manifold, instead, we take
$\gamma_{e}=\gamma_{f}/2$, as one would expect for the radiative decay of two
uncoupled two-level systems.

###  IV.2 Optimal states

The selective optimal two-photon wave functions (15) exciting either one of
the bipolariton states $\ket{f_{j}}$ are distinct. Their moduli square,
$|\tilde{\Phi}_{f_{j}}(\omega_{1},\omega_{2})|^{2}$, are all displayed in Fig.
1(a). If the excitation frequency of the target state is $\omega_{f}$, the
maxima of the wave functions are aligned along the antidiagonal
$\omega_{1}+\omega_{2}=\omega_{f}$ (black dotted line). Since the excitation
of $\ket{e_{2}}$ is suppressed ($\mu_{ge_{2}}\approx 0$, see table 1), any
bipolariton will be reached predominantly by either exciting $\ket{e_{1}}$ or
$\ket{e_{3}}$ first. Hence, we will get maximal population in $\ket{f_{1}}$,
for instance, with photon pairs of frequencies
$(\omega_{e_{1}},\omega_{f_{1}}-\omega_{e_{1}})$ and, to a lesser
degree444Since $\mu_{ge_{3}}\mu_{e_{3}f}<\mu_{ge_{1}}\mu_{e_{1}f}$.,
$(\omega_{e_{3}},\omega_{f_{1}}-\omega_{e_{3}})$. Consequently, these two
points in the two-photon frequency space is where the optimized wave
function’s density is peaked. As discussed at the end of Sec. II, the photons
couple identically to the matter, which means that either photon can excite
the $\ket{g}\rightarrow\ket{e}$ transition, while the other must complete the
two-photon absorption by exciting $\ket{e}\rightarrow\ket{f}$. As a
consequence, the optimal wave function is symmetric with respect to the two
frequencies $\omega_{1}$ and $\omega_{2}$.

The bipolariton population correspondingly excited by each such pulses is
given by the red (left) bars in Fig. 1(b). Note that the populations are given
in units of the state normalization $\mathcal{N}_{f}$ introduced in Eq. (14).
When columns are not visible, the populations are very small compared to the
dominating states, but are never exactly zero. From these histograms we see
that we can populate almost perfectly the bipolariton states $\ket{f_{1}}$ and
$\ket{f_{4}}$, which are well separated in energy from competing states. The
excitation targeting either $\ket{f_{2}}$ or $\ket{f_{3}}$, however, induces a
fraction of population also in the respective other state. This is consistent
with the two states having a non-negligible overlap for
$\gamma_{f}=0.01\omega_{0}$, as in this regime
$\omega_{f_{3}}-\omega_{f_{2}}\approx 3\gamma_{f}$.

As shown in Refs. Schlawin and Buchleitner, 2017; Carnio, Buchleitner, and
Schlawin, 2021, the optimal population of the bipolaritons is achieved by the
optimized coherent superposition of different frequency modes, such that the
atomic response to any frequency pair $(\omega_{1},\omega_{2})$ adds up
constructively. In general, this can only be accomplished if the incoming
field modes exhibit entanglement. The minimal set of modes required to
construct a given entangled state can be computed using the Schmidt
decompositionNielsen and Chuang (2009):

$\tilde{\Phi}(\omega_{1},\omega_{2})=\sum_{j=1}^{M}r_{j}\phi_{j}(\omega_{1})\psi_{j}(\omega_{2}),$
(19)

which is a weighted sum of $M$ orthonormal modes $\\{\phi_{j}\\}_{j}$ and
$\\{\psi_{j}\\}_{j}$. The weights $r_{j}$ can be chosen real and are listed in
decreasing order, by convention. The normalization of the state requires
$\sum_{j}r_{j}^{2}=1$.

The Schmidt modes are useful to construct the optimal _separable_ or
_classical_ state, which contains no quantum correlations:

$\tilde{\Phi}_{\text{cl}}(\omega_{1},\omega_{2})=\phi_{1}(\omega_{1})\psi_{1}(\omega_{2}).$
(20)

It can excite a fraction $r_{1}^{2}$ of the optimal population Schlawin and
Buchleitner (2017). This “classical” population, which we plot as blue (right)
columns in Fig. 1(b), can be excited in principle by conventionally shaping
the frequency spectrum of the individual broadband photons using, for
instance, a spatial light modulator to optimize the photonic pulse forms
Dudovich  _et al._ (2001); Dayan  _et al._ (2004); Brecht  _et al._ (2015). It
cannot however rely on entanglement, i.e. on the superposition of different
products of modes as in (19).

In Fig. 1(c) we show the optimal classical states derived from each optimal
wave function of row (a). The comparison between the two rows clarifies the
role of the coherent superposition of the modes in (19): only their
constructive or destructive interference can reproduce the profile, dominated
by the antidiagonal, of the optimal two-photon wave function. If we consider
again the case of $\ket{f_{1}}$, we see that the classical pulse excites the
$(\omega_{e_{1}},\omega_{f_{1}}-\omega_{e_{1}})$ and
$(\omega_{e_{3}},\omega_{f_{1}}-\omega_{e_{3}})$ transitions, which are
resonant with $\ket{f_{1}}$. In addition, however, there is also substantial
spectral weight at $(\omega_{e_{1}},\omega_{e_{1}})$ and
$(\omega_{f_{1}}-\omega_{e_{1}},\omega_{f_{1}}-\omega_{e_{1}})$, i.e. at
frequency combinations that are off-resonant for the given $\gamma_{f}$. In
the absence of quantum correlations, they cannot be suppressed. The
superposition of Schmidt modes in the entangled state enhances the first two
resonant frequency combinations and suppresses the off-resonant combinations,
thus enhancing the population by roughly a factor two compared to the
classical case [see Fig. 1(_b_)].

Since the natural linewidth $\gamma_{f}$ of the bipolaritons determines how
close to resonance the two-photon transitions are, Fig. 2 shows the optimal
two-photon wave functions, the target state population histograms and the
accordingly shaped classical pulses for bipolaritons with a linewidth enlarged
by a factor ten (i.e. shorter lifetimes),
$\gamma^{\prime}_{f}=0.1\omega_{0}=10\gamma_{f}$. As one would expect
intuitively, while the general structure of the optimal states can still be
recognized, given our knowledge of their structures for narrow linewidths as
depicted in Fig. 1, some details are washed out by the broadened resonances.
For the state $\ket{f_{1}}$, for instance, the square-like structure, visible
in Fig. 1, of the classical pulse, for $\gamma_{f}$, cannot be resolved
anymore for $\gamma^{\prime}_{f}$, because the single-photon frequency
distributions are much broader than the difference
$(\omega_{f_{1}}-\omega_{e_{1}})-\omega_{e_{1}}=\omega_{f_{1}}-2\omega_{e_{1}}$.
This explains why the population induced by the shaped classical state is
larger for $\gamma^{\prime}_{f}$ than for $\gamma_{f}$. Similarly, when we try
to optimally populate $\ket{f_{2}}$, we now obtain also a larger population in
$\ket{f_{3}}$, as well as some in $\ket{f_{4}}$. The latter is due to the
classical pulse having a significant peak around
$(\omega_{e_{3}},\omega_{e_{3}})$, which is broad enough to excite the
$(\omega_{e_{3}},\omega_{f_{4}}-\omega_{e_{3}})$ two-photon transition to
$\ket{f_{4}}$.

Figure 1: (a) Selective optimal two-photon wave functions (15) exciting, from
left to right, the bipolaritons $\ket{f_{1}}$ to $\ket{f_{4}}$, as indicated
in the top-right corner. Colors and fine-dotted contour lines indicate the
intensity of $|\tilde{\Phi}(\omega_{1},\omega_{2})|^{2}$ on a logarithmic
scale. Gray dashed lines indicate the resonance frequency of the states
$e_{j}$ in the single-excitation manifold, while black dotted lines delineate
the antidiagonal $\omega_{1}+\omega_{2}=\omega_{f}$. (b) Red (left) and blue
(right) bars: populations in each bipolariton, in units of $\mathcal{N}_{f}$
(14), excited by, respectively, the selective optimal state (15) depicted in
the panel directly above, and the corresponding classical pulse (20) depicted
in the panel directly below. The bars that are not visible are too small
compared to the scale, but are never exactly zero. (c) Classical pulses (20)
obtained from the pulses in row (a). Colors and fine-dotted contour lines
indicate the intensity of
$|\tilde{\Phi}_{\text{cl}}(\omega_{1},\omega_{2})|^{2}$ on a logarithmic
scale. Gray dashed lines indicate the resonance frequency of the states
$e_{j}$ in the single-excitation manifold, while black dotted lines delineate
the antidiagonal $\omega_{1}+\omega_{2}=\omega_{f}$. Simulations run with
$\gamma_{f}=0.01\omega_{0}$.  Figure 2: As in Fig. 1: (a) Selective optimal
two-photon wave functions (15) exciting, from left to right, the bipolaritons
$\ket{f_{1}}$ to $\ket{f_{4}}$, as indicated in the top-right corner. Colors
and fine-dotted contour lines indicate the intensity of
$|\tilde{\Phi}(\omega_{1},\omega_{2})|^{2}$ on a logarithmic scale. Gray
dashed lines indicate the resonance frequency of the states $e_{j}$ in the
single-excitation manifold, while black dotted lines delineate the
antidiagonal $\omega_{1}+\omega_{2}=\omega_{f}$. (b) Red (left) and blue
(right) bars: populations in each bipolariton, in units of $\mathcal{N}_{f}$
(14), excited by, respectively, the selective optimal state (15) depicted in
the panel directly above, and the corresponding classical pulse (20) depicted
in the panel directly below. The bars that are not visible are too small
compared to the scale, but are never exactly zero. (c) Classical pulses (20)
obtained from the pulses in row (a). Colors and fine-dotted contour lines
indicate the intensity of
$|\tilde{\Phi}_{\text{cl}}(\omega_{1},\omega_{2})|^{2}$ on a logarithmic
scale. Gray dashed lines indicate the resonance frequency of the states
$e_{j}$ in the single-excitation manifold, while black dotted lines delineate
the antidiagonal $\omega_{1}+\omega_{2}=\omega_{f}$. Simulations run with
$\gamma_{f}=0.1\omega_{0}$.

###  IV.3 Selectivity

The previous subsection showed that, if the linewidth of the states is much
larger than the frequency differences involved in the two-photon transitions,
(i) the quantum advantage due to the frequency entanglement in the optimal
pulses is reduced, (ii) the achievable populations in $\ket{f_{1}}$ and
$\ket{f_{4}}$ decrease (increase) with increasing $\gamma_{f}$ when using an
entangled (classical) two-photon state, and (iii) the target population in
$\ket{f_{2}}$ ($\ket{f_{3}}$) decreases with increasing $\gamma_{f}$, while
that of $\ket{f_{3}}$ ($\ket{f_{2}}$) increases. This latter point, in
particular, implies that the _selectivity_ of the optimal pulse (15), captured
by the contrast between the populations of $\ket{f_{2}}$ and $\ket{f_{3}}$,

$S=\frac{|p_{f_{2}}-p_{f_{3}}|}{p_{f_{2}}+p_{f_{3}}},$ (21)

worsens for broader linewidths.

Finally, an analysis of how the desired selectivity scales with the
intermediate and target states’ linewidths ($\gamma_{f}=2\gamma_{e}$), for
optimal quantum vs classical two-photon states, reveals the advantage of
performing the optimization presented in this work. In Fig. 3 we compare the
selectivity in those cases where the population of $\ket{f_{2}}$ is induced by
selective optimal entangled states (15) to that achieved by shaped classical
states of light (20). In addition, we compare these results to those achieved
with the indistinctive optimal quantum state (8), which, we remind from Sec.
III, is optimized to excite the $\ket{g}\rightarrow\ket{f_{2}}$ transition
most efficiently, without the additional constraint to minimize the excitation
of other nearby bipolaritons from the same manifold. For all three possible
injected two-photon states, the selectivity decays with $\gamma_{f}$, in a
qualitatively exponential fashion (top panel). Yet, _selective_ optimization
as conceived in Sec. III, for classical as well as for quantum light, slows
this decay down. For instance, at $\gamma_{f}=0.05\omega_{0}$, the
indistinctive optimal (quantum) pulse will hardly discriminate $\ket{f_{2}}$
against $\ket{f_{3}}$ any more, while the optimization procedure can still
achieve ca. $25$ % selectivity. This relative improvement is visualized by
plotting the ratio of selective vs. indistinctive yields in the figure’s
bottom panel: this ratio actually increases with increasing $\gamma_{f}$. This
has a very intuitive explanation: the possibility of constructing suitable
linear superpositions (15) of indistinctive optimal states (8) relies, in the
first place, on their overlap, which increases with $\gamma_{f}$.

Figure 3: Top: selectivity $S$ from Eq. (21) calculated for selectively
optimal [(15), red circles], shaped classical [(20), blue diamonds] and
indistinctive optimal pulses [(8), green squares]. Bottom: ratio of the
selectivity achieved with selectively optimal (red circles) or shaped
classical pulses (blue diamonds), to the selectivity achievable with an
indistinctive optimal pulse.

##  V Conclusions

We have investigated continuous-mode two-photon states to populate a target
matter state reachable by two-photon absorption. The method discussed in this
work optimizes the excitation’s selectivity, i.e. it maximizes the target
state population while minimizing residual population within the complement of
the target state’s manifold. The optimal quantum state of light to drive this
transition can be obtained by solving a generalized eigenvalue problem where
all matrices depend analytically on the spectral properties of the driven
system, and have a dimension given by that of the manifold the target state is
embedded in.

We have applied our method to the specific setting of two non-interacting
atoms dressed by a cavity mode, with the goal of driving the transition from
the ground state to one of the four bipolariton states in the two-excitation
manifold. If the bipolariton states are well separated in energy, our
procedure is equivalent to the optimization of a single pathway analyzed in
previous workSchlawin and Buchleitner (2017); Carnio, Buchleitner, and
Schlawin (2021). If they overlap, instead, the excitation of an individual
transition inevitably induces transitions to the nearby states, too. In this
case, however, we manage to obtain an appreciable selectivity even when it
would be impossible to otherwise discern closely neighboring resonances. The
Hamiltonian (17) considered in this study is formally similar to excitonic
models of molecular aggregatesMay and Kühn (2011) in physical chemistry, in
the sense that, in the parameter regime considered here, we obtain
energetically well-separated manifolds of states, where adjacent ones are
dipole-coupled. It seems highly likely that the excitation physics presented
here will carry over to entangled two-photon excitation of, e.g., molecular
aggregates Schlawin  _et al._ (2013); Raymer  _et al._ (2013). We note,
however, that the inevitable coupling of the electronic states to
environmental degrees of freedom will induce additional relaxation processes
such as incoherent electronic population transfer. These processes are not
captured by the non-Hermitian Hamiltonian considered here, and would rather
require an open systems description, where the material degrees of freedom
evolve according to a master equation.

While the solution of the selective optimization problem is general, its
benchmark inevitably depends on the specific structure of the target spectrum
considered. In particular, to enable a comparison with Gu and Mukamel, 2020,
we considered a sufficiently small coupling constant $g_{n}$ in (17), such
that the single- and double-excitation manifolds remain well separated in
energy. The selective excitation of transitions in a spectrum where said
manifolds mix, instead, will be a topic of future research.

###### Acknowledgements.

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.

## DATA AVAILABILITY

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

## References

  * Mukamel (1995) S. Mukamel,  _Principles of Nonlinear Optical Spectroscopy_ , Optical and Imaging Sciences Series (Oxford University Press, 1995).
  * Hamm and Zanni (2011) P. Hamm and M. Zanni,  _Concepts and Methods of 2D Infrared Spectroscopy_ (Cambridge University Press, 2011).
  * Boyd and Boyd (1992) R. Boyd and M. Boyd,  _Nonlinear Optics_ (Academic Press, 1992).
  * Dorfman, Schlawin, and Mukamel (2016) K. E. Dorfman, F. Schlawin, and S. Mukamel, Rev. Mod. Phys. 88, 045008 (2016).
  * Schlawin, Dorfman, and Mukamel (2018) F. Schlawin, K. E. Dorfman, and S. Mukamel, Acc. Chem. Res. 51, 2207 (2018).
  * Gilaberte Basset  _et al._ (2019) M. Gilaberte Basset, F. Setzpfandt, F. Steinlechner, E. Beckert, T. Pertsch, and M. Graefe, Laser Photonics Rev. 13, 1900097 (2019).
  * Mukamel  _et al._ (2020) S. Mukamel, M. Freyberger, W. Schleich, M. Bellini, A. Zavatta, G. Leuchs, C. Silberhorn, R. W. Boyd, L. L. Sánchez-Soto, A. Stefanov, M. Barbieri, A. Paterova, L. Krivitsky, S. Shwartz, K. Tamasaku, K. Dorfman, F. Schlawin, V. Sandoghdar, M. Raymer, A. Marcus, O. Varnavski, T. Goodson, Z.-Y. Zhou, B.-S. Shi, S. Asban, M. Scully, G. Agarwal, T. Peng, A. V. Sokolov, Z.-D. Zhang, M. S. Zubairy, I. A. Vartanyants, E. del Valle, and F. Laussy, J. Phys. B At. Mol. Opt. Phys. 53, 072002 (2020).
  * Szoke  _et al._ (2020) S. Szoke, H. Liu, B. P. Hickam, M. He, and S. K. Cushing, J. Mater. Chem. C 8, 10732 (2020).
  * del Valle  _et al._ (2012) E. del Valle, A. Gonzalez-Tudela, F. P. Laussy, C. Tejedor, and M. J. Hartmann, Phys. Rev. Lett. 109, 183601 (2012).
  * Gonzalez-Tudela  _et al._ (2013) A. Gonzalez-Tudela, F. P. Laussy, C. Tejedor, M. J. Hartmann, and E. del Valle, New J. Phys. 15, 033036 (2013).
  * Wientjes  _et al._ (2014) E. Wientjes, J. Renger, A. G. Curto, R. Cogdell, and N. F. van Hulst, Nat. Commun. 5, 4236 (2014).
  * Holdaway, Notararigo, and Olaya-Castro (2018) D. I. H. Holdaway, V. Notararigo, and A. Olaya-Castro, Phys. Rev. A 98, 063828 (2018).
  * Sánchez Muñoz and Schlawin (2020) C. Sánchez Muñoz and F. Schlawin, Phys. Rev. Lett. 124, 203601 (2020).
  * Krüger  _et al._ (2010) T. P. Krüger, V. I. Novoderezhkin, C. Ilioaia, and R. van Grondelle, Biophys. J. 98, 3093 (2010).
  * Hildner  _et al._ (2013) R. Hildner, D. Brinks, J. B. Nieder, R. J. Cogdell, and N. F. van Hulst, Science 340, 1448 (2013).
  * Kyeyune  _et al._ (2019) F. Kyeyune, J. L. Botha, B. van Heerden, P. Malỳ, R. van Grondelle, M. Diale, and T. P. J. Krüger, Nanoscale 11, 15139 (2019).
  * Dorfman and Mukamel (2012) K. E. Dorfman and S. Mukamel, Phys. Rev. A 86, 013810 (2012).
  * Schlawin, Dorfman, and Mukamel (2016) F. Schlawin, K. E. Dorfman, and S. Mukamel, Phys. Rev. A 93, 023807 (2016).
  * Dorfman and Mukamel (2018) K. E. Dorfman and S. Mukamel, Proc. Natl. Acad. Sci. U.S.A. 115, 1451 (2018).
  * Zhang  _et al._ (2018) Z. Zhang, P. Saurabh, K. E. Dorfman, A. Debnath, and S. Mukamel, J. Chem. Phys. 148, 074302 (2018).
  * Shi, Kumar, and Lee (2017) S. Shi, P. Kumar, and K. Lee, Nat. Commun. 8, 1934 (2017).
  * Li  _et al._ (2020a) F. Li, T. Li, M. O. Scully, and G. S. Agarwal, (2020a), arXiv:2007.05583 [physics.optics] .
  * Li  _et al._ (2020b) T. Li, F. Li, C. Altuzarra, A. Classen, and G. S. Agarwal, Appl. Phys. Lett. 116, 254001 (2020b).
  * Georgiades  _et al._ (1995) N. P. Georgiades, E. S. Polzik, K. Edamatsu, H. J. Kimble, and A. S. Parkins, Phys. Rev. Lett. 75, 3426 (1995).
  * Dayan  _et al._ (2004) B. Dayan, A. Pe’er, A. A. Friesem, and Y. Silberberg, Phys. Rev. Lett. 93, 023005 (2004).
  * Lee and Goodson (2006) D.-I. Lee and T. Goodson, J. Phys. Chem. B 110, 25582 (2006).
  * Upton  _et al._ (2013) L. Upton, M. Harpham, O. Suzer, M. Richter, S. Mukamel, and T. Goodson, J. Phys. Chem. Lett. 4, 2046 (2013).
  * Parzuchowski  _et al._ (2020) K. M. Parzuchowski, A. Mikhaylov, M. D. Mazurek, R. N. Wilson, D. J. Lum, T. Gerrits, C. H. C. J. au2, M. J. Stevens, and R. Jimenez, (2020), arXiv:2008.02664 [quant-ph] .
  * Landes  _et al._ (2020) T. Landes, M. Allgaier, S. Merkouche, B. J. Smith, A. H. Marcus, and M. G. Raymer, (2020), arXiv:2012.06736 [quant-ph] .
  * Tabakaev  _et al._ (2021) D. Tabakaev, M. Montagnese, G. Haack, L. Bonacina, J.-P. Wolf, H. Zbinden, and R. T. Thew, Physical Review A 103 (2021), 10.1103/physreva.103.033701.
  * Raymer  _et al._ (2020) M. G. Raymer, T. Landes, M. Allgaier, S. Merkouche, B. J. Smith, and A. H. Marcus, (2020), arXiv:2012.05375 [quant-ph] .
  * Raymer, Landes, and Marcus (2021) M. G. Raymer, T. Landes, and A. H. Marcus, (2021), arXiv:2103.02551 [quant-ph] .
  * Schlawin  _et al._ (2013) F. Schlawin, K. E. Dorfman, B. P. Fingerhut, and S. Mukamel, Nat. Commun. 4, 1782 (2013).
  * Raymer  _et al._ (2013) M. G. Raymer, A. H. Marcus, J. R. Widom, and D. L. Vitullo, J. Phys. Chem. B 117, 15559 (2013).
  * León-Montiel  _et al._ (2019) R. d. J. León-Montiel, J. c. v. Svozilík, J. P. Torres, and A. B. U’Ren, Phys. Rev. Lett. 123, 023601 (2019).
  * Ishizaki (2020) A. Ishizaki, J. Chem. Phys. 153, 051102 (2020).
  * Fujihashi, Shimizu, and Ishizaki (2020) Y. Fujihashi, R. Shimizu, and A. Ishizaki, Phys. Rev. Research 2, 023256 (2020).
  * Li  _et al._ (2017) H. Li, A. Piryatinski, J. Jerke, A. R. S. Kandada, C. Silva, and E. R. Bittner, Quantum Sci. Technol. 3, 015003 (2017).
  * Li  _et al._ (2019) H. Li, A. Piryatinski, A. R. Srimath Kandada, C. Silva, and E. R. Bittner, J. Chem. Phys. 150, 184106 (2019).
  * Castro, Appel, and Rubio (2019) A. Castro, H. Appel, and A. Rubio, The European Physical Journal B 92, 223 (2019).
  * Csehi  _et al._ (2019) A. Csehi, A. Vibók, G. J. Halász, and M. Kowalewski, Phys. Rev. A 100, 053421 (2019).
  * Oka (2018) H. Oka, Phys. Rev. A 97, 033814 (2018).
  * Rabitz  _et al._ (2000) H. Rabitz, R. de Vivie-Riedle, M. Motzkus, and K. Kompa, Science 288, 824 (2000).
  * Brif, Chakrabarti, and Rabitz (2010) C. Brif, R. Chakrabarti, and H. Rabitz, New J. Phys. 12, 075008 (2010).
  * Dubov and Rabitz (1996) V. Dubov and H. Rabitz, Phys. Rev. A 54, 710 (1996).
  * Assion (1998) A. Assion, Science 282, 919 (1998).
  * Meshulach and Silberberg (1998) D. Meshulach and Y. Silberberg, Nature 396, 239 (1998).
  * Meshulach and Silberberg (1999) D. Meshulach and Y. Silberberg, Phys. Rev. A 60, 1287 (1999).
  * Dudovich  _et al._ (2001) N. Dudovich, B. Dayan, S. M. Gallagher Faeder, and Y. Silberberg, Phys. Rev. Lett. 86, 47 (2001).
  * Shapiro and Brumer (2012) M. Shapiro and P. Brumer,  _Quantum Control of Molecular Processes_ (Wiley, Hoboken, USA, 2012).
  * Schlawin and Buchleitner (2017) F. Schlawin and A. Buchleitner, New J. Phys. 19, 013009 (2017).
  * Carnio, Buchleitner, and Schlawin (2021) E. G. Carnio, A. Buchleitner, and F. Schlawin, (2021), arXiv:2105.13876 [quant-ph] .
  * Gu and Mukamel (2020) B. Gu and S. Mukamel, J. Phys. Chem. Lett. 11, 8177 (2020).
  * Loudon (2000) R. Loudon,  _The quantum theory of light_ , 3rd ed. (Oxford University Press, Oxford, 2000).
  * Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Note (1) We remark here that, if we restrict our optimization problem to two-photon states, with one photon per mode, as we do here, the normalization constraint $\mathinner{\langle{\Phi|n_{1}n_{2}|\Phi}\rangle}$ used in Schlawin2017a, Carnio21, where $n_{j}$ is the photon number operator for beam $j$, corresponds to $\mathinner{\langle{\Phi|\Phi}\rangle}$ used here.
  * Bronshtein  _et al._ (2007) I. N. Bronshtein, K. A. Semendyayev, G. Musiol, and H. Mühlig,  _Handbook of Mathematics_ , 5th ed. (Springer Berlin Heidelberg, 2007).
  * Cohen-Tannoudji, Dupont-Roc, and Grynberg (1998) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Atom-Photon Interactions_ (Wiley, 1998).
  * Jaynes and Cummings (1963) E. Jaynes and F. Cummings, Proceedings of the IEEE 51, 89 (1963).
  * Tavis and Cummings (1968) M. Tavis and F. W. Cummings, Physical Review 170, 379 (1968).
  * Note (2) The notation indicates the states ($\mathinner{|{g}\rangle}$ or $\mathinner{|{e}\rangle}$) of the two two-level atoms, followed by the number of photons left in the cavity mode.
  * Mahan (2000) G. D. Mahan,  _Many-Particle Physics_, 3rd ed. (Springer Boston, 2000).
  * Note (3) Note that, as a peculiar feature of the dressed state scenario, the matrix elements also include transitions via sidebands, since  $\mu_{rs}=\mathinner{\langle{r|D|s}\rangle}=\sum@\slimits@_{\alpha\alpha^{\prime}m}c^{(s)}_{\alpha m}[c^{(r)}_{\alpha^{\prime}m}]^{*}\mathinner{\langle{\alpha^{\prime};m|D|\alpha;m}\rangle},$ (22) where $D=\sum@\slimits@_{n}\sigma^{\dagger}_{n}+\sigma_{n}$ and we expand the dressed states as $\mathinner{|{r}\rangle}=\sum@\slimits@_{\alpha m}c^{(r)}_{\alpha m}\mathinner{|{\alpha;m}\rangle}$, with $\alpha,\alpha^{\prime}\in\\{gg,ge,eg,ee\\}$. For the case considered here of small coupling $g_{n}$, however, these sideband transitions are orders of magnitude less likely than the resonant transitions we address.
  * Note (4) Since $\mu_{ge_{3}}\mu_{e_{3}f}<\mu_{ge_{1}}\mu_{e_{1}f}$.
  * Nielsen and Chuang (2009) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2009).
  * Brecht  _et al._ (2015) B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, Phys. Rev. X 5, 041017 (2015).
  * May and Kühn (2011) V. May and O. Kühn,  _Charge and Energy Transfer Dynamics in Molecular Systems_, 3rd ed. (John Wiley & Sons, Ltd, 2011).