How to optimize the absorption of two entangled photons

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.


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 quantumenhanced 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). 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. , the deviation of the excited state decay rates' ratio from two (again in units of γe), see (17), and, here of central interest, the frequency correlations (see panel (b)) inscribed into the incoming pulses. (b) Density |Tt(ω 1 , ω 2 )| 2 ∝ |Φ(ω 1 , ω 2 )| 2 of the optimal two-photon wave function (31) in the space of rescaled frequencies (ω j − ωe)/γe, and parametrised by detuning ∆ = 5 and deviation δ = −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 |Φ(ω 1 , ω 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 (|g → |e ) and second (|e → |f ) bare transition's (see panel (a)) rescaled frequencies 0 and ∆, respectively, with distinct widths determined by ∆ and δ, see (32,33). The widths of the local maxima give a qualitative impression of the nonclassical frequency correlation between both incoming photons.

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 f and H 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 ω [42]. Choosing the latter, for each beam j we obtain annihilation a j (ω) and creation a † j (ω) operators satisfying the commutation relation [a j (ω), a † l (ω )] = δ jl δ(ω − ω ). In this framework the photon number operator for beam j reads n j =´∞ 0 dω a † j (ω)a j (ω), while the total Hamiltonian for both fields is H f = j´∞ 0 dω ωa † j (ω)a j (ω). 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 where 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 ∆ω much smaller than their central frequency ω 0 , i.e. ∆ω ω 0 . We can then employ the narrow bandwidth approximation, where we extend the range of all above frequency integrals to (−∞, ∞), substitute ω with ω 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 Together with its adjoint a † j (t), it satisfies the commutation relation [a j (t), a † l (t )] = δ jl δ(t − t ). The electric field operator (1) for beam j is finally re-expressed as so that the total (positive-frequency) electric field seen by the atomic target reads (we neglect any geometry dependence in the coupling factors [23]) As shown in fig. 1(a), the three non-degenerate electronic energy eigenstates of our target are |g , |e , and |f , with increasing energy. We define the origin of the energy axis to coincide with the energy of |g , hence the excited state energies are ω e and ω f , respectively. With these labels, the matter Hamiltonian is H m = ω e |e e|+ ω f |f 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 f + H m , and for nearresonant perturbative driving of the atomic transitions -reads 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 where the dipole matrix transition elements µ ge and µ ef between |g and |e , and between |e and |f , respectively, can be chosen real-valued.

Two-photon absorption amplitude
Our objective is to identify the optimal two-photon field state |Φ that drives the matter degrees of freedom from its initial (at time t 0 ) state |g into the target state |f (at time t), by TPA via |e . This is tantamount of maximizing the transition probability where |0 indicates the vacuum state of both injected fields, |f (t) = e iω f (t−t0) |f exhibits the explicit time dependence of the interaction picture with respect to H m (whereas |g , at energy zero, remains unaffected), and U I (t) is the time evolution operator in the interaction picture of H f + H m , given by the Dyson series Since |Φ 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 acting solely on the field degrees of freedom. *

Matter response function
Given the faint incoming field implied by (7) upon fixing the impinging two-photon state |Φ , 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), Using (4) and (5), the integrand becomes where we dropped the terms a i (τ 1 )a i (τ 2 ) which result from E + (τ 2 )E + (τ 1 ) by (4), since the sought-after pulse |Φ has only one photon in each mode. Equation (2) and an exchange of the frequency and time integrations gives where T t;t0 (ω 1 , ω 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;t0 (ω 1 , ω 2 ) must also encode the symmetry of (11) under exchange of the photon from the first or the second pulse being absorbed first. This will become explicit in (14).
To compactify the explicit expressions for T t;t0 (ω 1 , ω 2 ), we introduce the line shape functions where the dipole matrix element µ s−1 s connects the state s to the next energetically lower lying state s − 1 (e.g., |e to |g ), and γ s phenomenologically accounts for finite decay rates of |f and |e (|g , being the ground state, cannot decay). With (10,11,13) we can thus extract from (12) (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].
where (ω 1 ↔ ω 2 ) in the last line indicates identical terms with the frequencies ω 1 and ω 2 of the two photons exchanged.

Infinitely extended pulses
Equation (14) 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 → ∞, and, as a consequence, the real exponential factors in (14) vanish due to the excited state lifetimes, and we obtain the simpler expression [40] 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 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 e in (6), which, by linearity, can be carried directly throughout all derivations above, to obtain: The structure of the matter response function in (15) directly reflects the absorption process in the matter: while the term L e (ω 1 )+L e (ω 2 ) describes the transition |g → |e induced by either one of the two photons, the term L f (ω 1 + ω 2 ) correlates the two photons' frequencies by requiring that their sum be resonant with the two-photon transition |g → |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 ω f = 2ω e and γ f = 2γ e ), which can be independently driven by one single-photon pulse each, we introduce the (dimensionless) detuning ∆ and the deviation δ as as indicated in figure 1(a).

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: |ψ . The frequency representation ψ(ω), or wave function, of the one-photon state |ψ is its amplitude in the mode ω: where we introduced the continuous-mode single-photon state [42] and the continuous-mode creation operator a † (ω) 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: 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 which is equivalent to writing, in the Hilbert space,

Fluctuation-constrained optimization
Our goal is to maximize the population p f (t) in state |f at the time t, by appropriate choice of the two-photon state |Φ . The latter is determined [40] by the extrema of the functional where the dependence of p f (t) on the state |Φ can be unveiled by inserting (9) into (7), to obtain The second term in (23) constrains the distribution of the injected two photons over the incoming fields via the Lagrange multiplier λ: the expectation value n 1 n 2 limits the search to two-photon states where each beam is populated by one photon * . 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[|Φ ], its functional derivative, (21,22), must vanish, With (24) in (23), this requirement results in the eigenvalue problem which, with the definition turns into In addition, requiring the variation of the functional (23) with respect to the Lagrange multiplier to vanish, i.e. δJ/δλ = 0, enforces the normalization of the two-photon state, such that we arrive at as the only solution * . The maximal population can thus be expressed directly in terms of N , by (29) in (24), together with Φ|Φ = 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 N in terms of the system parameters [40]: 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 γ −1 e , the optimization problem will yield the optimal population (30) for arbitrary choice of t 0 and t.

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 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 |f at time t is given by the time-reversed two-photon state emitted by the three-level system initially (at t) prepared in |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 ∆ and deviation δ. The concentration (15), and discussed in 2.4, that the total frequency of the injected photons be resonant with the two-photon transition |g → |f . The frequency sum distribution is readily extracted from |Φ(ω 1 , ω 2 )| 2 , by changing variables to ω ± = ω 1 ± ω 2 and integrating over ω − : * Since |T T | is a one-dimensional projector, it has only one non-trivial eigenstate.
i.e. a Lorentzian of width γ f = (2 + δ)γ e , centred at ω f = 2ω e + ∆γ e . In figure 1(c), instead, we plot the single-photon distribution, i.e. the marginal distribution of |Φ(ω 1 , ω 2 )| 2 with respect to either frequency (here the first), since the optimal two-photon wave function (31) is symmetric: where we can identify two distinct peaks near * the transitions |g → |e (frequency ω e and line width γ e ) and |e → |f (frequency ω f − ω e = ω e + ∆γ e and line width 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 |e , and the second completes it to |f . When δ → −2, that is in the limit of a final state with vanishing linewidth γ f → 0, p sum (ω + ) tends to δ(ω + − ω 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 |g → |e or the |e → |f transition. In this limit, p 1 (ω) exhibits two peaks of equal (unit) width at ω e and ω f − ω e , as discussed above.

Entanglement entropy
When ∆ = δ = 0, i.e., by (17), ω f = 2ω e and γ f = 2γ 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 ϕ * k (ω 1 ), ψ * k (ω 2 ) and their non-increasingly ordered, non-negative Schmidt coefficients r k , where normalization via N in (31) implies k r 2 k = 1 (0 ≤ r k ≤ 1), and complex conjugation is inherited from that of T t (ω 1 , ω 2 ), also in (31). The entanglement encoded in Φ(ω 1 , ω 2 ) can be quantified by the entanglement entropy [45], We now inspect how the degree of frequency entanglement of the incoming pulses correlates with enhanced absorption probabilities.

Quantum enhancement
By (30), the maximal population of the state |f is 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 * The peaks do not exactly match the transition frequencies because |Le(ω 1 ) + Le(ω 2 )| 2 = |Le(ω 1 )| 2 + |Le(ω 2 )| 2 , and interference terms contribute to p 1 (ω).  (44), as a function of δ, in the limit ∆ → ∞ (achieved by γe → 0).
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 2 1 N ≤ N , where equality is achieved only for ∆ = δ = 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 of those two excitation probabilities. According to (35), strong entanglement S requires L 1 and thus small r 1 , which implies strong quantum enhancement E 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 δ (panel (a)) and detuning ∆ (panel (b)), respectively, as well as the maximally achievable values of entanglement and enhancement, S ∞ and E ∞ , respectively, in the limit of very large ∆ (panel (c)), which we discuss separately in section 3.5. * In panels (a,b), we observe a monotonic increase of E q with S, and a plateau of the achievable entanglement and quantum enhancement emerges for large values of the detuning ∆, which we address in the next subsection. The increase of S with ∆ and as δ → −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].

Maximal quantum enhancement
Let us now discuss the quantum enhancement E q in the limit ∆ 1 (or, more physically, ω f − 2ω e γ 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 ϕ k , ψ k in (34) are equal for each k, possibly up to a phase. For ∆ 1, as in figure 3, said modes appear as orthogonal superpositions of non-overlapping, complex-valued line shapes α k (centred at ω e ) and β k (centred at ω f −ω e = ω e +∆γ e ). Moreover, mutually orthogonal linear combinations of the same modes, e.g. ϕ 1,2 = (α 1 ± β 1 )/ √ 2 in figure 3, are associated with Schmidt coefficients r 1 ≈ r 2 , at least for finite ∆.
Let us now understand the connection between the Schmidt modes ϕ k , ψ k and the line shapes α k , β k . The latter can be considered the Schmidt modes of a response function where the first photon is resonant with the |g → |e transition, while the second completes the two-photon transition by inducing |e → |f . The construction of such a response function yields The analogy to (15) becomes clear by noticing that that is, T t is the symmetrized * version of Q t . For visualization, Q t (ω 1 , ω 2 ) and Q t (ω 2 , ω 1 ) describe, respectively, the top-left and bottom-right peaks of T t (ω 1 , ω 2 ) in figure 1(b).
Let us now write the Schmidt decomposition of the asymmetric response function (38) as where the prefactor ensures k s 2 k = 1. The Schmidt coefficients must be independent of ∆, 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 ∆ 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 Going back to the example of figure 3, if we assume r 1 = r 2 and multiply out the combinations proposed for ϕ k (ω 1 ) and ψ k (ω 2 ), for k = 1, 2, we obtain precisely the products of α 1 and β 1 on the right-hand side of (41). Notice that these combinations are not arbitrary: as discussed earlier, because T t is a symmetric complex function, the modes ϕ k and ψ k might differ by a phase that ensures the positivity of the singular values. All in all, we must satisfy ϕ 1 (ω 1 )ψ 1 (ω 2 ) + ϕ 2 (ω 1 )ψ 2 (ω 2 ) = α 1 (ω 1 )β 1 (ω 2 ) + β 1 (ω 1 )α 1 (ω 2 ). The right-hand side of (41) is a valid Schmidt decomposition only in the limit ∆ → ∞, where the function bases {α k } k and {β k } k are also mutually orthogonal †, 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 ∆ → ∞, eq. (41) is a valid Schmidt decomposition and we can identify The pairwise appearance of Schmidt coefficients r 1 = s 1 / √ 2, by (42), implies that the enhancement E ∞ , when ∆ 1, is twice the enhancement E a obtainable with the asymmetric matter response function (38): For the entropy, instead, we can write where S a is the entanglement entropy of (40). For any value of δ, then, the quantum enhancement induced by the optimal pulse (31), and the entanglement between the two photons' frequencies must be bounded, respectively, by E ∞ and S ∞ , which are plotted in figure 2 (c). For δ → −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.

Realistic pulses
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 |Σ 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 |Σ .

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 |Σ . If we continue to denote with |Φ the optimal state (29) analyzed in section 3, we intuitively want M 1 M 2 |Σ to most closely resemble |Φ .
In the frequency representation of section 2.5, we write a pulse shaping function as ω j |M j |ν j = M j (ν j , ω j ).
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 ω 1 + ω 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 |Σ : The factor N stems from (29), and clarifies the meaning of this equation: if N is the maximal population achievable by using the optimal state |Φ , then the population in the case of a shaped realistic state M 1 M 2 |Σ is reduced by the overlap between these two initial states.
The new functional to optimize -instead of (23) -is then 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 (47) we constrain the search to local operators that do not change the number of photons of undetermined, and potentially unnormalized, single-photon states |ψ j as per (18). The summation term in (47) therefore limits our search to unitary local operators: The single-photon states |ψ j do not carry further relevance beyond imposing the constraint (48), 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 (A.7), whose derivation we defer to Appendix A.

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 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 (49) in the derivation of Appendix A, we obtain, in place of (A.7), The functional A[M 1 , M 2 ] is defined in (A.1), and we have also introduced the effective response function W t as the product of the matter response function (15,31) with the (frequency representation of the) input state |Σ : Using the effective response function, the final state population (46) achievable via the diagonal pulse shaping functions M 1 and M 2 of (49) reads Given the fact that M j is a unitary diagonal operator, see (48), we can write with F j (ω j ) a real function. Substitution of this ansatz in (50) yields If we also cast the effective response function in polar form, the identity (54) is satisfied by functions F 1 and F 2 such that and by a Lagrange multiplier ψ j with Substituting (53) and (55) in (52), and assuming (56) is satisfied, we find the final state population achievable upon manipulating ("s" for "shaped") the input state |Σ : On the other hand, in the absence of shaping operators, that is with 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 W t , i.e. by the phase difference between the optimal |Φ and the realistic state |Σ . These phase differences are precisely what the optimal * pulse shaping functions compensate, by transforming W t (ω 1 , ω 2 ) into |W t (ω 1 , ω 2 )| via (53) with (56). To capture the enhancement obtained by pulse shaping we therefore * Applying the Hölder inequality [49] to (52), and using the unitarity of the shaper functions (48), we can write the inequality Because equality holds for the shaper functions satisfying (56), these must define the optimal solution.
introduce the optimization ratio as the ratio between optimized and unoptimized populations: . (61) To conclude, in the following examples we will identify the effective response function W t . If we can determine, via its phase, the argument F j (ω j ) of the pulse shaping operators, such that (56) is satisfied, then the optimal final state population is directly given by (58).

Example I: shaping down-converted photons
The first example we consider is depicted in figure 4(a): a pump photon of frequency ω p is split via SPDC into two photons of frequencies ω 1 and ω 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 |Σ given by where ω p is the frequency of the pump photon creating the entangled pair, which we assume to be on resonance with the two-photon transition (ω p = ω 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 (62). We consider frequency-degenerate down-converted photons, i.e. the individual photon wave packets are described by a real function G symmetric around ω f /2 [51]. The two photons are therefore equally detuned from either transitions |g → |e and |e → |f , since ω f /2 = ω e + (∆/2)γ e = (ω f − ω e ) − (∆/2)γ e . For compactness of notation, then, we shift the frequencies as Ω j = ω j − ω f /2. The effective response function (51) for this example then reads In this expression we dropped all prefactors, because the optimal pulse shaping operator M is determined by the phase of W t (Ω 1 , Ω 2 ), as defined in (55,56), and therefore only * depends on the detuning ∆ (see (17)). It is worth anticipating here that the time dependence of T t from (15) becomes a trivial phase factor e −iω f (t−t0) , once the Dirac delta is integrated over, as we do later; for this reason, we do not spell it out explicitly in (63).
In the shifted frequencies Ω j above, and for identical SLMs, i.e. M 1 = M 2 = M , Eq. (50) becomes where we integrated over the variables involved in the Dirac deltas of (63). Mirroring equations (55)-(57), then, the optimal pulse shaping function reads with S given in (55), while These expressions presuppose that the effective response function (63) is symmetric in Ω, and that G(Ω) is real, as assumed above.
In this example, the optimization ratio (61) reads and we plot its values in figure 5(a) for the individual photons described by the Gaussian profile where σ 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 |g → |e and |e → |f , of frequencies, respectively, ω e and ω e + ∆γ e , see figure 4(c). We excite the transitions with two photons with frequency profiles centred on ω f /2 = ω e + (∆/2)γ e and width σ. As discussed above, the two photons are always detuned by (∆/2)γ e from the electronic transitions. The probability of exciting the individual transitions depends, therefore, on the width σ of the single-photon frequency distributions.
If ∆ = 0, i.e. in the case of no detuning, each photon can resonantly drive either transition |g → |e or |e → |f . In this case, then, there is nothing to optimize. This can be shown formally with little algebraic manipulation: W t (Ω, −Ω) from (63) becomes a real function, and hence carries no phase to compensate via (56). When σ ∆, instead, the photon frequencies are so narrowly distributed around ω f /2 that they are always off resonant (by ∆/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 ∆, E opt reaches a maximum and then saturates for σ ∆, 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 σ of the photon pulses to match the detuning ∆, and plot the maximal final state population p f achievable by both unshaped (62) and shaped (65) Gaussian wave packets. Remarkably, the latter departs rather slowly from the maximal value N (30) achievable by the optimal state (29). Hence, in this case, shaping can almost saturate the optimal bound E q allowed by quantum mechanics.
Finally, the remaining open question is how optimal pulse shaping (65) of the realistic entangled state |Σ compares to the classical limit set by optimal separable pulses (36). As in figure 5(b), in figure 6 we match the width σ of the Gaussian pulse (68) to the detuning ∆, to cover both transition frequencies. As discussed for (63), the deviation δ (see (17)) does not affect W t (Ω, −Ω); the dependence on δ 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 linethese are the areas where an enhancement E q ≤ 1 with respect to the optimal separable pulse (36) is achieved).

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 downconversion 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 This is an appropriate model, for instance, for photon pairs created by type-I downconversion [53,54,55]. Here, α(ω) is proportional to the amplitude of the pump laser driving the down-conversion, and β(ω) denotes the phase-matching function. To apply the optimization procedure of section 4.2, we change variables to ω ± = ω 1 ± ω 2 and write the effective response function (51) as 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 ω + = ω 1 + ω 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 α(ω + ), i.e. α(ω + ) → M 1 (ω + )α(ω + ). Consequently, in this example we consider M 1 (ω + ) = M (ω + ) and M 2 (ω − ) ≡ 1, to write (50) as with We can now apply our previous results: Adopting the ansatz (53), and mirroring equations (55)-(57), we write ξ(ω + ) = |ξ(ω + )|e iϑ(ω+) . The optimal pulse shaping function is then with the Lagrange multiplier Finally, the optimization ratio (61) becomes where G denotes the cumulative function of the standard normal distribution [49].
In figure 8(b,c) we fix the quadratic phase at φ = 1, and tune the parameters σ and ζ characterizing the incoming pulse (69) via, respectively, α(ω + ) and β(ω − ). The former determines the width of the pump pulse, and hence should be compared to the deviation δ, 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 ∆ between the transition frequencies |g → |e and |e → |f , see figure 1(c) and (33).
As in section 4.3, we adjust the pulse parameters (σ, ζ) to the matter characteristics (δ, ∆) by setting * σ = 3γ f = 3γ e (2 + δ) and ζ = γ e (2 + ∆). This allows us to explore in figure 8 the same parameter space (δ, ∆) that was considered in section 3. In panel (a) we show the enhancement E 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 (69) upon shaping the pump with (73), and without shaping. As shown in section 3.4, the ideal two-photon state (15) maximally populates the state |f , and hence yields the largest E q for any fixed value of ∆ and δ. In panels (b) and (c) we also observe E q ≥ 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 (73) perfectly counteracts the chirp. * To choose the specific functional dependence of σ and ζ we initially set σ = Aγ f and ζ = Bγe(2 + ∆), 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.  (36), achievable via (a) the optimal entangled two-photon wave function (15), (b) the unshaped Gaussian pulse of (69), and (c) the optimized Gaussian via the pulse shaping function of (73). For each (δ, ∆) defining the matter response function (15), we attune the Gaussian pulse by setting its widths as σ = 3γe(2 + δ) and ζ = γe(2 + ∆) (see text). The quadratic phase in (76) is φ = 1.
On the left of the dashed line the various two-photon wave functions perform better (Eq ≥ 1) than the optimal separable pulses (36). Panel (a) reproduces the numerical results of [40].

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 downconversion. 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.
We remind the reader that B j and C j are also functionals of ψ j (ω j ) and ψ * j (ω j ), the wave function (and corresponding conjugate) representing |ψ j . Since the latter play the role of Lagrange multipliers, however, we do not write them as arguments on the left-hand side of (A.2,A.3).
We can then rewrite (47) in the frequency representation as where the normalization C j of the |ψ j states, being a number, does not carry any dependence on the pulse shaping functions, over which we are optimizing.