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.^{1,2} In most situations, this approximation is extremely well justified owing to the weak nonlinearity of light–matter interactions in free space.^{3} Nevertheless, recent years have seen the rapid rise in theoretical investigations as well as first proof-of-concept experiments that challenge this convention and investigation on how quantum properties of light can be applied or exploited beneficially in spectroscopic applications.^{4–8} This includes the use of photon correlation measurements to analyze the light fields emitted by a sample^{9–13} in single-molecule spectroscopy^{14–16} or to exploit coincidence measurements to detect particular spectral features in the sample^{17–20} as well as the generation of photonic entanglement in fluorescent proteins.^{21}

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.^{22,23} Nonlinear optical signals, such as two-photon absorption, scale linearly with the photon flux,^{24–26} which could enable nonlinear spectroscopy of photosensitive samples at very low intensities. Following earlier experiments on two-photon absorption in biomolecules,^{26,27} a series of recent experiments have scrutinized the situation^{28–30} 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.^{31,32} 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^{33–35} and enable ultrafast spectroscopy in a cw setup.^{36} They could also be used for the generation of pseudo-sunlight to imitate natural conditions^{37} and provide sensitive probes for dynamical symmetry breaking^{38} and many-body correlations.^{39} The control of quantum correlations using temperature,^{35} photon statistics,^{40,41} or spectral shaping^{42} 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 optimized (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.^{43,44} In particular, optimal control of two-photon absorption was described and implemented experimentally in the late 1990s.^{45–49} 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. 50. However, we recently showed that this is not true for weak broadband fields, where quantum correlations of light can enhance excitation probabilities.^{51,52}

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. 53. 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*_{f} and *H*_{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*_{f} + *H*_{m} + *W*.

Both quantized pulses are injected along fixed and distinct spatial directions and are described by continuous-mode operators, *a*_{1}(*ω*) and *a*_{2}(*ω*), respectively.^{54} The fields’ Hamiltonian, ignoring the vacuum energy, is then $Hf=\u2211j\u222b0\u221ed\omega \u210f\omega aj\u2020(\omega )aj(\omega )$. The positive-frequency electric field operator acting on the Hilbert space of photon *j* (in the interaction picture with respect to *H*_{f}) reads

Here, *z* is the position along the propagation direction, *A* is the quantization area perpendicular to it, and *c* is 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 Δ*ω*, distributed around a central frequency *ω*_{0} ≫ Δ*ω*. 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*^{54} as

where $E0=(\u210f\omega 0/4\pi \u03f50cA)1/2$ approximates the field normalization of (1). The range of integration (−∞, ∞) will be assumed from now on.

As for the matter degrees of freedom, we consider a system with a ground state |*g*⟩, multiple intermediate states |*e*_{1}⟩, |*e*_{2}⟩, …, and multiple final states |*f*_{1}⟩, |*f*_{2}⟩, …. With each state |*s*⟩ we associate an energy *ℏω*_{s} (we set *ℏω*_{g} = 0), an inverse lifetime *γ*_{s} (*γ*_{g} = 0 since the ground state cannot decay), and a Lorentzian line shape $Ls(\omega )=(\omega \u2212zs)\u22121$, where *z*_{s} = *ω*_{s} − i*γ*_{s}. For eigenstates with a finite lifetime, we consider an effective, non-Hermitian^{55} matter Hamiltonian $Hm=\u2211j\u210fzej|ej\u3009\u3008ej|+\u2211k\u210fzfk|fk\u3009\u3008fk|$. Adjacent manifolds are dipole-coupled with dipole matrix elements (along the fields’ polarization) $\mu gej$ and $\mu ejfk$, respectively.

The light–matter interaction Hamiltonian in the rotating wave approximation—which is certainly valid in the present near-resonant perturbative regime where $\mu E0\u226a\u210f\omega 0$ (for any of the dipole matrix elements *μ* above)—and, as the subscript indicates, in the interaction picture with respect to *H*_{f} + *H*_{m}, is therefore given by^{54}

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

and $E+(\u2212)=E1+(\u2212)+E2+(\u2212)$ since we assume that each pulse couples identically to the matter degrees of freedom.

## III. OPTIMIZATION PROCEDURE

### A. Single final state

We seek a two-photon state |Φ_{f}⟩ of the incoming fields, which optimizes the two-photon transition from the ground state |*g*⟩ to a single final state |*f*⟩, via a manifold of *n*_{e} intermediate states |*e*_{j}⟩.

Mathematically, we consider an initial state |Φ_{f}⟩ ⊗ |*g*⟩ in the combined matter-field Hilbert space $H=Hf\u2297Hm$ (mirroring the subscripts of the Hamiltonians in Sec. II) and determine |Φ_{f}⟩ such that the time evolution operator *U*_{I}(*t*, *t*_{0}) in the interaction picture of *H*_{f} + *H*_{m} maximizes the population of the target state |0⟩ ⊗ |*f*⟩, where both photons have been absorbed in order to drive the matter degrees of freedom to the final state |*f*⟩. As described in detail in Ref. 52, we can find this transition amplitude perturbatively. Tracing out the matter degrees of freedom, the desired transition is mediated by an operator acting on $Hf$ alone,

with the explicit expression of the matter response function

The matter response function (6) is to be derived for the excitation of the population of |*f*⟩ at a fixed time *t*, with the interaction *W*_{I}(*t*) turned on at *t*_{0} → −∞. Without loss of generality, we can therefore set *t* = 0. The state |Φ_{f}⟩ can now be found by the variation of the functional^{51,52}

where $pf=|\u30080|T\u0302f|\Phi \u3009|2$ is the population in the target state |*f*⟩ created by the absorption of |Φ⟩, leaving the field in the vacuum state |0⟩. The Lagrange multiplier *λ* constrains the optimization to normalized states.^{56} Defining the (unnormalized) two-photon state $|Tf\u3009=T\u0302f\u2020|0\u3009$ via (5), the functional (7) is maximized^{51,52} by the state $|\Phi f\u3009=Nf\u22121/2|Tf\u3009$ whose associated *wave function* reads

given the appropriate normalization $Nf=\u27e8Tf|Tf\u27e9$, which we determine analytically in Subsection III B.

### B. 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 |*f*_{1}⟩ while minimizing the population of all other energetically near-degenerate states |*f*_{2}⟩, |*f*_{3}⟩, …, 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 functional (7) such that the target state’s population $pf1$ is maximized, while the excitation of any other states |*f*_{j≠1}⟩ is penalized,

Note that enforcing strictly vanishing populations in the states |*f*_{j}⟩ ≠ |*f*_{1}⟩ is prevented by the fact that the matter response functions (6) for different target states |*f*_{j}⟩, in general, exhibit finite overlap.

In analogy to Subsection III A, to maximize (9), we first define two-photon states $|Tfj\u3009=T\u0302fj\u2020|0\u3009$. We remind here that these states are not normalized, but their normalization will be taken care of in (12). By writing $pfj=|Tfj|\Phi |2$, functional (9) transforms into

The state $|\Phi \u0303\u3009$ that maximizes this functional is found by requiring the variational derivative with respect to the dual state^{52} ⟨Φ| to vanish,

with $Kfj=|Tfj\u3009\u3008Tfj|$. 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 ${|Tfj\u3009}j$. We then obtain the generalized eigenvalue problem^{57} (GEP)

where the matrix *M* is given by the overlaps

and the normalization constants

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

The *selective* optimal two-photon wave function $\Phi \u0303(\omega 1,\omega 2)$ is determined by the eigenvector $v\u0303$ 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 fj$ from (8),

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

where we remind the reader that the eigenvectors of a GEP are orthonormal with respect to the scalar product induced by the *M* matrix: $(vj,vk)M=vj*\u22c5M\u22c5vk=\delta jk$.

ω_{r}/ω_{0}
. | −2.04 × 10^{−2}
. | 0.698 . | 0.979 . | 1.26 . | 1.58 . | 1.97 . | 2.00 . | 2.37 . |
---|---|---|---|---|---|---|---|---|

μ_{rs}
. | g
. | e_{1}
. | e_{2}
. | e_{3}
. | f_{1}
. | f_{2}
. | f_{3}
. | f_{4}
. |

G | 0 | 1.09 | 1.61 × 10^{−3} | 0.939 | 0 | 0 | 0 | 0 |

e_{1} | 1.09 | 0 | 0 | 0 | 0.891 | 0.475 | 0.705 | 0.125 |

e_{2} | 1.61 × 10^{−3} | 0 | 0 | 0 | 0.757 | −7.95 × 10^{−2} | −7.34 × 10^{−2} | 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 × 10^{−2} | 0.515 | 0 | 0 | 0 | 0 |

f_{3} | 0 | 0.705 | −7.34 × 10^{−2} | 0.706 | 0 | 0 | 0 | 0 |

f_{4} | 0 | 0.125 | 0.688 | −0.779 | 0 | 0 | 0 | 0 |

ω_{r}/ω_{0}
. | −2.04 × 10^{−2}
. | 0.698 . | 0.979 . | 1.26 . | 1.58 . | 1.97 . | 2.00 . | 2.37 . |
---|---|---|---|---|---|---|---|---|

μ_{rs}
. | g
. | e_{1}
. | e_{2}
. | e_{3}
. | f_{1}
. | f_{2}
. | f_{3}
. | f_{4}
. |

G | 0 | 1.09 | 1.61 × 10^{−3} | 0.939 | 0 | 0 | 0 | 0 |

e_{1} | 1.09 | 0 | 0 | 0 | 0.891 | 0.475 | 0.705 | 0.125 |

e_{2} | 1.61 × 10^{−3} | 0 | 0 | 0 | 0.757 | −7.95 × 10^{−2} | −7.34 × 10^{−2} | 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 × 10^{−2} | 0.515 | 0 | 0 | 0 | 0 |

f_{3} | 0 | 0.705 | −7.34 × 10^{−2} | 0.706 | 0 | 0 | 0 | 0 |

f_{4} | 0 | 0.125 | 0.688 | −0.779 | 0 | 0 | 0 | 0 |

## 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 $pf1$. However, this information must be encoded in the coefficients $v\u0303j$ 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 mode.^{58} After defining the model, we study how the structure of the optimal two-photon states reflects the different *excitation pathways* |*g*⟩ → |*e*_{j}⟩ → |*f*_{k}⟩, i.e., the possible transitions, through the intermediate manifold of |*e*_{j}⟩ states, which maximize the yield of the target two-photon transition |*g*⟩ → |*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.

### A. Dressed-state Hamiltonian

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

with *ℏ* ≡ 1. We consider *N* = 2 two-level atoms, each with its individual frequency *ω*_{n} and raising (lowering) Pauli operators *σ*_{n} ($\sigma n\u2020$). Each atom couples with the strength *g*_{n} to a monochromatic field, i.e., to a quantized harmonic oscillator with frequency *ω*_{0} and annihilation (creation) operator *b* (*b*^{†}). When the rotating wave approximation is applied, i.e., when *g*_{n} ≪ |*ω*_{0} − *ω*_{n}| ≪ *ω*_{n}, ∀*n*, the counter-rotating terms *σ*_{n}*b* and $\sigma n\u2020b\u2020$ can be ignored to obtain the Jaynes–Cummings^{59} (for *N* = 1) or Tavis–Cummings^{60} (*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 states^{61} |*gg*; 2⟩, |*ge*; 1⟩, |*eg*; 1⟩, and |*ee*; 0⟩, which 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 |*f*_{1}⟩, …, |*f*_{4}⟩, where the two central states |*f*_{2}⟩ and |*f*_{3}⟩ are degenerate with energy 2*ω*_{0}; however, when *g*_{n} ∼ |*ω*_{0} − *ω*_{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.^{62} 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 Ref. 53. The transition was excited by an entangled biphoton state, created by a cw pump laser, with a frequency sum matching the excitation energy *ω*_{1} + *ω*_{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”^{53}) 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 functional (9). We therefore gather all spectral information on the (dressed) matter degrees of freedom by diagonalizing 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. 53: *ω*_{0} = 1 eV, *ω*_{1} = 0.8*ω*_{0}, *ω*_{2} = 1.2*ω*_{0}, and *g*_{1} = *g*_{2} = 0.14 eV. Of the polaritonic spectrum, we consider the non-degenerate ground state |*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 *ω*_{g}, although, for completeness, we report the unshifted energy levels of the dressed states in Table I. To characterize the available excitation pathways from the ground state to the two-excitation manifold, we compute the dipole matrix elements entering (4). The element $\mu rs=\u2211n\u3008r|\sigma n\u2020+\sigma n|s\u3009$ is taken along the polarization direction of the electric field and connects the eigenstates |*r*⟩ and |*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 I.^{63}

The last quantity we need to calculate the overlaps in (13) 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\u2248(\omega f3\u2212\omega f2)/3$. This implies a significant overlap between the indistinctive optimal pulses (8) that excite each bipolariton individually, as one finds directly with (13). For the single-excitation manifold, instead, we take *γ*_{e} = *γ*_{f}/2 as one would expect for the radiative decay of two uncoupled two-level systems.

### B. Optimal states

The selective optimal two-photon wave functions (15) exciting either one of the bipolariton states |*f*_{j}⟩ are distinct. Their moduli square, $|\Phi \u0303fj(\omega 1,\omega 2)|2$, are all displayed in Fig. 1(a). If the excitation frequency of the target state is *ω*_{f}, the maxima of the wave functions are aligned along the antidiagonal *ω*_{1} + *ω*_{2} = *ω*_{f} (black dotted line). Since the excitation of |*e*_{2}⟩ is suppressed ($\mu ge2\u22480$; see Table I), any bipolariton will be reached predominantly by either exciting |*e*_{1}⟩ or |*e*_{3}⟩ first. Hence, we will get maximal population in |*f*_{1}⟩, for instance, with photon pairs of frequencies $(\omega e1,\omega f1\u2212\omega e1)$ and, to a lesser degree,^{64} $(\omega e3,\omega f1\u2212\omega e3)$. 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 |*g*⟩ → |*e*⟩ transition, while the other must complete the two-photon absorption by exciting |*e*⟩ → |*f*⟩. As a consequence, the optimal wave function is symmetric with respect to the two frequencies *ω*_{1} and *ω*_{2}.

The bipolariton population correspondingly excited by each such pulse is given by the red (left) bars in Fig. 1(b). Note that the populations are given in units of the state normalization $Nf$ 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 |*f*_{1}⟩ and |*f*_{4}⟩, which are well separated in energy from competing states. The excitation targeting either |*f*_{2}⟩ or |*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 *γ*_{f} = 0.01*ω*_{0}, as in this regime $\omega f3\u2212\omega f2\u22483\gamma f$.

As shown in Refs. 51 and 52, 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 (*ω*_{1}, *ω*_{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 following Schmidt decomposition:^{65}

which is a weighted sum of *M* orthonormal modes ${\varphi 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 $\u2211jrj2=1$.

The Schmidt modes are useful in constructing the optimal *separable* or *classical* state, which contains no quantum correlations,

It can excite a fraction $r12$ of the optimal population.^{51} 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.^{25,49,66} It cannot, however, rely on entanglement, i.e., on the superposition of different products of modes as in (18).

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 (18): 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 |*f*_{1}⟩, we see that the classical pulse excites the $(\omega e1,\omega f1\u2212\omega e1)$ and $(\omega e3,\omega f1\u2212\omega e3)$ transitions, which are resonant with |*f*_{1}⟩. In addition, however, there is also substantial spectral weight at $(\omega e1,\omega e1)$ and $(\omega f1\u2212\omega e1,\omega f1\u2212\omega e1)$, i.e., at frequency combinations that are off-resonant for the given *γ*_{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 of 2 compared to the classical case [see Fig. 1(b)].

Since the natural linewidth *γ*_{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 of 10 (i.e., shorter lifetimes), $\gamma f\u2032=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 |*f*_{1}⟩, for instance, the square-like structure, visible in Fig. 1, of the classical pulse, for *γ*_{f}, cannot be resolved anymore for $\gamma f\u2032$ because the single-photon frequency distributions are much broader than the difference $(\omega f1\u2212\omega e1)\u2212\omega e1=\omega f1\u22122\omega e1$. This explains why the population induced by the shaped classical state is larger for $\gamma f\u2032$ than for *γ*_{f}. Similarly, when we try to optimally populate |*f*_{2}⟩, we now also obtain a larger population in |*f*_{3}⟩, as well as some in |*f*_{4}⟩. The latter is due to the classical pulse having a significant peak around $(\omega e3,\omega e3)$, which is broad enough to excite the $(\omega e3,\omega f4\u2212\omega e3)$ two-photon transition to |*f*_{4}⟩.

### C. Selectivity

Subsection IV B 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 |*f*_{1}⟩ and |*f*_{4}⟩ decrease (increase) with increasing *γ*_{f} when using an entangled (classical) two-photon state, and (iii) the target population in |*f*_{2}⟩ (|*f*_{3}⟩) decreases with increasing *γ*_{f}, while that of |*f*_{3}⟩ (|*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 |*f*_{2}⟩ and |*f*_{3}⟩,

worsens for broader linewidths.

Finally, an analysis of how the desired selectivity scales with the intermediate and target states’ linewidths (*γ*_{f} = 2*γ*_{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 |*f*_{2}⟩ is induced by selective optimal entangled states (15) to that achieved by shaped classical states of light (19). 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 |*g*⟩ → |*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 *γ*_{f}, in a qualitatively exponential fashion (top panel). Yet, *selective* optimization, as conceived in Sec. III, for classical and for quantum light, slows this decay down. For instance, at *γ*_{f} = 0.05*ω*_{0}, the indistinctive optimal (quantum) pulse will hardly discriminate |*f*_{2}⟩ against |*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 bottom panel of Fig. 3: this ratio actually *increases* with increasing *γ*_{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 *γ*_{f}.

## 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 work.^{51,52} 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 aggregates^{67} 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 molecular aggregates.^{33,34} 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 system 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 Ref. 53, 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 the said manifolds mix, instead will be a topic of future research.

## ACKNOWLEDGMENTS

This work was 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 support the findings of this study are available within the article.

## REFERENCES

The notation indicates the states (|*g*⟩ or |*e*⟩) of the two two-level atoms, followed by the number of photons left in the cavity mode.

*α*,

*α*′ ∈ {

*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.

Since $\mu ge3\mu e3f<\mu ge1\mu e1f$.