{'fraction_non_alphanumeric': 0.04180488779123861, 'fraction_numerical': 0.01866167852810943, 'mean_word_length': 4.480192205018686, 'pattern_counts': {'":': 1, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 81, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear optical signals and spectroscopy with quantum light

Conventional nonlinear spectroscopy uses classical light to detect matter properties through the variation of its response with frequencies or time delays. Quantum light opens up new avenues for spectroscopy by utilizing parameters of the quantum state of light as novel control knobs and through the variation of photon statistics by coupling to matter. We present an intuitive diagrammatic approach for calculating ultrafast spectroscopy signals induced by quantum light, focusing on applications involving entangled photons with nonclassical bandwidth properties - known as"time-energy entanglement". Nonlinear optical signals induced by quantized light fields are expressed using time ordered multipoint correlation functions of superoperators. These are different from Glauber's g- functions for photon counting which use normally ordered products of ordinary operators. Entangled photon pairs are not subjected to the classical Fourier limitations on the joint temporal and spectral resolution. After a brief survey of properties of entangled photon pairs relevant to their spectroscopic applications, different optical signals, and photon counting setups are discussed and illustrated for simple multi-level model systems.


I. INTRODUCTION
Nonlinear optics is most commonly and successfully formulated using a semiclassical approach whereby the matter degrees of freedom are treated quantum mechanically, but the radiation field is classical (Boyd, 2003;Scully and Zubairy, 1997). Spectroscopic signals are then obtained by computing the polarization induced in the medium and expanding it perturbatively in the impinging field(s) (Ginsberg et al., 2009;Hamm and Zanni, 2011;Mukamel, 1995). This level of theory is well justified in many applications, owing to the typically large intensities required to generate a nonlinear response from the optical medium, which can only be reached with lasers. Incidentally, it was shortly after Maiman's development of the ruby laser that the first nonlinear optical effect was observed (Franken et al., 1961).
Recent advances in quantum optics extend nonlinear signals down to the few-photon level where the quantum nature of the field is manifested, and must be taken into account: The enhanced light-matter coupling in cavities (Raimond et al., 2001;Schwartz et al., 2011;Walther et al., 2006), the enhancement of the medium's nonlinearity by additional driving fields (Chen et al., 2013;Peyronel et al., 2012), large dipoles in highly excited Rdberg states (Gorniaczyk et al., 2014;He et al., 2014), molecular design (Castet et al., 2013;Loo et al., 2012), or strong focussing (Faez et al., 2014;Pototschnig et al., 2011;Rezus et al., 2012) all provide possible means to observe and control nonlinear optical processes on a fundamental quantum level. Besides possible technological applications such as all-optical transistors (Shomroni et al., 2014) or photonic quantum information processing (Braunstein and Kimble, 2000;Franson, 1989;Jennewein et al., 2000;Knill et al., 2001;Kok et al., 2007;O'Brien et al., 2003;U'Ren et al., 2003), these also show great promise as novel applications as spectroscopic tools. Parameters of the photon field wavefunction can serve as control knobs that supplement classical parameters such as frequencies and time delays. This review surveys these emerging applications, and introduces a systematic diagrammatic perturbative approach to their theoretical description.
One of the striking features of quantum light is photon entanglement. This occurs between two beams of light (field amplitudes (Van Enk, 2005)) when the quantum state of each field cannot be described in the individual parameter space of that field. Different degrees of freedom can become entangled. Most common types of entanglement are: their spin (Dolde et al., 2013), polarization (Shih and Alley, 1988), position and momentum , time and energy (Tittel et al., 1999), etc. Entangled photon pairs constitute an invaluable tool in fundamental tests of quantum mechanics -most famously in the violation of Bell's inequalities (Aspect et al., 1982a(Aspect et al., , 1981(Aspect et al., , 1982b or in Hong, Ou and Mandel's photon correlation experiments (Hong et al., 1987;Ou and Mandel, 1988;Shih and Alley, 1988). Besides, their nonclassical bandwidth properties have long been recognized as a potential resource in various "quantum-enhanced" applications, where the quantum correlations shared between the photon pairs are thought to offer an advantage. For example when one photon from an entangled pair is sent through a dispersive medium, the leading-order dispersion is compensated in photon coincidence measurements -an effect called dispersion cancellation (Abouraddy et al., 2002;Franson, 1992;Larchuk et al., 1995;Minaeva et al., 2009;Steinberg et al., 1992a,b). In the field of quantum-enhanced measurements, entanglement may be employed to enhance the resolution beyond the Heisenberg limit (Giovannetti et al., 2004(Giovannetti et al., , 2006Mitchell et al., 2004). Similarly, the spatial resolution may be enhanced in quantum imaging applications Pittman et al., 1995), quantum-optical coherence tomography (Abouraddy et al., 2002;Esposito, 2016;Nasr et al., 2003), as well as in quantum lithographic applications (Boto et al., 2000;D'Angelo et al., 2001).
However, it is now recognized that many of these applications may also be created in purely classical settings: Some two-photon interference effects originally believed to be a hallmark of quantum entanglement can be simulated by postselecting the signal (Kaltenbaek et al., 2008). This had enabled quantum-optical coherence tomography studies with classical light . Similarly, quantum imaging can be carried out with thermal light (Valencia et al., 2005), albeit with reduced signal-to-noise ratio. When proposing applications of quantum light, it is thus imperative to carefully distinguish genuine entanglement from classical correlation effects.
A first clear signature of the quantumness of light is the scaling of optical signals with light intensities: Classical heterodyne χ (3) signals such as two photon absorption scale quadratically with the intensity, and therefore require a high intensity to be visible against lower-order linear-scaling processes. With entangled photons, such signals scale linearly (Dayan et al., 2004Friberg et al., 1985;Georgiades et al., 1995;Javanainen and Gould, 1990). This allows to carry out microscopy  and lithography (Boto et al., 2000) applications at much lower photon fluxes. The different intensity scaling with entangled photons has been first demonstrated in atomic systems by (Georgiades et al., 1995) and later by (Dayan et al., 2004, as well as in organic molecules (Lee and Goodson, 2006). An entangled two-photon absorption (TPA) experiment performed in a porphyrin dendrimer is shown in Fig. 1a). The linear scaling can be rationalized as follows: entangled photons come in pairs, as they are generated simultaneously. At low light intensity, the different photon pairs are temporally well separated, and the two photon absorption process involves two entangled photons of the same pair. It thus behaves as a linear spectroscopy with respect to the pair. At higher intensities, it becomes statistically more plausible for the two photons to come from different pairs, which are not entangled, and the classical quadratic scaling is recovered [ Fig. 1a)]. This is related to the Poisson distribution of photons.
The presence of strong time and frequency correlations of entangled photons is a second important feature, which we shall discuss extensively in the course of this review. Fig. 1b) shows the two-photon absorption signal of entangled photons in Rubidium vapor (Dayan et al., 2004). In the left panel, a delay stage is placed into one of the two photon beams, creating a narrow resonance as if the TPA resonance was created by a 23 fs pulse. However, as the frequency of the pump pulse which creates the photons is varied in the right panel, the resonance is also narrow in the frequency domain, as if it was created by a ns pulse. This simultaneous time and frequency resolution along non Fourier conjugate axes is a hallmark of the time-energy entanglement, and its exploitation as a spectroscopic tool offers novel control knobs to manipulate the excited state distribution, and thereby enhance or suppress features in nonlinear spectroscopic signals.
In a different line of research, the seminal photon coincidence counting experiments were turned into a spectroscopic tool by placing a sample into the beam line of one of the two entangled photons, and recording the change of the coincidence count rate (Kalachev et al., 2007(Kalachev et al., , 2008Kalashnikov et al., 2014;Yabushita and Kobayashi, 2004). We will examine related schemes for utilizing entanglement in nonlinear spectroscopy.
This review is structured as follows: In section I.A, we briefly give a the background and introduce the superoperator notation used in the following. In section II, we discuss FIG. 1 a) Linear entangled TPA rate and quadratic nonlinear random TPA rate in a porphyrin dendrimer at different entanglement times. The inset shows the dominant effect of the entangled TPA at low input flux of correlated photons. b) Two-photon absorption 5S 1/2 → 4D 5/2,3/2 in atomic Rubidium (Dayan et al., 2004) vs. the time delay between the two photons (left panel), and vs. the pump frequency (right panel).
properties of entangled photons, and present their impact on excited state distributions upon their absorption in a complex quantum system. In section III we provide general superoperator expressions for nonlinear optical signals, and review the proposed available setups. Finally, section IV presents a general classification of nonlinear optical processes induced by quantum light, clarifying under which conditions its quantum nature may play a role.

A. Liouville space superoperator notation
The calculation of nonlinear optical signals becomes most transparent by working in Liouville space (Mukamel, 1995;Mukamel and Nagata, 2011), i.e. the space of bounded operators on the combined matter-field Hilbert space. It offers a convenient bookkeeping device for matter-light interactions where signals are described as time-ordered products of superoperators.
Here, we introduce the basic notation. With each operator A acting on the Hilbert space, we associate two superoperators Mukamel, 2009b, 2010) A L X ≡ AX, which represents the action from the left, and the action from the right. We further introduce the linear combinations, the commutator superoperator and the anti-commutator This notation allows to derive compact expressions for spectroscopic signals. At the end of the calculation, after the time-ordering is taken care of, we can switch back to ordinary Hilbert space operators. The total Hamiltonian of the field-matter system is given by H 0 describes the matter. The radiation field Hamiltonian is given by Here we introduced the creation (annihilation) operators for mode s which satisfy bosonic commutation relations [a s , a † s ] = δ s,s , [a s , a s ] = [a † s , a † s ] = 0. In many situations, we will replace the discrete sum over modes by a continuous integral s → V (2π) 3 dω sD (ω s ) withD(ω s ) being the density of states that we assume to be flat within the relevant bandwidthsD(ω s ) D .
We assume a dipole light-matter interaction Hamiltonian, where V ν = V ν + V † ν the dipole operator of molecule ν with the summation running over all the molecules contained in the sample. V (V † ) is the excitation lowering (raising) part of the dipole. = E + E † denotes the electric field operator and E (E † ) its positive (negative) frequency components which can be written in the interaction picture with respect to the field Hamiltonian (Loudon, 2000) as which is written in the slowly-varying envelope approximation. Unless specified otherwise, E will denote the sum of all relevant field modes. In the following applications, we will neglect rapidly oscillating terms, by employing the dipole Hamiltonian in the rotating wave approximation, where

B. Diagram construction
We adopt the diagram representation of nonlinear spectroscopic signals which can be found, for instance, in . It bears close similarity to analogous methods in quantum electrodynamics (Cohen-Tannoudji et al., 1992). We will employ two types of diagrams which are associated with calculating the expectation value of an operator A(t) using the density matrix or the wavefunction. For details see appendix A. First, we evaluate it by propagating the density matrix ρ(t), Eq. (10) can be best analyzed in Liouville space (3): We write the time evolution of the joint matter plus field density matrix using a time-ordered exponential which can be expanded as a Dyson series, where H int− is a superoperator (3) that corresponds to the interaction Hamiltonian (7). The time-ordering operator T is defined by its action on a product of arbitrary superoperators A(t) and B(t), T orders the following products of superoperators so that time increases from right to left. The perturbative expansion of Eq. (11) to n-th order in H int generates a number of pathways -successions of excitations and deexcitations on both the bra or the ket part of the density matrix. These pathways are represented by double-sided ladder diagrams, which represent convolutions of fully time-ordered superoperator nonequilibrium Green's functions (SNGF) of the form ∞ 0 dt 1 · · · ∞ 0 dt n V (n) νn...ν1 (t − t n . . . , t − t n · · · − t 1 ) × E (n) νn...ν1 (t − t n . . . , t − t n · · · − t 1 ) .
The material SNGF's in Eq. (13) are similarly defined: where denotes the propagator of the free evolution of the matter system, ν j = L, R. Similarly SNGF may be obtained by replacing some V † by V and E by E † . This representation further allows for reduced descriptions of open systems where bath degrees of freedom are eliminated.
As an example, the set of fourth-order pathways for the population of state f in a three-level scheme shown in Fig. 8 are given in Fig. 2b. We will refer to these pathways repeatedly in the course of this review. This ladder diagram representation is most suitable for impulsive experiments involving sequences of short, temporally well-separated pulses ranging from NMR to the X-ray regimes (Hamm and Zanni, 2011;Mukamel, 1995). In such multidimensional experiments, the time variables used to represent the delays between successive pulses ) t 1 , t 2 , t 3 , ... serve as the primary control parameters. Spectra are displayed vs. the Fourier conjugatesΩ 1 ,Ω 2 ,Ω 3 , ... of these time variables.
As indicated in Fig. 2b), each interaction with a field also imprints the phase φ onto the signal. Filtering the possible phase combinations ±φ 1 ± φ 2 ± φ 3 ± φ 4 of the signal, known as phase cycling, allows for the selective investigation of specific material properties Keusters et al., 1999;Krčmář et al., 2013;Scheurer and Mukamel, 2001;Tan, 2008;Tian et al., 2003;Zhang et al., 2012a,b). Acousto-optical modulation discussed in Fig. 15 offers the possibility to achieve this selectivity even at the single-photon level. Phase cycling techniques have been successfully demonstrated as control tools for the selection of fixed-phase components of optical signals generated by multiwave mixing (Keusters et al., 1999;Tan, 2008;Tian et al., 2003;Zhang et al., 2012a,b). Phase cycling can be easily implemented by varying the relative inter-pulse phases using a pulse shaper, which is cycled over 2π radians in a number of equally spaced steps (Keusters et al., 1999;Tian et al., 2003).
Alternatively, rather than propagating the density matrix Eq. (10), we can follow the evolution of the wave function by expanding with Keeping track of the wave function results in different pathways which can be represented by loop diagrams. Rather than propagating of both the bra and the ket, we then place the entire burden of the time evolution on the ket, and write where and T −1 denotes the ati-time-ordering operator. Here the ket first evolves forward and then backward in time, eventually returning to the initial time (Keldysh, 1965;Schwinger, 1961). Back propagation of the ket is equivalent to forward propagation of the bra. The resulting terms are represented by loop diagrams, as is commonly done in many-body theory (Hansen and Pullerits, 2012;Mukamel, 2008;Rahav and Mukamel, 2010;Rammer, 2007). As can be seen from Fig. 2, this representation yields a more compact description (fewer pathways) of the signals, since the relative time ordering of ket and bra interactions is not maintained. In this example, the f-state population is given by the single loop diagram, which is the sum of the three ladder diagrams (as well as their complex conjugates). It is harder to visualize short pulse experiments in this representation, and due to the backward time propagation the elimination of bath degrees of freedom in an open system is not possible. Nevertheless, this representation proves most useful and compact for frequency domain techniques involving long pulses where time ordering is not maintained anyhow  and for many-body simulations that are usually carried out in Hilbert space (Dalibard et al., 1992).
The double sided (ladder) and the loop diagrams are two book keeping devices for field-matter interactions. The loop diagrams suggest various wavefunction based simulation strategies for signals (Dorfman et al., 2013a). The first is based on the numerical propagation of the wavefunction, which includes all relevant electronic and nuclear (including bath) degrees of freedom explicitly. A second protocol uses a Sum Over States (SOS) expansion of the signals. In the third, semiclassical approach a small subsystem is treated quantum mechanically and is coupled to a classical bath, which causes a time dependent modulation of the system Hamiltonian. The third approach for a wave function is equivalent to the stochastic Liouville equation for the density matrix (Tanimura, 2006), which is based on ladder diagram book keeping.

II. QUANTUM LIGHT
In this section, we first discuss classical light which is the basis for the semiclassical approximation of nonlinear spectroscopy. We then briefly survey the main concepts from Glauber's photon counting theory for a single mode of the FIG. 2 Diagrammatic construction of the double excitation probability f f in a three level model a). b) The single loop diagram representing the evolution of the wavefunction. On the left hand side, the vertical line represents ket and the right hand it corresponds to the bra. The wave function in this closed time path loop diagram propagates forward in a loop from the bottom branch of the diagram along the ket branch to the top of the diagram. It then evolves backward in time from the top to the bottom of the right branch (bra). This forward and backward propagation is similar to the Keldysh contour diagram rules (Keldysh, 1965). Horizontal arrows represent field-matter interactions. c) The corresponding three ladder diagrams for the evolution of the density matrix. Here the density matrix evolves forward in time upward from the bottom to the top of the diagram. In order to achieve doubly excited state population starting from the ground state the density matrix has to undergo four interactions represented by the absorption of two photons represented by four inwardly directed arrows. The three diagrams represent different time orderings of ket vs bra interactions. These are lumped together in the single loop diagram b) electromagnetic field. Finally, we discuss in detail the multimode entangled states of light, which will be repeatedly used in this review.

A. Classical vs. Quantum Light
Clearly, our world is governed by quantum mechanics, so on a fundamental level light is always quantum. Yet, in many situations a classical description of the light field may be sufficient (Boyd, 2003;Scully and Zubairy, 1997). For pedagogical reasons, before we describe properties of quantum light, we first discuss in some detail how the classical description of the field emerges from quantum mechanics, and -more importantly -under which circumstances this approximation may break down.
The coherent state of the field, i.e. the eigenstate of the photon annihilation operator at frequency ω, is generally considered as "classical". In this state the annihilation operator, and hence the electric field operator -has a non-vanishing expectation value, a(ω) = 0. In general, we can write a multimode coherent state as where the mode amplitudes α are given by In a normally ordered correlation function (all a † are to the right of all a), we may then simply replace the field operators in the correlation functions by by classical amplitudes (Mandel and Wolf, 1995), But a coherent state is not yet sufficient for the field to be classical: As we have seen in Eq. (13), the light field enters into our formalism through its multipoint correlation function E (n) νn...ν1 (t − t n . . . , t − t n · · · − t 1 ). These correlation functions can always be rewritten as the sum of a normally ordered term and lower order normal correlation functions multiplied with commutator terms For the semiclassical limit to hold, we must be able to neglect all the terms containing field commutators. This is typically the case when the intensity of the coherent state, i.e. its mean photon number is very large. More generally, probability distributions of coherent states are also classical (Mandel and Wolf, 1995): States, whose density matrix has a diagonal Glauber-Sudarshan Prepresentation, are considered classical. This includes classical stochastic fields, which show non-factorizing field correlation functions (for instance, for Gaussian statistics with E = 0 and E † E = 0). Such classical correlations can also be employed in spectroscopy (Turner et al., 2013), but are not the focus of the present review.
According to this strict criterion, any other quantum state must therefore be considered nonclassical. However, that does not mean that such state will also show nonclassical features, and the elucidation of genuine quantum effects in a given experiment is often very involved; see e.g. the discussion around dispersion cancellation with entangled photons (Franson, 1992;Steinberg et al., 1992b), or quantum imaging Valencia et al., 2005).
The interest in quantum light for spectroscopy is twofold: As stated in (Walmsley, 2015), "the critical features of quantum light [. . .] are exceptionally low noise and strong correlations". In the following section, we first discuss the noise properties of quantum light, i.e. quantum correlations within an electromagnetic mode. The rest of the chapter is devoted to the strong time-frequency correlations, i.e. quantum correlations between different modes, which is the main focus of this review.

B. Single mode quantum states
In confined geometries, the density of states of the electromagnetic field can be altered dramatically. In a cavity, it is often sufficient to describe the field by a single field mode with a sharp frequency ω 0 , which may be strongly coupled to dipoles inside the cavity. In this section, we shall be concerned with this situation, where we can write the cavity field Hamiltonian simply as where a (a † ) describes the photon annihilation (creation) operator of the cavity mode. We introduce the dimensionless cavity field quadratures which represent the real and imaginary part of the electric field, respectively. A whole family of quadratures X(θ) and P (θ) may be obtained by rotating the field operators as a → a e −iθ . Quantum mechanics dictates the Heisenberg uncertainty for any quantum state. Coherent states form minimal uncertainty states, in which the lower bound in Eq. (28) applies This coincides with the quantum fluctuations of the vacuum state. Any quantum state that features fluctuations in one of its quadratures below this fundamental limit 1/2 is therefore said to be squeezed. Note that due to the Heisenberg uncertainty (28) the conjugate quadrature must show larger fluctuations.
The simplest example thereof is the squeezed vacuum state where we have Thus, a squeezed vacuum state is also a minimum uncertainty state. These fluctuations also show up in nonlinear signals, as we will discuss in section III.A. In the remainder of this section, we are interested in multimode analogues of such squeezed states, and show how new correlations arise due to the multimode structure of the field. This is known as time-energy entanglement.

C. Photon Entanglement in multimode states
Before presenting some specific models for quantum light, we first discuss time-energy entanglement in a more general setting. An exhaustive review of continuous variable entanglement theory can be found in a number of specialized reviews (Kok et al., 2007;Lvovsky and Raymer, 2009). Here, we briefly introduce some important quantities of entangled photon pairs, in order to later distinguish genuine entanglement from other bandwidth properties, which could also be encountered with classical light sources.
We shall be concerned with pairs of distinguishable photons that live in separate Hilbert spaces H (1) field and H (2) field , respectively. They may be distinguished, e.g., by their polarization, their frequencies, or their wavevectors and propagation direction. A two-photon state may be generally expanded as whereΦ(ω a , ω b ) is the two-photon amplitude.
To analyze the properties of this state, it is convenient to map it into onto a discrete basis. This can be achieved by the singular value decomposition of the two-photon amplitude (Law and Eberly, 2004;Law et al., 2000; where {ψ k } and {φ k } form orthonormal bases, known as the Schmidt modes, andr k are the mode weights. These are obtained by solving the eigenvalue equations with and with We then recast Eq. (33 using Eq. (34) The Schmidt decomposition (34) may be used to quantify the degree of entanglement between photon pairs. Rewriting the singular values asr k = √ Bλ k , where B denotes the amplification factor of the signal and λ k the normalized set of singular values of the normalized two-photon state, with k λ 2 k = 1. A useful measure is provided by the entanglement entropy (Law et al., 2000) In quantum information applications, E(ψ) represents the effective dimensionality available to store information in the state. We call the state (33) separable, if the two-photon amplitude factorizes into the product of single-photon amplitudes . This implies that λ 1 = 1, and the entanglement entropy is E = 0. In this situation, no correlations exist between the two photons: Measuring the frequency of photon 1 will not alter the wavefunction of the other photon. Otherwise, the state is entangled. We then have λ i < 1, and correspondingly E > 0; measuring the frequency of photon 1 now reduces photon 2 to a mixed state described by the density matrix, the measurement outcome ω (0) a thus influences the quantum state of photon 2.

D. Entangled photons generation by parametric down conversion
Entangled photon pairs are routinely created, manipulated, and detected in a variety of experimental scenarios. These include decay of the doubly excited states in semiconductors (Edamatsu et al., 2004;Stevenson et al., 2007), four-wave mixing in optical fibers (Garay-Palmett et al., 2007 or cold atomic gases (Balić et al., 2005;Cho et al., 2014a). Here, we focus on the oldest and most established method for their production -parametric downconversion (PDC) in birefringent crystals (Kwiat et al., 1995;Wu et al., 1986).
We will introduce the basic Hamiltonian which governs the PDC generation process, and discuss the output fields it creates. These will then be applied to calculate spectroscopic signals. We restrict our attention to squeezed vacua, where the field modes are initially in the vacuum state, and only become populated through the PDC process (see also Eqs. (251) -(255)).
In the PDC setup, a photon from a strong pump pulse is converted into an entangled photon pair by the interaction with the optical nonlinearity in the crystal (as will be discussed in detail in section IV). A birefringent crystal features two ordinary optical axes (o), and an extraordinary optical axis (e) in which the group velocity of optical light is different. The PDC process can be triggered in different geometries: One distinguishes type-I (e → oo) and type-II (o → eo) downconversion (Shih, 2003), and more recently type-0 (Abolghasem et al., 2010;Lerch et al., 2013). The two photons show strong time and frequency correlations stemming from the conservation of energy and momentum: Since they are created simultaneously, the two photons are strongly correlated in their arrival time at the sample or detector, such that each individual photon wavepacket has a very broad bandwidth. At the same time, the sum of the two photon frequencies has to match the energy of the annihilated pump photon, which may be more sharply defined than the individual photons (pump photon bandwidth is typically below 100 MHz in visible range).
The created photon pairs are entangled in their frequency, position and momentum degree of freedom (Walborn et al., 2010). The setup may be exploited to control the central frequencies of the involved fields (Grice and Walmsley, 1997). For simplicity, we consider a collinear geometry of all fields. This greatly simplifies our notation, while retaining the timefrequency correlations, which are most relevant in spectroscopic applications. After passing through the nonlinear crystal, the state of the light field is given by (Christ et al., 2013(Christ et al., , 2011 The propagator U PDC is given by the effective Hamiltonian which creates or annihilates pairs of photons, whose joint bandwidth properties are determined by the two-photon amplitude Here, A p is the normalized pump pulse envelope, sinc(∆kL/2) denotes the phase-matching function, which originates from wavevector mismatch inside the nonlinear crystal, ∆k(ω a , denotes the wavevector of either the pump or beams 1 or 2 at frequency ω. The prefactor α which is proportional to the pump amplitude determines the strength of the PDC process (and hence the mean photon number). Hence, in contrast to ]P hi in the enangled state (33), the two-photon amplitude Φ is not normalized, but increases with the pump intensity. Typically, the wavevector mismatch ∆k(ω a , ω b ) depends very weakly on the frequencies ω a and ω b . It is therefore possible to expand it around the central frequencies ω 1 and ω 2 of the two downconverted beams. For type-I downconversion, the group velocities dk 1 /dω 1 and dk 2 /dω 2 are identical [unless the PDC process is triggered in a strongly non-degenerate regime (Kalachev et al., 2008)], and the expansion around the central frequencies yields (Joobeur et al., 1994;Wasilewski et al., 2006) The first two terms of Eq. (45) create correlations between the two photon frequencies, while the third determines the bandwidth of the individual photons.
In type-II downconversion, the group velocities of the two beams differ, and the wavevector mismatch may be approximated to linear order (Keller and Rubin, 1997;Rubin et al., 1994) ∆k where the two time scales T 1 = L(dk p /dω p − dk 1 /dω 1 ) and T 2 = L(dk p /dω p − dk 2 /dω 2 ) denote the maximal time delays the two photon wavepackets can acquire with respect to the pump pulse during their propagation through the crystal. Without loss of generality, we assume T 2 > T 1 . In the following applications we adopt type-II phase matching (46). The corresponding wavefunctions may be controlled to maintain frequency correlations, as we will elaborate in section II.G.
The following applications will focus on the weak downconversion regime, in which the output light fields are given by entangled photon pairs, Eq. (47) now takes on the form of the entangled two-photon state, Eq. (33), discussed above, whereΦ denotes the normalized two-photon amplitude. We will review the quantum correlations of the created fields upon excitation by cw pump, which simplifies the discussion, or finite-bandwith pump pulses. The two cases require different tools.

E. Narrowband pump
One important class of time-frequency entangled photon pairs is created by pumping the nonlinear crystal with a narrow bandwidth laser, where the pump spectral envelope in Eq. (44) may be approximated as Using Eqs. (46) and (48), the phase mismatch in Eq. (44) can be expressed as (Peřina et al., 1998) sinc where T ≡ T 2 −T 1 is the entanglement time, which represents the maximal time delay between the arrival of the two entangled photons. For a cw pump laser, the first photon arrives at a completely random time, but the second photon necessarily arrives within the entanglement time. This property has been exploited in several proposals to probe ultrafast material processes using a cw photon pair source (Raymer et al., 2013;Roslyak and Mukamel, 2009a). As mentioned before, entangled photons affect optical signals via their multi-point correlation functions (Mukamel and Nagata, 2011;Roslyak et al., 2009b). The relevant field quantity in most applications discussed below -describing the interaction of pairs of photons with a sample -is the four-point correlation function. Using Eqs. (48) and (49), for the entangled state (47) this correlation function can be factorized as Here, is known as the two-photon wavefunction (Rubin et al., 1994). Upon switching to the time domain, we obtain where rect(x)= 1 for 0 ≤ x ≤ 1 and zero otherwise, and N , N denote the normalization of the two-photon wavefunction. The physical significance of the entanglement time is now clear: it sets an upper bound for the arrival of the second photon, given the absorption of the first one. Note that Eq. (52) is symmetric with respect to t 1 and t 2 because each interaction occurs with the entire field E = E 1 + E 2 . In a situation where, e.g., 0|E 2 (t 2 )E 1 (t 1 )|ψ twin is measured the two-photon wavefunction is not symmetric.

F. Broadband pump
We next consider entangled light created by a pump pulse with a normalized Gaussian envelope, This is a realistic model in many experimental scenarios. The analysis of Eq. (42) becomes most transparent by switching to a basis obtained by the Schmidt decomposition (34) of the two-photon wavefunction. The mode weights r k are positive and form a monotonically decreasing series, such that in practical applications the sum in Eq. (34) may be terminated after a finite number of modes (For a cw pump defined above, an infinite number of Schmidt modes is required to represent the delta-function). In appendix B, we present approximate analytic expressions for the eigenfunctions {ψ k } and {φ k }. The following analysis is restricted to a PDC regime in which six-(and higher) photon processes may be neglected. Such corrections can affect the bandwidth properties at very high photon numbers (Christ et al., 2013). The linear to quadratic intensity crossover of signals will be further discussed in section II.J.3. We next introduce the Schmidt mode operators and which inherit the bosonic commutation relations from the orthonormality of the eigenfunctions {ψ k } and {φ k }. The trans-formation operator U PDC now reads (Christ et al., 2013(Christ et al., , 2011) The output state |ψ out is thus a multimode squeezed state with squeezing parameters r k , in which fluctuations of the collective quadrature of the two The mean photon number in each of the two fields is given by We will be mostly concerned with the time-frequency correlations in the weak downconversion regime, i.e. n ≤ 1, when the output state is dominated by temporally well separated pairs of time-frequency entangled photons. The multi-point correlation functions of state (57), which are the relevant quantities in nonlinear spectroscopy, may be most conveniently evaluated by switching to the Heisenberg picture, in which the Schmidt mode operators become (Christ et al., 2013(Christ et al., , 2011 The four-point correlation function then reads ) with and The first line in Eq. (62) shows the same structure as Eq. (50), in that the two absorption events at frequencies ω a and ω b (and at ω a and ω b ) are correlated. Indeed, in the weak pump regime, when r k 1, Eq. (63) reduces to the two-photon wavefunction of the pulsed entangled pairs, and sinh(r k ) cosh(r k ) r k , h 12 denotes the two-photon contribution to the correlation function , which should be distinguished from the autocorrelation contributions g 1 and g 2 . The ratio of the inverse entanglement time T −1 defined above and the pump bandwidth σ p determines the frequency correlations in Eq. (63): As shown in Fig. 3a), for σ p T −1 2 we recover the cw regime with strong frequency anticorrelations. In the opposite regime, when σ p T −1 2 [panel c)], the two photons show positive frequency correlations, and in between, there is a regime shown in panel b), in which the photonic wavefunction factorizes, and no frequency correlations exist. The field correlation function as depicted in Fig. 3 can be measured experimentally (Kim and Grice, 2005).
In spectroscopic applications, two-photon events involving uncorrelated photons from different pairs become more likely at higher pump intensities. These events are described by the functions g 1 and g 2 , which at low intensities scale as ∼ r 2 k , so that Eq. (63) with h 12 ∼ r k dominates the signal. As the pump intensity is increased, events involving photons from different pairs must be taken into account as well.
The various contributions to the correlation function behave differently with increasing photon number: They depend nonlinearly on the mode weights r k which in turn depend linearly on the pump amplitude. Thus, with increasing pump amplitude (and photon number) the few largest eigenvalues get enhanced nonlinearly compared to the smaller values, and fewer Schmidt modes contribute to Eqs. (63) and (64). This is shown in Fig. 4 where the two correlation functions are plotted for different mean photon numbers,n = 0.1, · · · 100. As the number of participating Schmidt modes is decreased, the frequency correlations encoded in h 12 are weakened, and h 12 broadens. Conversely, the bandwidth of the beams g 1 is reduced with increasingn.

G. Shaping of entangled photons
The ability to manipulate amplitude and phase of the ultrashort pulses allows to coherently control matter information in chemical reactions and other dynamical processes (Wollenhaupt et al., 2007). Pulse shaping allows to drive a quantum system from an initial state to a desired final state by exploiting constructive quantum-mechanical interferences that build up the state amplitude, while avoiding undesirable final states through destructive interferences (Silberberg, 2009). The most common experimental pulse shaping configuration is based on spatial dispersion and is often based on back-toback optical grating spectrometer which contains two gratings (see the top part of Fig. 5a). The first grating disperses the spectral components of the pulse in space, and the sec-ond grating packs them back together, while a pixelated spatial light modulator (SLM) applies a specific transfer function (amplitude, phase, or polarization mask), thereby modifying the amplitudes, phases, or polarization states of the various spectral components. Originally developed for strong laser beams, these pulse shaping techniques have been later extended to single photon regime (Bellini et al., 2003;Carrasco et al., 2006;Defienne et al., 2015;Pe'er et al., 2005;Zäh et al., 2008) allowing to control the amplitude and phase modulation of entangled photon pairs, thereby providing additional spectroscopic knobs.
An example of a pulse shaping setup is shown in Fig. 5a). A symmetric phase profile in SLM yields a single unshaped Gaussian pulse [ Fig. 5b) and c)]. The phase step results in a shaped pulse that mimics two pulses with ∼ 400 fs delay (Fig.  5d,e). By employing more complex functions in the SLM, one may produce a replica of multiple well separated pulses or more complex shapes. This may be also achieved using e.g. the Franson interferometer with variable phases and delays in both arms of the interferometer as proposed in (Raymer et al., 2013) (see section III.B.4). Beam splitters in the two arms can create four pulses, when starting with a single entangled photon pair.

H. Polarization entanglement
So far we have ignored the photon polarization degrees of freedom. In type-II PDC the two photons are created with orthogonal polarizations. Using a suitable setup, this allows for the preparation of Bell states of the form (Pan et al., 2012) where (|H , |V ) denote the horizontal (vertical) polarization, respectively. The fidelity for this state preparation is maximal, when the two photon wavepackets factorize. Similarly, type-I PDC allows for the creation of states of the form (|HH ± |V V )/ √ 2 (Pan et al., 2012). The polarization degrees of freedom offer additional control knobs which may be used to suppress or enhance the signal from (anti-)parallel or orthogonal dipoles in a sample system -in a quantum mechanical extension of polarized photon echo techniques (Voronine et al., 2006(Voronine et al., , 2007.

I. Matter correlations in noninteracting two-level atoms induced by quantum light
The entangled photon correlation functions, may be used to prepare desired distributions of excited states in matter. In an insightful article (Muthukrishnan et al., 2004), which triggered other work (Akiba et al., 2006;, it was argued that using time-ordered entangled photon pairs, two-body two-photon resonances, where two noninteracting particles are excited simultaneously, can be observed in twophoton absorption. The surprising consequence is that the The photon pairs are detected in the inverse process of PDC -sum frequency generation (SFG). The SFG photons are subsequently counted in a single-photon counting module. In order to demonstrate two-photon interference oscillations, a Mach-Zehnder interferometer is placed between the last prism and the SFG crystal. b) The SFG counts (circles) and the calculated second-order correlation function (line) of the unperturbed wavefunction as a function of the signalidler delay. c) the same for the shaped wavefunction. nonlinear response is non additive and does not scale as number of atoms N . Such cooperative (nonadditive) response is not possible with classical or coherent light fields. Arguments were made that the cooperatively is induced in two-photon absorption by the manipulation of the interference among pathways. One consequence of the interference in nonlinear response is that the fluorescence from one atom can be enhanced by the presence of a second atom, even if they do not interact. If true this effect could be an interesting demonstration of quantum nonlocality and the Einstein-Podolsky-Rosen (EPR) paradox. In this section, we shall calculate this two photon process with quantum light and show that quantum locality is not violated.
FIG. 6 (Color online) a) Scheme of two noninteracting two-level atoms A and B and corresponding many-body state diagram where ground state g corresponds to both atom in the ground state, a -atom A is excited and B is in the ground state, b -atom B is in the excited state and atom A is in the ground state, and ab -both atoms are in the excited state. Arrow directed upward represent two-photon excitation. b) Relevant set of diagrams corresponding to atom A being in excited state to ∼ |µA| 2 |µB| 2 in field-matter interactions. α, β run over a and b to account for all possible permutations in the excitation pathways.

Collective resonances induced by entangled light
We use superoperator formalism (Marx et al., 2008) to investigate how the two-atom excitation cross section depends on the properties of the photon wave function. Consider two noninteracting two-level atoms A and B coupled to the radiation field (see Fig. 6)a. We assume that the total field-matter density matrix is initially in a factorizable form: where the ρ A,0 (ρ B,0 ) corresponds to the density matrix of the atom A (B) and ρ ph,0 is the density matrix of the field. The time-dependent density matrix is given by Eq. (11), which in the present case reads If the radiation field is classical then the matter density matrix factorizes and atoms A and B remain uncorrelated at all times: with This result remains valid for quantum fields as long as all relevant normally ordered field modes are in a coherent state, and cooperative spontaneous emission is neglected so that all field modes behave classically (Glauber, 1963;Marx et al., 2008). Eq. (70) does not hold for a general quantum state. We define the reduced matter density matrix in the joint space w = tr ph (ρ). Upon expanding Eq. (69) order by order in the field operators and tracing over the field modes, we obtain for the matter density matrix where ν is summed over all possible pathways. Pathway ν has n νṼ A interactions and m νṼ B interactions. ρ ν A (ρ ν B ) are time-ordered products of system A (system B) operators and F ν (τ 1 , ...τ nν , τ 1 , ..., τ mν ) are time-ordered field correction functions. In each order of this perturbative order-byorder calculation, all the correlation functions are factorized between the three spaces. The factorization Eq. (70) no longer holds in general, and atoms A and B may become correlated or even entangled. Eq. (73) will be used in the following. Note that pathways with n ν = 0 or m ν = 0 are single-body pathways, where all interactions occur either with system A or with B. Our interest is in the two-body pathways, where both n ν and m ν contribute to collective response.

Excited state populations generated by nonclassical light
So far, we have not discussed how the overall excitation probability by entangled light sources compares to classical light with similar photon flux. This is relevant to recent demonstration of molecular internal conversion (Oka, 2012) and two-photon absorption with the assistance of plasmonics (Oka, 2015). This is a very important practical point, if quantum spectroscopy is to be carried out at these very low photon fluxes, in the presence of additional noise sources. We now discuss the correlations on the excitation probability induced by the interaction with nonclassical light for the present of two noninteracting two-level systems.
The doubly excited state population is given -to leadingorder perturbation theory in the interaction Hamiltonian (9) -by the loop diagram in Fig. 2a) which may be written as a modulus square of the corresponding transition amplitude with the transition amplitude between initial state |ψ and final state |ψ , In the following, we only consider a unique final state and will drop the subscript ψ .

Classical vs entangled light
We first evaluate Eq. (74) for a classical field composed of two modes α and β , which is switched on at t = t 0 : T (c) ab (t) is then given by (Muthukrishnan et al., 2004) T (c) where A νm = µ m E ν e i(ων − m )t0−γt0 − e i(ων − m)t−γt and γ → 0. A straightforward calculation of Eq. (75) shows several terms, where each of the pathway contains single-particle resonances a,b = ω α,β as well as a collective resonance a + b = ω α + ω β . While combining these terms the collective resonance disappear due to destructive interference and the final result factorizes as a product of two amplitudes that contain only single-particle resonances j = ω µ , j = a, b and µ = α, β and no two-photon resonances a + b = ω α + ω β (see ).
We first describe the properties of the entangled source before discussing the spectroscopy applications. Consider a field made of entangled photon pairs of a cascade state |ψ ent depicted in Fig. 7 which was considered by (Muthukrishnan et al., 2004). This state can be prepared when an atom is promoted to the doubly excited state which consequently decays spontaneously back to the ground state by emitting a cascade of two photons and described by the wave function: It is important that photon with momentum p comes first and interacts with upper transition whereas photon with q comes later and interacts with lower β -transition. The two-photon frequency ω p +ω q is narrowly distributed around ω α +ω β with a width γ α , the lifetime of the upper level, whereas the singlephoton frequencies ω p , ω q are distributed around ω β (ω α ) with a width of γ β (γ α ), the lifetime of the intermediate (highest) level. Maximum entanglement occurs for γ β γ α . Using Eq. (79), and assuming that the atoms A and B have the same distance from the cascade source, so that t R is the time retardation (with t R = |r R |/c) we obtain (Zheng et al., 2013) T ent where The first term in Eq. (80) contains single-particle resonances, where the two systems are separately excited. The second term represents collective two-photon resonances a + b −ω α −ω β which can only be distinguished for non identical atoms a = b .
Note, that the classical amplitude (77) scales quadratically with the field amplitude A νm , whereas the entangled amplitude (80) scales linearly. This reflects the fact that at low intensity (entangled case) two photon absorption is effectively a linear absorption of the entangled pair as discussed earlier.
More detailed analysis of the scaling is presented in Section II.J.3.
Comparing Eqs. (77) and (80), we see that two-photon resonances are induced by the lack of time ordering in the photonic field. To explain this, we calculate the marginal probability and the conditional probability for the two photon absorption process displayed in Fig. 7 (Muthukrishnan et al., 2004). These probabilities have the following meaning: the absorption of α is turned on at τ α = 0 and decays slowly at the rate γ α , while the absorption of β turns on at τ β = τ α and decays rapidly at the rate γ β . Thus, the two photons arrive in strict succession, α followed by β, with the time interval between the absorption going to zero when γ β γ α , i.e. in the limit of large frequency entanglement. Another related effect, which is discussed later in this section and which eliminates one-particle observables, is based on the lack of time ordering of the absorption of the two systems, while this interference effect originally described in (Muthukrishnan et al., 2004) is based on a lack of time ordering of the two photons. Only the single-body single-photon resonances remain.
We now turn to the spectroscopy carried out using the model of entangled light given by Eq. (79). To that end we turn to the excited state population of atom A, which is given by where the two terms correspond to the final state of atom B: p a0 and p ab represent B being in the ground and excited state, respectively. We shall only calculate the ∼ µ 2 A µ 2 B contributions to p a , which are relevant to our discussion. The relevant set of diagrams is depicted in Fig. 6b. Diagram i contains three interactions with ket-and one with bra-and represents p a0 whereas diagram ii corresponds to p ab . The exact formulas read off these diagrams are presented in . Unlike the p a1 (t) which contains only normally ordered field operators [see Eq. (74)], p a0 (t) contains also non normally ordered contributions. These can be recast as a normally ordered correlation plus a term that includes a commutator: The second term involves a commutator of the field which is a c-number. For certain type of states of the field, the commutator remains dependent upon the state of the field and this term has to be evaluated exactly (for instance, in the case of coherent state this effect is responsible for the revival of damped Rabi oscillations discussed below (Rempe et al., 1987). For other types of states of the field, the commutator becomes independent of the external field (e.g. Fock state) and therefore represents spontaneous emission. The spontaneous emission pathways introduce a coupling between the two systems, since a photon emitted by system B can be absorbed by system A. This coupling has both real (dipole-dipole) and imaginary (superradiance) parts. These couplings will obviously result in collective signals which involve several atoms.
In the case of a classical field, one can neglect the spontaneous contributions and only include the stimulated ones which results in p a0 (t) = −p ab (t). The population p a (t) therefore vanishes, and is not affected by the collective resonances. We thus do not expect any enhanced fluorescence from A. The two-photon absorption signal is a product of individual excitation probabilities for atoms A and B, and shows no collective resonances.

Collective two-body resonances generated by illumination with entangled light
We shall focus on the change of the photon number for the two-body part ∆n ph . This will depend on the following probabilities: the excitation probability of p ab ,which means that two photons are absorbed (counts twice), and the excitation probabilities of only system A p a and of only system B p b (we assume that both photons have equal frequency and are resonant to the two-photon absorption): Now, in the stimulated emission and TPA pathway absorption, we know that p ab (t) and p a0 (t) cancel, and that p ab (t) and p b0 (t) cancel. So for the stimulated pathways, Eq. (88) vanishes. Since the field-matter interaction Hamiltonian connects the photon number with the excitation probability number, the photon number itself is a single-particle observable like the population of state |a or state |b , and therefore vanishes. However, collective resonances in p ab (t) can be revealed in two-photon counting (Hanbury-Brown-Twiss measurements) (Brown, 1956). For our entangled photon state, Eq. (79), the change in the photon-photon correlation ∆ n ph is attributed to any buildup of probability, that the either atom A, or atom B is excited or both atoms are excited, which will cause a reduction of the photon-photon correlation: Since p ab (t) enters twice [unlike Eq. (88)] the stimulated contributions can only cancel with one of the two other contributions p a0 (t) or p b0 (t), and we have The interference mechanism, which caused the cancellation for the stimulated signal and the photon number, does not lead to a full cancellation, and two-photon absorption involving both systems might be observed.
a. Two-photon two-atom problem with non entangled quantum states We now discuss whether entangled light is essential to the collective resonances. Whether or not photon entanglement is essential for these properties is the subject of current debate. In (Zheng et al., 2013), the authors presented a detailed analysis of the cross section created by entangled pulses, and compared this with "correlated-separable" states in which the entanglement is replaced by classical frequency correlations. They concluded that it is the frequency anticorrelations [see Fig. 4a)], and not the entanglement per se, which is responsible for the enhancement of the cross section. The effect of enhanced two-photon absorption probability has been later shown to come from the entangled matter/field evolution that occurs with any quantum light (which is not necessarily entangled) (Zheng et al., 2013) which is consistent with earlier demonstrations of (Georgiades et al., 1997). Starting from entangled pure quantum state having a density matrix ρ 0 = |ψ ent ψ ent | of matrix elements ρ kk qq = 1 k , 1 q |ρ 0 |1 k 1 q one can construct others that have the same mean energy and the same single-photon spectrum, and hence that would give the same transition probabilities for a single-photon resonance. We choose a special case of the states that originates from entangled state (79) that will allow a quantitative evaluation of the role of correlations. It is defined as It is the diagonal part of ρ 0 . It has lost any temporal field coherence and is time-independent. It is actually a correlated separable state (Duan et al., 2000), in which the quantum correlations are replaced by a purely classical frequency distribution. It gives rise, however, to correlations between its two parties. Using the entangled state (79), Eq. (91) reads The state (92) corresponds to an atomic cascade for which the starting time is random, thereby averaging to zero all the offdiagonal time-dependent terms in the density matrix. It gives rise to the following transition probability: At exact two-photon two-atom resonance, we have p 1 p 0 γ β γ −1 α c 2 t 2 /(∆L) 2 . Note that P 1 depends on time, as can be expected in a situation where the detecting atoms, which have an infinite lifetime, are submitted to a stationary quantum state, and therefore to cw light. In order to compare p 1 to corresponding probability (83) for entangled light, which is induced by a pulse of light, we need to fix an interaction time t. One then obtains at t = L/c and at exact resonance We thus find that a correlated-separable state like ρ 1 can induce the two-photon two-atom transition similar to the entangled cascade state. Note that even though ρ 1 is not entangled, it has quantum properties, being a mixture of single-photon states which are highly nonclassical.
In summary, the two-photon absorption probability of noninteracting atoms can be enhanced compared to classical light by using special types of quantum light. This quantum light does not require entanglement, rather it is necessary to have spectral anticorrelations which can be achieved in e.g. correlated separable states. This is a consequence of the stationary system of two noninteracting atoms which has no dynamical properties.
In the following section we consider more complex material system that is subject to various relaxation channels, e.g. transport between excited states therefore adding temporal scale to the problem. We show that in this case, entanglement is required in order to achieve both high spectral and high temporal correlations which is not possible by classical correlated light.

Model system
The standard calculation of the nonlinear response to classical light assumes that the matter is made up of N noninteracting particles (atoms or molecules) in the active zone, such that the matter Hamiltonian may be written as the sum over the individual particles [see Eq. (7)] The individual nonlinear susceptibilities or response functions of these atoms then add up to give the total response. The |gi |e 1 i |e 2 i Level scheme of the multilevel model employed in this review: Two singly excited states e1 and e2 with energies 11, 000 cm −1 and 11, 500 cm −1 are coupled to doubly excited states f1 and f2, as indicated. Furthermore, e2 decays to e1 within 1/k 30 fs.
nonlinear response becomes then a single-body problem and no cooperative resonances are expected. It is not obvious how to rationalize the ∼ N scaling for noninteracting atoms had we chosen to perform the calculation in the many-body space. Massive cancellations of most ∼ N (N − 1) scaling lightmatter pathways recover in the end the final ∼ N signal scaling (Spano and Mukamel, 1989). When the atoms are coupled, the calculation must be carried out in their direct-product many-body space whose size grows exponentially with N (∼ n N dimensions for n-level atoms). The interatomic coupling can be induced by the exchange of virtual photons leading to dipole-dipole and cooperative spontaneous emission, or superradiance . In molecular aggregates, the dipole interaction between its constituents creates inherent entanglement on the level of the quasiparticles (Mukamel, 2010), which shifts the doubly excited state energies, and redistributes the dipole moments, in which the individual Hamiltonians H ν are (Abramavicius et al., 2009) where B m (B † m ) describes an excitation annihilation (creation) operators for chromophore m. These excitations are hard-core bosons with Pauli commutation rules (Lee et al., 1957): To describe two level sites, which cannot be doubly excited according to the Pauli exclusion we set ∆ m → ∞. In condensed matter physics and molecular aggregates the many particle delocalized states are called excitons. The model Hamiltonian (96) can represent e.g. Rydberg atoms in optical lattice or Frenkel excitons in molecular aggregates. For considering (100), prepared by the absorption of entangled photons with pump frequency ωp = 22, 160 cm −1 . b) the same, with ωp = 24, 200 cm −1 . [from (Schlawin et al., 2012)] four-wave mixing processes, we have truncated the Hamiltonian H ν in Eq. (96) at the doubly excited level, such that it diagonalization reads where g indexes the ground state, e the singly, and f the doubly excited states, which will be the central concern. We will consider processes in which the signals factorize into ones stemming from individual constituents, as well as collective signals from two or more constituents. We illustrate basic properties of the interaction of entangled photons with complex quantum systems at the hand of the simple multilevel model depicted in Fig. 8. It consists of two excited states e 1 and e 2 , a two doubly excited states f 1 and f 2 . The higherenergy excited state e 2 decays to e 1 within few tens of femtoseconds. We shall focus on doubly-excited states created by the absorption of pairs of photons.
Similarly to the noninteracting case (75), the population in the final state f at time t is given by the loop diagram in Fig. 2, which can be written as the modulus square of a transition amplitude. Tracing out the matter degrees of freedom, the amplitude can be written as a nonlinear field operator, which may be evaluated for arbitrary field input states. Specif-ically, using Eq. (52), we obtain where we used N = N A p which is linear with respect to classical pump amplitude. For comparison, the classical probability governed by excitation by two classical fields with amplitudes A 1 and A 2 is given by which is quadratic in the field amplitudes, unlike Eq. (100). The second line denotes the same quantity, but with the beam frequencies ω 1 and ω 2 interchanged. We will further discuss more on the intensity scaling. Eq. (100) reflects the entangled photon structure described earlier: The pump frequency ω p = ω 1 + ω 2 is sharply defined such that the ω f -resonance is broadened only by the state's lifetime broadening and pure dephasing rate γ f . The strong time correlations in the arrival time creates a resonance of the form (exp(iωT ) − 1)/ω, which, for very short entanglement times is independent of the frequency: (exp(iωT ) − 1)/ω T + O(T 2 ). This implies that the intermediate e-states effectively interact with broadband light, whereas the f -states interact with cw light. Put differently, thanks to their large bandwidth, the entangled photons can induce all possible excitation pathways through the e-manifold leading to a specific selected f -state (Schlawin et al., 2012). Tuning the pump frequency ω p allows to select the excitation of specific wave packets: Fig. 9 depicts the density matrices for doubly excited states induced by the absorption of entangled photon pairs with different pump frequencies. This behavior remains the same for pulsed excitation in the strong frequency anticorrelations regime (Schlawin, 2015).

Control of energy transfer
When the material system is coupled to an environment which causes relaxation among levels, the distributions of excited states populations may no longer be described by the wavefunction and transition amplitudes, and the density matrix must be used instead. As shown in Fig. 2, the loop diagram, which represent the transition amplitudes needs to be broken up into three fully time-ordered ladder diagrams (and their complex conjugates). These are given in the time domain , Here, Γ collectively denotes the set of control parameters of the light field. The matter correlation functions in Eqs. (102) -(104) are given by Liouville space superoperator correlation functions, which translate for the integration variables in . The total f -population is normalized to unity at each frequency, i.e. f p f (t) = 1.
The entangled photon field correlation functions are given by the Fourier transform of Eq. (62). We first restrict our attention to the weak downconversion regime,n 1, in which the autocorrelation functions g 1,2 , Eqs. (64) and (65), may be neglected.
The essential properties of entangled photon absorption may be illusttrated by using the simple model system introduced in section II.J.1 where the correlation functions (102) -(104) can be written as the sum-over-states expressions given in appendix C. To excite state f 2 faithfully, the intermediate decay process e 2 → e 1 needs to be blocked (see Fig. 8), and one needs to select the final state f 2 spectrally. This can be achieved with entangled photons as shown in Fig. 10, where we depict the relative population in each final state vs. the pump frequency ω p . At each frequency the total population is normalized to unity, i.e. p f1 + p f2 = 1. By choosing a spectrally narrow pump bandwidth σ p = 100 cm −1 in combination with a short entanglement time T = 10 fs, almost 90 % of the total f -population can be deposited in the state f 2 at ω p = 22500 cm −1 (and ∼ 95 % in f 1 at ω p = 21800 cm −1 ).
The same degree of state selectivity may not be achieved with classical laser pulses: For comparison, the f -manifold populations created by classical laser pulses with bandwidths σ = 100 cm −1 (same spectral selectivity) and with σ = 1000 cm −1 (same time resolution in the manifold of singly excited states) are shown as well in Fig. 10. In the former case, the intermediate relaxation process limits the maximal yield in f 2 to ∼ 35 % (with ca. 65 % population in f 1 at ω p = 22500 cm −1 ), in the latter case, the lower spectral resolution limits the achievable degree of control over the fpopulations.
For very short entanglement times, transport between various excited states in e-manifold may be neglected, and the f -populations may as well be calculated using transition amplitudes -thus greatly reducing the computational cost (Schlawin, 2015). At the same time, by varying the entanglement time becomes possible to probe subsets of transport pathways via the selection of specific f -states in the detected optical signal.
The properties described above may also be observed in more complex systems such as molecular aggregates , where the number of excited states is much higher, but the described physics is essentially the same: Ultrafast relaxation processes limit the classically achievable degree of selectivity due to the trade-off between spectral and time resolution for every single absorption process. With entangled photons, this trade-off only limits the overall twophoton process, but not each individual transition. This may be understood in the following way: In the absorption of entangled photons, the light-matter system becomes entangled in between the two absorption events, but it remains separable at all times for classical pulses so that the energy-time uncertainty only applies to the entire system.
Finally it is worth noting an additional advantage of using entangled photons for two-photon excitation in molecules pointed out by (Raymer et al., 2013). In many aggregates the doubly-excited state undergoes rapid nonradiative decay to the singly excited states. This may occur on the same time scale of the excited-state transport (Van Amerongen et al., 2000). For weakly anharmonic aggregates where the transition frequencies f → e and e → g are close, it is hard to discriminate between the two-photon-excited and single photon excited fluorescence. Entangled light can help isolate the doubly-excited state population by monitoring the transmitted single-photon pathways with a high-quantum-efficiency detector, that can partially rule out one-photon absorption.

Scaling of two photon absorption with pump intensity
Eq. (101) shows that TPA probability |T f g | 2 scales quadratically with pump intensity A 2 1 A 2 2 for classical light but Eq. (100) shows linear scaling A 2 p with weak entangled light. As the pump pulse intensity is increased, the two-photon state (33) no longer represents the output state (57), and the autocorrelation contributions g 1,2 to Eq. (62) must be taken into account. To understand how this affects the excited state distributions, we first simulate in Fig. 11 the relative contribution of the two-photon term h 12 in Eq. (62) to the population in f 2 , when the pump frequency is set at resonance, ω p = 22500 cm −1 , and its variation with the mean photon number (59). Clearly, in the weakly entangled case this contribution drops rapidly, and already atn = 0.5 the autocorrelation contribution (i.e. excitations by uncorrelated photons from different downconversion events) dominates the signal. In contrast, in the strongly entangled case the contribution drops much more slowly, and -as can be seen in the insetstill accounts for roughly a third of the total absorption events even at very high photon numbern 100. As described in , this may be attributed to the strong time correlations in the former case: For T = 10 fs, one can fit ten times as many photon pairs into a time period compared to T = 100 fs, before different pairs start to overlap temporally. In addition, the strong frequency correlations discussed in the previous section imply that, on resonance with the final state f 2 , the coherent contribution is enhanced even further.
The situation becomes even more striking in panel b when the intermediate state e is shifted away from the entangled photon frequencies), where we set it to ω e2 = 18, 500 cm −1 , and only the two-photon transition g → f 2 is resonant. In the weakly entangled case (solid line), little has changed, and parity between the two-photon and the autocorrelation contribution is reached again for n ∼ 2. In the strongly entangled regime, however, the coherent contribution remains close to unity, and dominates the signal even atn 500, as can be seen in the inset.
The strong frequency correlations of entangled photon pairs lead to the enhancement of the two-photon contribution h 12 compared to the autocorrelation contributions g 1,2 -even for very high photon fluxes. Consequently, The nonclassical excited states distributions created by entangled photons can still dominate the optical signal even when the mean photon number in each beam greatly exceeds unity (Schlawin, 2015). The linear to a quadratic crossover of the two-photon absorption rate with the pump intensity ( Fig. 1a) which is traditionally regarded as a crossover from a quantum to a classical regime (Lee and Goodson, 2006), is not necessarily a good indicator for this transition after all, as pointed out by (Georgiades et al., 1997). The time-frequency correlations of entangled pairs may be harnessed at much higher photon fluxes -as long as the coherent contribution h 12 dominates the incoherent contributions g 1,2 .

III. NONLINEAR OPTICAL SIGNALS OBTAINED WITH ENTANGLED LIGHT
We now revisit the excited state populations created in matter following the absorption of entangled photon pairs in section II, and present optical signals associated with these distributions. We first derive compact superoperator expressions for arbitrary field observables, which naturally encompass standard expressions based on a semiclassical treatment of the field (Mukamel, 1995).
Using Eq. (11), the expectation value of an arbitrary field operator A(t) is given by Here, the first term describes the expectation value in the absence of any interaction with the sample, which we set to zero. The second term describes the influence of the sample on the expectation value. It may be represented more compactly, when the interaction Hamiltonian is not written in the RWA approximation (Schlawin, 2015), where we recall ( Eq. (111) may be used to calculate arbitrary optical signals induced by entangled photons, as will be demonstrated in the following.

A. Stationary Nonlinear signals
Here, we consider the situation described in section II.B, where the matter system is either coupled to a single cavity mode, and where only a single quantum mode of the light is relevant. This is most prominently the case in the strong coupling regime of cavity quantum electrodynamics (Walther et al., 2006), or when one is interested in steady state solutions, where only a single (or few) frequency(ies) of the field is(are) relevant.
The intricate connection between nonlinear optical signals and photon counting was first made explicit by (Mollow, 1968), who connected the two-photon absorption rate of stationary fields with their G (2) -function 2 is the coupling element evaluated at the central frequency ω 0 of the stationary light field. The transition rate loosely corresponds to the time-integrated, squared two-photon transition amplitude (99) we derived earlier. If the lifetime of the final state γ f is much smaller than the bandwidth of the field ∆ω, γ f ∆ω, one may even neglect the time integration, and replace This implies, as pointed out in (Gea-Banacloche, 1989), that strongly bunched light can excite two-photon transitions more efficiently than classical light with identical mean photon number. It further implies that the two-photon absorption rate scales linearly in the low gain regime, even though the single mode squeezed state does not show time-energy entanglement in the sense of section II.C. An experimental verification thereof is reported in (Georgiades et al., 1995(Georgiades et al., , 1997. A more recent proposal (López Carreño et al., 2015) aims to employ nonclassical fluctuations contained in the fluorescence of a strongly driven two-level atom: As is well known, the driven two-level atom's fluorescence develops side peaks -known as the Mollow triplet (Scully and Zubairy, 1997). By driving polaritons -strongly coupled light-matter states in a cavity -with this light, it is predicted to allow for precise measurements of weak interactions between polaritons "even in strongly dissipative environments".

B. Fluorescence detection of nonlinear signals
Fluorescence is given by the (possibly time-integrated) intensity A(t) = E † (t)E(t), when the detected field mode is initially in the vacuum. It is obtained by expanding Eq. (111) to second order in the interaction Hamiltonian with the vacuum field (Mukamel, 1995), We first investigate the time-integrated f → e fluorescence signal (113), following excitation by either entangled photons or classical pulses. Using Eqs. (113), and (102)-(104), we may readily evaluate the signal as 2 We change the notation to match the remainder of this review.

Logarithmic cross section
Energy (eV) FIG. 12 Entangled virtual-state spectroscopy: TPIF signal vs. the Fourier transform of an additional time delay between the two photons (see text) in hydrogen [taken from ].
In a different scenario, the two-excitation state may rapidly decay nonradiatively -e.g. via internal conversion -and the e − g fluorescence is detected, Fluorescence signals are proportional to excited state populations. Therefore, they are closely related to the excited state distributions discussed in sections II.I to II.J.3. This is not necessarily the case for absorption measurements, as will be shown in section III.D.

Two-photon absorption vs. two-photon-induced fluorescence
Before reviewing fluorescence signals induced by entangled photons, we comment on some ambiguity in the nomenclature of nonlinear signals such as fluorescence that depend on the doubly excited state population p f . Since population is created by the absorption of two photons, the signal is often termed two-photon absorption. However, it does not pertain to a χ (3) -absorption measurement, which will be discussed in section III.D. We will refer to signals measuring p f as twophoton-induced fluorescence and to two-photon resonances in χ (3) as two-photon absorption.
In section III.C.1, we will further show that for off-resonant intermediate state(s) e -as is the case in most experimental studies to date (Dayan et al., 2004Guzman et al., 2010;Harpham et al., 2009;Lee and Goodson, 2006;Upton et al., 2013) -the two signals carry the same information.

Two-photon induced transparency and entangled-photon virtual-state spectroscopy
Entangled virtual-state spectroscopy, proposed in , suggests a means to detect far off-resonant intermediate states in the excitation of f -states by employing entanglement-induced two-photon transparency (Fei et al., 1997): The transition amplitude (100) of entangled photons created by a cw-source oscillates as a function of the entanglement time T . For infinite excited-state lifetimes (γ e = 0), the excitation of f is even completely suppressed whenever (ω 1 − ω e )T = n × 2π ∀n ∈ N -the medium becomes transparent. Saleh et al. had proposed to turn this counterintuitive effect into a spectroscopic tool by sending one of the two photons through a variable delay stage τ , which simply amounts to evaluating the two-photon wavefunction (52) 0|E(τ 2 + τ )E(τ 1 )|ψ twin . For degenerate downconversion the transition amplitude (100) then reads where we introduced the detuning ∆ e = ω p /2 − ω e . Fourier transform of the TPIF signal ∝ |T f g (τ )| 2 with respect to τ reveals different groups of resonances shown in Fig. 12 for the TPIF signal in hydrogen: Group A denotes resonances of differences between intermediate states, ∆ e − ∆ e , group B of detunings ∆ e , and group C of the sum of detunings ∆ e + ∆ e . Such signals were simulated in molecular system (Kojima and Nguyen, 2004), and similar resonances can be identified in absorption TPA measurements, as we will discuss in section III.C.1. However, it turns out that the effect depends crucially on the tails of the spectral distribution in Eq. (49) (de J Leon- Montiel et al., 2013). The sinc-function in the two-photon wavefunction (51) has a long Lorentzian tail ∼ 1/ω, which covers extremely far off-resonant intermediate states. When the sinc is replaced by Gaussian tails, the resonances vanish. The effect is thus caused by the long spectral tails, and is not intrinsically connected to entanglement.

Fluorescence from multi-level systems
The TPIF signal, Eq. (114), directly reflects the doubly excited state population distributions created by the absorption of entangled photon pairs. This may be seen in simulations of the TPIF signal from Frenkel excitons in a molecular aggregates (Schlawin et al., 2012). In Fig. 13a), we present the TPIF signal resulting from the decay of the f -state distributions in Fig. 10. Clearly, the signal peaks whenever a two-excitation state is on resonance with the pump frequency ω p . The signal from the two f -states has approximately the same strength. In contrast, the TPIF signal created by classical pulses shown in panel b) has a strong resonance only when the light is on resonance with the state f 1 , whereas f 2 can only be observed as a weak shoulder of the main resonance. As discussed in chapter II.I, the fast decay of the intermediate state the state can hardly be observed, even though it has the same dipole strength as f 1 .
These nonclassical bandwidth features of entangled photon pairs may be further exploited to control vibronic states, as reported in (Oka, 2011a,b). They were also investigated in semiconductor quantum wells. In this system, the absorption of excited state states competes with the absorption into a continuum of excited electronic states of free electrons and holes. By interpolating between negative and positive frequency correlations (see Fig. 3), it was shown in (Salazar et al., 2012) how the absorption of the excited states may be either enhanced or suppressed, as shown in Fig. 14.

Multidimensional signals
In Fig. 2, we had explained how each light-matter interaction event also imprints the light phase φ onto the matter response. We now exploit this fact to create multidimensional spectroscopic signals. Phase-cycling is essential for partially non-collinear or collinear geometry 2D spectroscopic experiments (Keusters et al., 1999;Krčmář et al., 2013;Scheurer and Mukamel, 2001;Tan, 2008;Tian et al., 2003;Yan and Tan, 2009;Zhang et al., 2012a). In phase-cycling protocol, 0.5 0.6 0.7 Photon energy !/! g TPA (a.u.) FIG. 14 Two-photon absorption (TPIF in our notation, see section. III.B.1) in a semiconductor quantum well [taken from (Salazar et al., 2012)]: excited state resonances within a continuum of delocalized states may be enhanced with frequency anti-correlations (black) or suppressed with positive frequency correlations (red).
the desired 2D signals are retrieved by the weighted summation of data collected using different interpulse phases, φ 21 , which are cycled over 2π radians in a number of equally spaced steps. In the pump-probe configuration, Myers et al. (Myers et al., 2008) have shown that other than the pure absorptive signal, the rephasing and nonrephasing contributions may also be retrieved. This phase difference detection can be defined as a two-step phase-cycling scheme. For the phase difference detection, the Ogilvie group needs to collect signals with interpulse phases of φ 21 = 0, π, π/2 and 3π/2. The same signal can be obtained with chopping, but this is only half as intense compared to the phase difference method. In the context of this article, what they have performed is similar to a four-step phase-cycling scheme, where four sets of data need to be collected and linearly combined. Phase cycling provides a means to post-select signals with a certain phase signature, as will be shown in the next section. Using additional delay stages -like in the entangled-photon virtual state spectroscopy discussed above -these signals can be spread to create multidimensional frequency correlation maps, which carry information about couplings between different resonances or relaxation mechanisms (Ginsberg et al., 2009). (Raymer et al., 2013) had exploited the formal similarity between the TPIF and a photon coincidence signal to propose the setup shown in Fig. 15a). Here, the successive absorption events promoting the molecule into the f -state are modulated by phase cycling and delay stages in both photon beams, creating interference effects between absorption events in which the photon takes the short path through the interferometer, and those in which it takes the long path. The Fourier transform creates two-dimensional signals which are shown in panel b). In contrast to the classical signal shown in the upper panel, the proposed scheme only shows resonances at the cross-peak between the two electronic states. This could provide a useful tool to study conformations of aggregates. We next discuss this strategy in more detail.

Loop (LOP) vs ladder (LAP) delay scanning protocols for multidimensional fluorescence detected signals
Since the loop and the ladder diagrams involve different time variables they suggest different multidimensional signals obtained by scanning the corresponding time delays. We consider the TPIF signal created by a train of four pulses centered at times T 1 , T 2 , T 3 , and T 4 with phases φ 1 , φ 2 , φ 3 , and φ 4 (Pestov et al., 2009;Tekavec et al., 2007). We first analyze signals obtained by the LOP/LAP delay scanning protocols with classical light. The two protocols highlight different resonances and processes. Below we demonstrate what type of information can be extracted from each protocol for excited electronic states in a model molecular aggregate. We then show some benefits of LOP protocols when applied to entangled light.
For the system model shown in Fig. 2a) the signal (114) is given by the single loop diagram in Fig. 2b). a, b, c, d denote the pulse sequence ordered along the loop (but not in real time); a represents "first". on the loop etc. Chronologicallyordered pulses in real time will be denoted 1, 2, 3, 4 which are permutations of a, b, c, d, as will be shown below. One can scan various delays T α − T β , α, β = 1, · · · 4 and control the phases ±φ 1 ± φ 2 ± φ 3 ± φ 4 . In standard multidimensional techniques the time variables represent the pulses as they interact with sample in chronological order (Mukamel, 1995). These are conveniently given by the ladder delays t i shown in Fig. 16c). The LAP maintains complete time-ordering of all four pulses. The arrival time of the various pulses in chronological order is T 1 < T 2 < T 3 < T 4 . The ladder delays are defined as t 1 = T 2 − T 1 , t 2 = T 3 − T 2 , t 3 = T 4 − T 3 . One can then use phase cycling to select the rephasing and nonrephasing diagrams shown in Fig. 16c, and read off the signal in a sum-over-states expansion The LAP signals (117) -(118) factorize into a product of several Green's functions with uncoupled time arguments. This implies that the corresponding frequency domain signal will also factorize into a product of individual Green's functions, each depending on a single frequency argumentΩ j , j = 1, 2, 3 which yields . (120) This factorization holds only in the absence of additional correlating mechanisms the frequency variables caused by e.g. dephasing (bath) or the state of light.
In the LOP the time ordering of pulses is maintained only on each branch of the loop but not between branches. To realize the LOP experimentally the indices 1 to 4 are assigned as follows: first by phase cycling we select a signal with phase φ 1 + φ 2 − φ 3 − φ 4 . The two pulses with positive phase detection are thus denoted 1, 2 and with negative phase -3, 4. In the 1, 2 pair pulse 1 comes first. In the 3, 4 pair pulse 4 comes first. The time variables in Fig. 16b) With this choice s 1 and s 3 are positive whereas s 2 can be either positive or negative. This completely defines the LOP experimentally. When the electronic system is coupled to a bath, it cannot be described by a wave function in the reduced space where the bath is eliminated. As described in section II.J.2, the loop diagram must then be broken into several ladder diagrams shown in Fig. 2b) which represent the density matrix. The full set of diagrams and corresponding signal expressions are given in (Dorfman and Mukamel, 2014b). In Fig. 16c) we present simplified expressions for the rephasing k II and nonrephasing k III signals in the limit of well separated pulses: where I denotes the imaginary part and G αβ (t) = (−i/ )θ(t)e −[iω αβ +γ αβ ]t is the Liouville space Green's function [see Eq. (16)]. Note that the loop delays s j , j = 1, 2, 3 are coupled and enter e.g. in the Green's function G eg in Eq.
(121). Due to the Heaviside-theta function in this Green's function, the delays s j are not independent but rather couple the dynamics of the system during these delay times, which eventually result in cross-resonances in multidimensional spectra. To see the effect on the mixing of the frequency variables that yield these cross-peaks we take a Fourier transform of Eqs. (121) -(122) with respect to loop delay variable s 1 and s 3 keeping s 2 = 0 and obtain the resonant component of the signal analogous to the frequency-domain LAP signals (119) -(120) where R denotes real part. Eqs. (123)  in the frequency mixing of arguments. Fig. 17 compares the LOP and LAP signals for the model dimer parameters of (Raymer et al., 2013) under two photon excitation by classical light. The LOP spectra for rephasing k II = − 1 + 2 + 3 , nonrephasing k III = + 1 + 2 − 3 and their sum are shown in Fig. 17a)-c), respectively. The corresponding LAP spectra are shown in Fig. 17d-f. We see that the scanning protocol makes a significant difference as seen by the two columns. The LOP resonances are narrow and clearly resolve the e 1 and e 2 states whereas the corresponding LAP signals are broad and featureless. This is a consequence of the display variables chosen in each protocols. LOP variables are coupled in a very specific fashion that allows to extract the intraband dephasing rate γ ee in (Ω 1 , Ω 3 ) plot with higher precision compared to LAP case. Of course, one can extract similar information using LAP if displayed using (Ω 1 ,Ω 2 ) or (Ω 2 ,Ω 3 ) which we will discuss below in the context of entangled light. This difference in two protocols has been originally attributed to entanglement by (Raymer et al., 2013). However, it was explicitly demonstrated in (Dorfman and Mukamel, 2014b) supported by Fig. 17 that they calcu-lated entangled signals using LOP whereas they used LAP for classical signals.
a. Entangled vs classical light So far, we presented the LOP and the LAP delay scanning protocols for classical light. We now turn to the LOP protocol with entangled light. Similarly to Eqs. (121) -(122) we calculate the signals for the CWpump model (51) discussed earlier, and obtain where j = k II , k III and the display function for both rephasing and nonrephasing contributions reads and matter responses are Similarly we obtain for LAP signals Eqs. (117) where the display function for the rephasing signal is and for nonrephasing signal The complete set of signals along with frequency domain signals for entangled light can be found in Section 2 of (Dorfman and Mukamel, 2014b). The Fourier transform of the signal (125) was simulated using the LOP protocol and compared it with the standard fully time ordered LAP protocol given by Eq. (128) for a model trimer with parameters discussed in Section 5 of (Dorfman and Mukamel, 2014b). Fig. 18 shows the simulated S LOP (Ω 1 , τ 2 = 0, Ω 3 ) for a trimer using classical light (top row) and entangled light (bottom row). Fig. 18a shows classical light which gives a diagonal cross peak e = e and one pair of weak side peaks parallel to the main diagonal at (e, e ) = (e 2 , e 3 ). The remaining two pairs of side peaks at (e, e ) = (e 1 , e 2 ) and (e, e ) = (e 1 , e 3 ) are too weak to be seen. Fig. 18d) depicts the signal obtained using entangled photons where we observe two additional strong side cross peak pairs with (e, e ) = (e 1 , e 3 ) and (e, e ) = (e 1 , e 2 ). The weak peak at (e, e ) = (e 2 , e 3 ) is significantly enhanced as well.
As has been shown in Fig. 17 the LAP signal displayed vs (Ω 1 ,Ω 3 ) does not effectively reveal intraband dephasing. This can, however, be done by the LAP signal displayed vs (Ω 2 ,Ω 3 ). Fig. 18c) reveals severalΩ 3 = ω f e peaks that overlap inΩ 3 axes due to large dephasing γ f e . With entangled light, the LAP protocol yields some enhancement in several peaks around ±0.01 eV but overall the LOP yields much cleaner result. The pathways for the density matrix (LAP) and the wavefunction (LOP) are different and suggest different types of resonances.

C. Heterodyne detection of nonlinear signals
So far we have considered fluorescence (homodyne) detection. Homodyne and heterodyne are two complementary detection schemes for nonlinear optical signals. In terms of classical fields, if the sample radiation is detected for different modes than the incident radiation, the signal is proportional to |E| 2 . This is known as homodyne detection whereby the intensity is the square modulus of the emitted field itself. If the emitted field coincides with a frequency of the incident radiation E in , then the signal intensity is proportional to |E + E in | 2 . Consequently, the detected intensity contains a mixed interference term, i.e., E * E in + c.c.. This is defined as the heterodyne signal, since the emitted field is mixed with another field. In terms of quantized fields, the signal denoted homodyne if detected at a field mode that is initially vacant and heterodyne when detected at a field mode that is already occupied. Note that the current definition, which is commonly used in nonlinear multidimensional spectroscopy, is different from the definition used in quantum optics and optical engineering where homodyne and heterodyne refer to mixing with a field with the same or different frequency of the signal. In this review, we will use the spectroscopy terminology (Potma and Mukamel, 2013). Fluorescence detection is often more sensitive than heterodyne as the latter is limited by the pulse duration so there are fewer constraints on the laser system. In addition the low intensity requirements for biological samples limit the range of heterodyne detection setups. This have been demonstrated by (Fu et al., 2007;Tekavec et al., 2007;Ye et al., 2009) even in single molecule spectroscopy (Brinks et al., 2010). Historically, Ramsey fringes constitute the first example of incoherent detection (Cohen-Tannoudji and Guéry-Odelin, 2011;Ramsey, 1950;. Information similar to coherent spectroscopy can be extracted from the parametric dependence on various pulse sequences applied prior to the incoherent detection (Mukamel and Richter, 2011a;Rahav and Mukamel, 2010). Possible in- coherent detection modes include fluorescence (Barkai, 2008;Elf et al., 2007), photoaccoustic (Patel and Tam, 1981), AFM (Mamin et al., 2005;Poggio et al., 2009;Rajapaksa and Kumar Wickramasinghe, 2011;Rajapaksa et al., 2010) or photocurrent detection (Cheng and Xie, 2012;Nardin et al., 2013).

Heterodyne Intensity measurements -Raman vs. TPA pathways
With A(t) = E † (t)E(t), Eq. (111) yields the rate of change of the transmitted photon number, The semiclassical signal may be obtained from Eq. (131) by simply replacing the field operator E † (t) by a classical field amplitude A * (t). Similarly, by spectrally dispersing the timeintegrated intensity, which amounts to measuring A(ω) = E † (ω)E(ω), we obtain The third-order contribution of the time-integrated absorp-tion signal (131) is given by the four loop diagrams in Fig. 19, An identical expansion of Eq. (132) yields the frequencyresolved third-order signal (its sum-over-state expansion is shown in Appendix D) It is instructive to relate the four diagrams corresponding to Eqs. (D1)-(D4) to the transition amplitudes between the initial and final matter states to gain some intuition for this signal. This is only possible by the total photon counting signal, dω S(ω; Γ)/(2π), which represents the full energy exchanged between the light field and the matter system. We now define the transition amplitude operators When we assume unitary time evolution, we can replace the dephasing rates in Eqs. (D1)-(D4) by infinitesimal imaginary factors γ → . This allows us to use the identity 1/(ω + i ) = PP1/ω +iπδ(ω), and carry out the remaining frequency integrals. We arrive at The details of sum-over-state expansion are presented in Appendix D. These results clarify the statement we made at the end of section III.B.1: The χ (3) -absorption signal comprises matter transitions from the ground to the f -state -just like the TPIF signal -, but it also contains transitions to eand g-states.
The absorption measurement and the TPIF signal contain the same information, only when transitions to the eand g-states can be neglected, as is the case when e is off-resonant.
2. Heterodyne detected four-wave mixing; the double-quantum-coherence technique Time domain two dimensional (2D) spectroscopic techniques (Mukamel, 2000), provide a versatile tool for exploring the properties of molecular systems, such as photosynthetic aggregates (Abramavicius et al., 2008;Engel et al., 2007) or coupled (Hybrid-)nanostructures to semiconductor quantum wells (Pasenow et al., 2008;Vogel et al., 2009;Zhang et al., 2007). These techniques use sequences of short coherent pulses in comparison to the dephasing times of the system.
Earlier we had presented different delay scanning protocols (LOP and LAP) for multidimensional spectroscopy with entangled photons. These protocols can be experimentally realized using entangled photon pulse shaping, using collinear geometry and precise control of the phase, phase cycling. These protocols allow to extract inter-and intraband dephasing with high resolution by exciting doubly excited state distributions. We now present a different multidimensional timedomain spectroscopy that involve higher electronic states manifold, but does not require excited populations. This is so called double-quantum-coherence (DQC) signals (Kim et al., 2009;Mukamel et al., 2007;Palmieri et al., 2009; where the system evolves in the coherence between ground and doubly excited state |f g| rather than in a population |f f |. This technique renders the energies of single and doubly excited state energies accessible and reveals the correlations between single and doubly excited states . We show how pulsed entangled photons affect the two photon resonances. Some bandwidth limitations of classical beams are removed and selectivity of quantum pathways is possible. a. The DQC signal In the following we use the LAP delay scanning protocol which monitors the density matrix. A pulse shaper creates a sequence of four well separated chronologically ordered pulses described by field operator E j (t) = dω 2π E j (ω)e −iω(t−Tj ) , j = 1, 2, 3, 4. The control parameters are their central times T 1 < T 2 < T 3 < T 4 and phases φ 1 , φ 2 , φ 3 , and φ 4 (see Fig. 20a). The DQC signal selects those contributions with the phase signature φ 1 + φ 2 − φ 3 − φ 4 . The signal is defined as the change in the time-integrated transmitted intensity in component φ 4 , which is given by We thus have a similar configuration to an impulsive experiment with four short well separated classical fields. Introducing the LAP delays t 3 = T 4 − T 3 , t 2 = T 3 − T 2 and t 1 = T 2 − T 1 we can calculate Eq. (150) by expanding perturbatively in the dipole field-matter interaction Hamiltonian (7). The two contributions to the signal are represented by the ladder diagrams shown in Fig. 20b. The corresponding signal (150) can be read off these diagrams and is given by S (LAP ) We have introduced the complex frequency variables ξ ij = ω ij + iγ ij , where ω ij = ε i − ε j are the transition frequencies and γ ij are the dephasing rates. The signal may be depicted by its variation with various parameters of the field wavefunction. These are denoted collectively as Γ. Various choices of Γ lead to different types of 2D signals. These will be specifiedbelow (see Eqs. (162) and (163)).
b. The field correlation function for entangled photon pairs We consider the pulsed entangled photon pairs described in section II.F. Specifically, we here consider the wavefunction introduced in (Keller and Rubin, 1997), where the two-photon wavefunction Φ( i , k j ) takes the form where A P (· · · ) is the pulse envelope and Ω p the central frequency of the pump pulse used to generate the pairs. The correlation functions of the entangled fields reads (Keller and Rubin, 1997) where V 0 is given by ... It describes the pulsed counterpart of the cw-correlation function (52) we had described earlier, where the step-function rect (t) is amended by the finite pulse amplitude A p .
c. Simulated 2D Signals Below we present simulated DQC signal for a model trimer system with parameters given in Ref. (Richter and Mukamel, 2010). By inserting Eq. (155) in Eq. (151) -(152) and carrying out all integrations, assuming that γ σ, we arrive at the final expression for the two contribu-tions to the signal: The control parameters Γ now include the delay times (t 1 , t 2 , t 3 ) and the entanglement time T . For comparison, we present the same signal obtained with four impulsive classical pulses with envelopes E j (·) j = 1, . . . , 4 E 2 (·), E 3 (·) and E 4 (·) and carrier frequency Ω 0 p : With Eqs. (157) and (158), we can compare the entangled photon and classical DQC signals. We first note that Eq. (157) scales with the intensity of the generating pump pulse, in contrast with the intensity square of the classical case Eq. (160).
In the classical case the signal is limited by the bandwidths of the four pulses (cf. Eqs. (160-161)), which control the four transitions (ω eg , ω f e , ω e g , ω f e ) in the two photon transitions inside the pulse bandwidth ). In Eq.(157) bandwidth limitations of the envelopes are only imposed through the bandwidth of entangled photon pair A p (·) and the limitation is only imposed on the two photon transition ω f g , leading to a much broader bandwidth for the ω eg and ω f e transitions, if the ω f g transition is within the generating pump pulse bandwidth.
We illustrate this effect in Fig. 20a for the following 2D signal: ω eg resonances are seen alongΩ 1 and ω eg and ω f e on axis Ω 3 . As the bandwidth is reduced, we only get contributions from the doubly excited state resonant to the generating pump pulse. This results in four identical patterns along the Ω 1 axis. All peaks are connected to the same doubly excited states. More precisely we get four contributions along theΩ 1 axis connected to the transitions ω f3e1 , ω f3e2 overlapping with ω e3g , ω e2g overlapping with ω f3e3 and ω e1g . The remaining transitions are not affected by the reduced pump bandwidth.
Here the narrow bandwidth of the pump can be used to select contributions in the spectra connected to a specific doubly excited state. In Fig. 20b we display a different signal: ω eg resonances are now seen alongΩ 1 and ω f g inΩ 2 . This is similar to Fig. 20f except that here we see the singly excited state contributions ω e1g ,ω e2g , ω e3g to the selected doubly excited state f 3 along a single row.
Bandwidth limitations on the of the singly excited state transitions ω eg and ω f e are only imposed indirectly by the factors in Eqs. (157-158) which depend on the entanglement time T . These become largest for ω f g = 2ω f e .
The factors in Eq.(157-158) which depend on the entanglement time T contain an interference term of the form (e i(ω−γ)T − 1), where ω is a material frequency (see sections II.J and III.B.2). If we now vary the entanglement time, some resonances will interfere destructively for values of the entanglement time which match the period of ω. One can therefore use the entanglement times to control selected resonances. This holds only as long the entanglement times are not much bigger than the dephasing time, since in this case the signal will be weak. The frequencies ω can be in the contributing diagrams ω f g /2 − ω eg or ω f g /2 − ω f e (which is differs from the first frequency only by a sign) for different combination of the states e, e and f . By varying T , we expect an oscillation of the magnitudes of resonances with different frequencies.The details of manipulation of entanglement time as a control parameters have been studied in Refs. (Richter and Mukamel, 2010; Previously, to utilize entangled photons for DQC signal we used two pulsed entangled photon pairs (Richter and Mukamel, 2010) (k 1 , k 2 ) and (k 3 , k 4 ). With the pulse shaping described in section II.G, a single shaped entangled photon pair is sufficient to realize any time-domain four-wave mixing signals.

D. Multiple photon counting detection
Another class of multidimensional signals is possible by detecting sequences of individual photons emitted by an optically driven system. With proper gating, each photon j can be characterized by a frequency ω j and a time t j . By detecting N photons we thus obtain a 2N dimensional signal parametrized by ω 1 , t 1 , ..., ω N , t N . Unlike coherent multidimensional heterodyne signals which are parametrized by delays of the incoming fields, these incoherent signals are parametrized by the emitted photons. The photon time t j and frequency ω j are not independent and can only be detected to within a Fourier uncertainty ∆ω j ∆t j > 1. This poses a fundamental limit on the temporal and spectral resolution. We derive these signals and connect them to multipoint 2N dipole correlation functions of matter. The Fourier uncertainty is naturally built in by a proper description of photon gating and need not be imposed in an ad hoc matter as is commonly done.
A semiclassical formalism for photon counting was first derived by (Mandel, 1958(Mandel, , 1959. The full quantum mechanical description of the field and photon detection was developed by Glauber (Glauber, 2007b). The theory of the electromagnetic field measurement through photoionization and the resulting photoelectron counting has been developed by Kelley and Kleiner (Kelley and Kleiner, 1964). The experimental application to normal and time ordered intensity correlation measurements was given in the seminal work of Kimble et al. (Kimble et al., 1977). According to these treatments free-field operators, in general, do not commute with sourcequantity operators. This is the origin of the fact that the normal and time ordering of the measured field correlations, according to the Kelley-Kleiner theory (Kelley and Kleiner, 1964), are transformed into normal and time ordered source quantities occurring inside the integral representations of the filtered source-field operators. This constitutes the back reaction of the detector on the field state (Cohen-Tannoudji et al., 1992). An ideal photon detector is a device that measures the radiation field intensity at a single point in space. The detector size should be much smaller than spatial variations of the field. The response of an ideal time-domain photon detector is independent of the frequency of the radiation.
The resolution of simultaneous frequency and time domain measurements is limited by the Fourier uncertainty ∆ω∆t > 1. A naive calculation of signals without proper time and frequency gating can work for slowly varying spectrally broad optical fields but otherwise may yield unphysical negative signals (Eberly and Wodkiewicz, 1977). In Ref. (Mukamel et al., 1996), the mixed time-frequency representation for the coherent optical measurements with interferometric or autocorrelation detection were calculated in terms of a mixed material response functions and a Wigner distribution for the incoming pulses, the detected field and the gating device. Multidimensional gated fluorescence signals for single-molecule spectroscopy have been calculated in Ref. (Mukamel and Richter, 2011b).
The standard Glauber's theory of photon counting and correlation measurements (Glauber, 2007a;Mollow, 1972;Scully and Zubairy, 1997) is formulated solely in the radiation field space (matter is not considered explicitly). Signals are related to the multi-point normally-ordered field correlation function, convoluted with time and frequency gating spectrograms of the corresponding detectors. This approach assumes that the detected field is given. Thus, it does not address the matter information and the way this field has been generated. Temporally and spectrally resolved measurements can reveal important matter information. Recent single photon spectroscopy of single molecules (Fleury et al., 2000;Lettow et al., 2010;Rezus et al., 2012) call for an adequate microscopic foundation where joint matter and field information could be retrieved by a proper description of the detection process.
A microscopic diagrammatic approach may be used for calculating time-and-frequency gated photon counting measurements (Dorfman and Mukamel, 2012a). The observed signal can be represented by a convolution of the bare signal and a detector spectrogram that contains the time and frequency gate functions. The bare signal is given by the product of two transition amplitude superoperators ) (one for bra and one for ket of the matter plus field joint density matrix), each creating a coherence in the field between states with zero and one photon. By combining the transition amplitude superoperators from both branches of the loop diagram we obtain the photon occupation number that can be detected. The detection process is described in the joint field and matter space by a sum over pathways each involving a pair of quantum modes with different time orderings. The signal is recast using time ordered superoperator products of matter and field. In contrast to the Glauber theory that uses normally ordered field operator, the microscopic approach of (Dorfman and Mukamel, 2012a) is based on time-ordered superoperators. Ideal frequency domain detection only requires a single mode (Mukamel and Richter, 2011b). However, maintaining any time resolution requires a superposition of several field modes that contain the pathway information. This information is not directly accessible in the standard detection theory that operates in the field space alone (Glauber, 2007b).

Photon correlation measurements using gated photon number operators
Time-and-frequency gated N -th order photon correlation measurement performed at N detectors characterized by central time t j and central frequency ω j , j = 1, ..., N is defined as where ... T = Tr[...ρ T (t)] and ρ T (t) represents the total density matrix of the entire system in joint field plus matter space and contains information about system evolution prior to the detection (e.g. photon generation process, etc.). Γ j , j = 1, ...N represents other parameters of the detectors such as e.g. bandwidth (σ j T and σ j ω are the time gate, and frequency gate bandwidths, respectively). The time-and-frequency gated photon number superoperator is given bŷ n t,ω = dt dτ D(t, ω; t , τ )n(t , τ ).
Here D(t, ω, t , τ ) is a detector time-domain spectrogram (the ordinary function, not an operator) which takes into account the detector's parameters, and is given by where F t and F f are time and frequency gating functions that are characterized by central time t and frequency ω and detection bandwidths σ T and σ ω , respectively. Note that in Eq. (166) the time gate is applied first, followed by the frequency gate. Similar expression can be written if the order in gating is reversed.n(t, t ) is a bare photon number superoperator defined in terms of the bare field operators aŝ Details may be found in Appendix E.

Photon counting and matter dipole correlation functions
To connect the photon coincidence counting (PCC) signals to matter properties one needs to expand the density operator in Eq. (167) in perturbative series over field-matter interactions. We first calculate the time-and-frequency resolved emission spectra n t,ω = dt dτ D(t, ω; t , τ )n(t , τ ), where the bare photon number n(t , τ ) ≡ Tn(t , τ ) T is an expectation value of the bare photon number operator with respect to the total density matrix. The leading contribution is coming from the second-order in field matter interactions with vacuum modes [see diagram in Fig. 21b] where D(ω) = 1 2πD (ω) is a combined density of states evaluated at the central frequency of the detector ω for smooth enough distribution of modes. The corresponding detected signal (165) is given by One can similarly calculate the second-order bare correlation function Tn t1,ω1nt2,ω2 T = dt 1 dτ 1 D (1) (t 1 ω 1 ; t 1 , τ 1 ) The bare PCC rate Tn(t 1 , τ 1 )n(t 2 , τ 2 ) T can be read off the diagram shown in Fig. 21c. The leading contribution requires a fourth-order expansion in field-matter interactions After tracing over the vacuum modes we obtain The gated coincidence signal (172) is finally given by Therefore, the fundamental material quantity that yields the emission spectra (171) is a two-point dipole correlation function given in Eq. (170), and for the coincidence g (2)measurement (175) it is the four-point dipole correlation function in Eq. (174). Sum-over-state expansion of these expressions is given in appendix E.

Connection to the physical spectrum
An early work (Eberly and Wodkiewicz, 1977), had argued that detector gating with finite bandwidth must be added to describe the real detector. In recent work (González-Tudela et al., 2015; del Valle et al., 2012) used a two-level model detector with a single parameter Γ that characterize both timeand-frequency detection is comsidered. This was denoted the physical spectrum (178) which we shall derive in the following. It can be recovered from our model by removing the time gate F t = 1 and using a Lorentzian frequency gate Using the physical spectrum, the time-and-frequency resolved photon coincidence signal is given by is the gated field. This model provides a simple benchmark for finite-band detection, which has several limitations. First, the time and frequency gating parameters are not independent, unlike the actual experimental setup, where frequency filters and avalanche photodiodes are two independent devices. Second, this method does not address the generation and photon bandwidth coming from the emitter, as the analysis is performed solely in the field space. Finally, the multi-photon correlation function presented in (González-Tudela et al., 2015) is stationary. For instance, the four-point bare correlation function  (179) is recast using C B (ω 1 , ω 2 , t 2 − t 1 ), which only depends on the time difference t 2 − t 1 , which is an approximation for stationary fields. This model also works if t Γ −1 , which means that Γ cannot approach zero (perfect reflection in Fabri Perot cavity). It also works when Γτ 0 1 where τ 0 is the scale of change in the field envelope. For comparison, the photon coincidence counting (PCC) (164) for N = 2 reads g (2) (t 1 , ω 1 , Γ 1 ; t 2 ; ω 2 , Γ 2 ) Tn t1,ω1nt2,ω2 Tn t1,ω1 Tn t2,ω2 , which depends on two time t 1 , t 2 and two frequency ω 1 , ω 2 arguments. The theory summarized above which gives rise to Eq. (180) has several advantages compared to the physical spectrum used by (González-Tudela et al., 2015;del Valle et al., 2012). First, independent control of time and frequency gates (with guaranteed Fourier uncertainty for the time and frequency resolution) along with the fact that the bare photon number operator depends on two time variablesn(t, τ ) allows to capture any dynamical process down to the very short scale dynamics in ultrafast spectroscopy applications. Second, the gating (165) provides a unique tool that can capture nonequlibrium and non-stationary states of matter which can be controlled by gating bandwidths. In this case a series of frequency ω 1 , ω 2 correlation plots (keeping the central frequencies of the spectral gates as variables) for different time delays t 1 − t 2 yields a 2D spectroscopy tool capable of measuring ultrafast dynamics. Third, the superoperator algebra allows to connect the gated field correlation function with the bare correlation function (179), with time-andfrequency gates (arbitrary, not necessarily Lorentzian) as well as material response that precedes the emission and detection of photons. The superoperator expressions require timeordering, and can be generalized to correlation functions of field operators that are not normally ordered. Superoperators provide an effective bookkeeping tool for field-matter interactions prior to the spontaneous emission of photons. We next apply it to the detection of photon correlations. Finally, as we show in the next section, PCC can be recast in terms of matter correlation functions by expanding the total density matrix operator in a perturbation series, and tracing the vacuum modes. This way, photon counting measurements can be related to the matter response which is the standard building block of nonlinear spectroscopy.

E. Interferometric detection of photon coincidence signals
Photon coincidence signals also known as biphoton signals (Kalachev et al., 2007;Scarcelli et al., 2003;Slattery et al., 2013;Yabushita and Kobayashi, 2004) became recently available as a tool for nonlinear spectroscopy. In a typical setup, a pair of entangled photons denoted as E s , and E r generated by PDC are split on a beam splitter [see Fig. 22a)]. One photon E s is transmitted through the molecular sample and then detected in coincidence with E r . In order to use it as a spectroscopic tool, a frequency filter can be placed in front of one of the detectors which measures the spectrum. This type of signal shows a number of interesting features: First, coincidence detection improves the signal-to-noise ratio (Kalashnikov et al., 2014). Second, it is possible to make the two detectors operate in very different spectral regions and at different spatial locations (Kalachev et al., 2007). For example, to measure the spectroscopic properties of a sample in vacuum ultraviolet (VUV) range, it is not necessary to set a spectrometer in a vacuum chamber and control it under the vacuum condition. Instead, provided that a entangled photon pair consists of a VUV and a visible photon, only the latter is to be resolved by a spectrometer under usual conditions. Another advantage is when spectroscopic measurements are to be performed in infrared range. The power of the light source must often be very low to prevent possible damage of a sample, but an infrared photodetector is usually noisy. Photon coincidence measurements involve the lowest intensities of lightsingle photons, and can overcome the noise.
In the following, we present photon coincidence version of three signals: linear absorption, pump-probe, and Femtosecond Stimulated Raman Signals (FSRS).  (Kalachev et al., 2007) with coincidence detection. The entangled photon pair in beams Es and Er are split on a beam splitter. Es is employed as a probe transmitted through the sample, while Er is detected in coincidence -(a).The absorption spectrum of Er 3+ ion in YAG crystal around 650 nm obtained by single photon counting (blue filled squares, left Y-axis) and coincidence counting (red open squares, right Y-axis) at various values of signal-to-noise ratio in the channel with the sample: 1/2 (b), 1/5 (c), and 1/30 (d)

Coincidence detection of linear absorption
Linear absorption is the most elementary spectroscopic measurement. Combined with photon coincidence detection [see Fig. 22a)], it yields interesting results: Below we present the simplest intuitive phenomenological approach. In following sections we present more rigorous microscopic derivation for pump-probe and Raman signals.
In the case of the linear signal, the joint detection of two entangled photons provides linear absorption information if one of the photons is transmitted through a molecular sample. If the coincidence gate window accepts counts for a time T , then the joint detection counting rate, R c , between detectors r and s is proportional to where ψ is a two-photon entangled state (33). The gating window T is typically much larger than the reciprocal bandwidth of the light or the expected dispersive broadening Denoting the spectral transfer functions of the sample and monochromator H S (ω) and H M (ω), respectively. In this case, and provided that the signal and idler photons are separated by a 50/50 beam splitter. For narrowband down conversion Φ(ω 1 , ω 2 ) = F ω1−ω2 2 δ(ω 1 + ω 2 − 2ω 0 ). The coincidence counting rate is then given by Now, assume that H M (ω) is much narrower than H S (ω) and Φ(ω 1 , ω 2 ), we can set H M (ω) = δ(ω − ω M ), and the frequency ω M does not exceed the frequency range in which function F (ω) is essentially nonzero. Dividing the coincidence counting rate with a sample, R c,sample , by one without the sample, R c , we obtain the absorption spectrum where subscript ILA marks the interferometric nature of the photon coincidence detection combined with linear absorption measurement. Thus, the joint detection counting rate reproduces the spectral function of the sample, which is reversed in frequency with respect of the pump frequency, provided that the line-width of the pumping field as well as bandwidth of the monochromator are narrow enough for resolving the absorption spectrum features.
In (Kalachev et al., 2007) the coincidence signals were used to measure the spectroscopic properties of YAG:Er 3+ crystal. In order to demonstrate the advantage of coincidence detection in the presence of an enhanced background noise. Figs. 22b)-d) show the central part of the absorption spectrum measured in two ways: using the coincidence counting as described before, and using the single-photon counting, when the sample was placed above the monochromator in the idler channel. The signal-to-noise ratio in the channel was changed from 1/2 to 1/30 in both cases. It is evident from these experimental data that the standard classical method which suffers from a high noise level does not allow one to obtain any spectroscopic information, but the coincidence counting measurement in contrast does not undergo the reduction in resolution and is resistant to noise. The example of YAG:Er 3+ has been later extended to plasmonic nanostructures (Kalashnikov et al., 2014).

Coincidence detection of pump-probe signals
Consider the setup depicted in Fig. 23a). Unlike the linear absorption case here the probe photon is sent through the sample, which has previously been excited by a classical ultrafast laser pulse, and then detected .
The interferometric pump probe (IPP) signal with photon coincidence detection is dispersed spectrally by placing spectral filters in front of both detectors, and our signal is given by the change to this two-photon counting rate that is governed by a four-point correlation function of the field (Cho et al., 2014a;Mosley et al., 2008) Here, ω/ω r denotes the detected frequency of the respective spectral filter, and the brackets · · · represent the expectation value with respect to the transmitted fields. To obtain the desired pump-probe signal in a three-level system, we used a third order perturbation theory . The first two interactions are with the classical pump pulse, which is taken to be impulsive, E a (t) = E a δ(t), and the third is with the probe centered at t = t 0 . Assuming E s to be far offresonant from the e − g transition, we obtain only the single diagram shown in Fig. 23c), which reads We have defined the matter correlation function, where µ ge and µ ef denote the dipole moments connecting ground state with the singly excited states manifold, as well as single with doubly excited state manifold, respectively. · · · env denotes the average with respect to environmental degrees of freedom, obtained from tracing out the bath. Here, we employ a stochastic Liouville equation (Tanimura, 2006) which represents the two-state jump (TSJ) model: A ground state g is dipole-coupled to an electronic excited state e, which is connected to two spin states ↑ and ↓ undergoing relaxation (Šanda and . We additionally consider a doubly excited state f , which is dipole-coupled to both |e, ↑ and |e, ↓ [see Fig. 23b)]. The electronic states are damped by a dephasing rate γ. We assume the low-temperature limit, where only the decay from ↑ to ↓ is allowed (Dorfman et al., 2013b). The decay is entirely incoherent, such that the description may be restricted to the two spin populations | ↑ ↑ |=(1, 0) T and | ↓ ↓ |=(0, 1) T . The field correlation function (190) is then given by where δ the energy difference between the two spin states, and ω ± = ω f e ±δ. Note that, since we monitor the f −e transition, the detected frequency will increase in time, from ω − to ω + . In the following, we use these results to first simulate the classical pump-probe signal, and then the two-photon counting signal with entangled photons. For the former case, we consider a classical Gaussian probe pulse We chose the following system parameters: ω f e = 11, 000 cm −1 , δ = 200 cm −1 , k = 120 cm −1 , and γ = 100 cm −1 . The center frequency ω 0 is fixed at the transition frequency ω f e , and we vary the probe bandwidth. Panel a) shows the signal for σ = 1, 000 cm −1 . Two peaks at ω f e ± δ correspond to the detected frequency, when the system is either in the upper state (at ω f e − δ), or in the lower state (ω f e + δ). Due to the spectrally dispersed detection of the signal, the resonance widths are given by the linewidth γ, and not the much broader probe pulse bandwidth σ. For very short time delays t 0 , both resonances increase, until the probe pulse has fully passed through the sample. Then the resonance at 10, 800 cm −1 , i.e. the state |e, ↑ , starts to decay rapidly, while the resonance at 11, 200 cm −1 peaks at longer delay times due to its initial population by the upper state. For longer delays, both resonance decay due to the additional dephasing.
The two-photon counting signal with entangled photons offers novel control parameters: The dispersed frequency ω of beam 1, the pump frequency ω p and its bandwidth σ p loosely correspond to the classical control parameters, i.e. the central frequency ω 0 and bandwidth σ. In addition, we may vary the entanglement time T and the detected frequency of the reference beam ω r . Fig. 24 depicts the signal (189) obtained from entangled photons with T = 90 fs for different time delays t 0 . For comparison, we show the classical pump-probe with bandwidth σ = 1, 000 cm −1 in the left column, and with 100 cm −1 in the right column. The TPC signals are normalized with respect to the maximum value of the signal at t 0 = 3 fs and FIG. 23 (Color online) a) The IPP setup: the entangled photon pair in beams Es and Er are split on a beam splitter. A classical, actinic, ultrafast pulse Ea excites the sample, and Es is employed as a probe in a pump-probe measurement, while Er is detected in coincidence. b) The level scheme for the two-state jump (TSJ) model considered in this work. The g − e transition is far off-resonant from the spectral range of the entangled photon wavepacket, which only couples to the e − f transition. c) Diagram representing the pump-probe measurement. [taken from ] ω r = 11, 400 cm −1 . The classical signal is normalized to its peak value at zero time delay, and the TPC signals to the signal with ω r = 11, 400 cm −1 at zero delay. As becomes apparent from the figure, a broadband classical probe pulse (left column) cannot excite specific states, such that the two resonances merge into one band. A narrowband probe (right column), on the other hand, cannot resolve the fast relaxation at all, and only shows the unperturbed resonance at ω f e . Interferometric signals, however, can target the relaxation dynamics of the individual states.

Coincidence detection of Femtosecond Stimulated Raman Signals
So far, we have demonstrated how coincidence detection can enhance linear absorption and pump-probe signals. We now demonstrate the power of this interferometric detection for stimulated Raman signals commonly used to probe molecular vibrations. Applications include probing time-resolved photophysical and photochemical processes (Adamczyk et al., 2009;Kukura et al., 2007;Kuramochi et al., 2012;Schreier et al., 2007), chemically specific biomedical imaging (Cheng et al., 2002), remote sensing (Arora et al., 2012;Pestov et al., 2008). Considerable effort has been devoted to increasing the sensitivity and eliminating off-resonant background, thus improving the signal-to-noise ratio and enabling the detection of small samples and even single molecules. Pulse shaping (Oron et al., 2002;Pestov et al., 2007) and the combination of broad and narrow band pulses (technique known as Femtosecond Stimulated Raman Spectroscopy (FSRS) (Dietze and Mathies, 2016)) have been employed. Here, we present an Interferometric FSRS (IFSRS) technique that combines quantum entangled light with interferometric detection (Kalachev et al., 2007;Scarcelli et al., 2003;Slattery et al., 2013;Yabushita and Kobayashi, 2004) in order to enhance the resolution and selectivity of Raman signals . The measurement uses a pair of entangled photons, one (signal) photon interacts with the molecule, and acts as the broadband probe, while the other (idler) provides a reference for the coincidence measurement. By counting photons, IF-SRS can measure separately the gain and loss contributions to the Raman signal (Harbola et al., 2013) which is not possible with classical FSRS signals that only report their sum previousy (i.e. the net gain or loss). We showed how the entangled twin photon state may be used to manipulate two-photon absorption ω 1 + ω 2 type resonances in aggregates (Dorfman and Mukamel, 2014b;Lee and Goodson, 2006;Saleh et al., 1998; but these ideas do not apply to Raman ω 1 − ω 2 resonances. In FSRS, an actinic resonant pulse E a first creates a vibrational superposition state in an electronically excited state (see Fig. 25a,b). After a variable delay τ , the frequency resolved transmission of a broadband (femtosecond) probe E s in the presence of a narrowband (picosecond) pump E p shows excited state vibrational resonances generated by an off-resonant stimulated Raman process. The FSRS signal is given by (Dorfman et al., 2013a) where α is the electronic polarizability, I denotes the imaginary part, and E s = E s is expectation value of the probe field operator with respect to classical state of light (hereafter E denotes classical fields and E stands for quantum fields). H − is the Hamiltonian superoperator in the interaction picture which, for off resonance Raman processes, can be written as where V is the dipole moment, α is the off resonant polarizability. Formally, this is a six-wave mixing process. Expand- ing the signal (193) to sixth order in the fields ∼ E 2 s E 2 p E 2 a . we obtain the classical FSRS signal The two terms correspond to the two diagrams in Fig. 25c). All relevant information is contained in the two four point correlation functions where the retarded Green's function G(t) = (−i/ )θ(t)e −iHt represents forward time evolution with the free-molecule Hamiltonian H (diagrams (1, 1)a, (1, 1)b) and G † represents backward evolution. F i involves one forward and two backward evolution periods and F ii contains two forward followed by one backward propagation. F i and F ii differ by the final state of the matter (at the top of each diagram). In F i (F ii ) it is different (the same) as the state after preparation by the actinic pulse.
a. Photon correlation measurements. In IFSRS, the probe pulse E s belongs to a pair of entangled beams generated in degenerate type-II PDC. The polarizing beam splitter (BS) in Fig.  25d) then separates the orthogonally polarized photons. The horizontally polarized beam E s is propagating in the s arm of the interferometer, and interacts with the molecule. The vertically polarized beam E r propagates freely in the r arm, and serves as a reference. IFSRS has the following control knobs: the time and frequency parameters of the single photon detectors, the frequency of the narrowband classical pump pulse ω p , and the time delay T between the actinic pulse E a and the probe E s .
As discussed in section III.D.1, the joint time-andfrequency gated detection rate of N s photons (N s = 0, 1, 2) in detector s and a single photon in r when both detectors have narrow spectral gating is given by where Γ i stands for the incoming light beam parameters. In the standard Glauber's approach (Glauber, 2007a), the correlation function is calculated in the field space using normallyordered field operators. The present expressions in contrast are given in the joint field and matter degrees of freedom, and the bookkeeping of the fields is instead solely based on the time ordering of superoperators. Normal ordering is never used.
To make sure that there is no photon at detector s we had integrated over its entire bandwidth, thus eliminating the dependence on detector parameters. For the Raman gain N s = 2 signal (i.e. two photons in the s-arm, one photon in the r-arm), when both s and r detectors have narrow frequency gates, we get S (2,1) Finally, the N s = 1 signal (single photon in each arm) is given by b. Photon counting detection window for the molecular response The input two-photon state has a single photon in each the smode and single photon in the r -mode. and is described where a † ωs (a † ωr ) is the creation operator of a horizontally (vertically) polarized photon and the two-photon amplitude is given by where is the classical pump T j = (1/v p − 1/v j )L, j = 1, 2 is the time delay between the j − th entangled and the classical pump beam after propagation through the PDC crystal. T = T 2 − T 1 is the entanglement time. In A p (ω), the sumfrequency ω s + ω r is centered around 2ω 0 with bandwidth σ 0 . For a broadband classical pump, the frequency difference ω s − ω 0 becomes narrow with bandwidth T −1 j , j = 1, 2. The output state of light in mode s may contain a varying number of photons, depending on the order of the field-matter interaction.
As discussed in section II.C, the twin photon state Eq. (204) is not necessarily entangled. Using the Schmidt decomposition (34), we obtain for the two-photon amplitude in Eq. (205) the rich spectrum of eigenvalues shown in Fig. 26l) shows that the state is highly entangled with up to 20 modes contributing with appreciable weight. As will be shown later [Eq. (208)], this entanglement is reflected in the violation of the Fourier uncertainty ∆ω∆t ≥ 1 in two-photon transitions.
We next turn to how entanglement affects the Raman resonances. Both FSRS and IFSRS signals are governed by a fourpoint matter correlation function (two polarizabilities α ac and two dipole moments V ag as depicted by the loop diagrams shown in Fig. 25c) and e), respectively. Depending on the number of photons detected, this four-point matter correlation function is convoluted with different field correlation functions. For N s = 0, and N r = 1 Eq. (200) is given by the four-point correlation function for a quantum field [red arrows in Fig. 25e)]. For a twin photon state, we recall from Eq. (50) that the field correlation function can be factorized as The We next compare the different field spectrograms which represent the temporal and spectral windows created by the fields. Fig. 26i depicts the time-frequency Wigner spectrogram for the classical probe field E s : The Fourier uncertainty ∆ω∆t ≥ 1 limits the frequency resolution for a given time resolution. The corresponding Wigner two-photon spectrogram for the entangled twin photon state is depicted in Fig. 26k. For the same temporal resolution of the FSRS, the spectral resolution of IFSRS is significantly higher since the time and frequency variables for entangled light are not Fourier conjugate variables . The high spectral resolution in the entangled case is governed by T −1 j , j = 1, 2 which is narrower than the broadband probe pulse such that ∆ω∆t 0.3. Fig. 26j) demonstrates that entangled window function R (Ns,1) q for N s = 1, 2 (see Eqs. (212), (213)) that enters the IFSRS Signal (210) yields a much higher spectral resolution than the classical R c in Eq. (215).
The molecular information required for all three measurement outcomes (N s = 0, 1, 2) is given by two correlation functions F i and F ii (see Fig. 25c,), which are then convoluted with a different detection window for FSRS and IFSRS. F i and F ii may not be separately detected by FSRS. However, in IFSRS the loss S IF SRS depends on F ii (initial and final states are identical). Interferometric signals can thus separately measure the two matter correlation functions.
c. IFSRS for a vibrational mode in a tunneling system. To demonstrate the effect of entanglement in interferometric measurements, we show the calculated signals for the threelevel model system depicted in Fig. 25a. Once excited by the actinic pulse, the initial state with vibrational frequency ω + = ω ac + δ can tunnel through a barrier at a rate k and assume a different vibrational frequency ω − = ω ac − δ (see Ref. (Dorfman et al., 2013b)). The probability to be in the initial state with ω + decreases exponentially P + (t) = e −kt , whereas for ω − it grows as P − (t) = 1−e −kt . This model represents Kubo's two-state jump model in the low temperature limit (Dorfman et al., 2013b;Kubo, 1963) which we had also discussed before in section III.E.2. The absorption lineshape is given by . (209) This shows two peaks with combined width governed by dephasing γ a and tunneling rate k. The corresponding IFSRS signal S (Ns,1) IF SRS with N s = 0, 1, 2 is given by where ν = − for N s = 0, 2 and ν = + for N s = 1, µ = − for N s = 1, 2 and µ = + for N s = 0. The Raman response R (Ns,1) q which depends on the window created by the quantum field for different photon numbers N s is given by For comparison, we give the classical FSRS signal (193) where is the Raman response gated by the classical field.  Fig. 26a)] has two well-resolved peaks at ω ± . The classical FSRS shown in Fig. 26b) has one dominant resonance at ω + which decays with increase of delay T , whereas the ω − peak slowly builds up and dominates at longer T . This signal contains both blue-and red-shifted Raman resonances relative to the narrowband pump frequency: ω − ω p = ±ω ± . If the modulation and dephasing rates are comparable to the level splitting k, γ a ∼ δ, then the ω ± resonances in the absorption [ Fig. 26e)] and the classical FSRS [ Fig. 26f)] broaden, and become less resolved.
We next compare this with IFSRS. For small modulation and long dephasing, S (1,1) IF SRS is similar to the classical FSRS [see Fig. 26c)]. However, both temporal and spectral resolution remains high, even when the modulation is fast and the dephasing width is large, as seen in Fig. 26g). The same applies to the S (2,1) IF SRS signal depicted for slow - Fig. 26d) and fast - Fig. 26h) tunneling.
Apart from the different detection windows, there is another important distinction between IFSRS (Eq. (210)) and the classical FSRS (214) signals. In the latter, the gain and loss contributions both contain red-and blue-shifted features relative to the narrow pump. The FSRS signal can contain both Stokes and anti Stokes components. FSRS can only distinguish between red and blue contributions. In contrast, the interferometric signal can measure separately the gain S d. The role of entanglement. We now show that the enhanced resolution of Raman resonances may not be reproduced by classically shaped light and entanglement is essential. To this end, we calculate the IFSRS signals (210) for the correlatedseparable state (Zheng et al., 2013) described by the density matrix ρ cor = ∞ −∞ dω s dω r |Φ(ω s , ω r )| 2 |1 ωs , 1 ωr 1 ωs , 1 ωr |.

(216)
This is a diagonal part of the density matrix corresponding to state Eq. (204) with amplitude Eq. (205), which results from the disentanglement of the twin state. This state yields the same single-photon spectrum and shows strong frequency correlations similar to entangled case, and is typically used as a benchmark to quantify entanglement in quantum information processing (Law et al., 2000). We further compare this with signals from the fully separable uncorrelated Fock state given by Eq. (204) with , k = s, r with parameters matching the classical probe pulse used in FSRS. Fig. 26m-p illustrates S (1,1) IF SRS for these two states of light. The separable correlated state shown in Fig. 26m has high spectral and no temporal resolution, as expected from a cw time-averaged state in which the photons arrive at any time (Zheng et al., 2013). The separable uncorrelated state (see Fig. 26o) yields slightly better resolution than in the classical FSRS signal in Fig. 26f. Similar results can be obtained for the S (2,1) IF SRS (see Fig. 26n and p respectively).

IV. ENTANGLED LIGHT GENERATION VIA NONLINEAR LIGHT-MATTER INTERACTIONS; NONCLASSICAL RESPONSE FUNCTIONS
The light generation process is typically based on the assumption that the light field is given, or that different optical fields that participate in generation of the desired light state interact with other fields. This way the generation process can be described by an effective Hamiltonian in the field space alone. The material quantities that assist the light conversion are typically set to be constant parameters governed by nonlinear semiclassical susceptibilities.
Superoperator non-equilibrium Green's functions are useful for calculating nonlinear optical processes involving any combination of classical and quantum optical modes. Closed correlation-function expressions based on superoperator timeordering may be derived for the combined effects of causal (response) and non-causal (spontaneous fluctuations) correlation functions (Roslyak and Mukamel, 2009b).
Below we survey several wave-mixing schemes for generating quantum light by using a combination of classical and quantum modes of the radiation field. Homodyne-detected sum frequency generation (SFG) (Shen, 1989) and difference frequency generation (DFG) (Dick and Hochstrasser, 1983;Mukamel, 1995) involve two classical and one quantum mode. Parametric down conversion (PDC) Klyshko, 1988;Louisell et al., 1961;Mandel and Wolf, 1995) involves one classical and two quantum modes and is one of the primary sources of entangled photon pairs (Edamatsu, 2007;Gerry and Knight, 2005;U'Ren et al., 2006). All of these are coherent measurements, and scale as N (N − 1) for N active molecules the signals (Marx et al., 2008).
We further consider incoherent ∼ N -scaling signals. These include heterodyne detected SFG and DFG, which involve three classical modes, and two types of two photon fluorescence (Denk et al., 1990): two photon induced fluorescence with one classical and two quantum modes (TPIF) (Callis, 1993;Rehms and Callis, 1993;Xu and Webb, 1996) and two photon emitted fluorescence with two classical and one quantum make (TPEF). The list of the different measurement is summarized in the Table I. Finally we present a more detailed microscopic theory of entangled light generation in two schemes, which describe light-matter interactions that involve quantum field. The first is based on type-I PDC in a cascade three-level scheme, and in the second scheme an entangled photon pair is generated by two remote molecules assisted by an ideal 50:50 beam splitter.

A. Superoperator description of n-wave mixing
So far we mostly used L and R representation for describing signals. Here, we consider various nonlinear signals which involve two photon resonances (Roslyak and Mukamel, 2009b) (see sketch in Fig. 27) and describe them using ± representation. For coherent optical states, all field superoperator nonequilibrium Green's functions (SNGF) in the L, R representation are identical E L = E R (the superoperator index makes no difference since all operations commute). In the +, − representation "minus" field indices are not allowed, since E − = 0. The general m wave mixing signals are given by 2 m products of material and corresponding optical field SNGF's of m th order: where t m+1 , t m , . . . , t 1 are the light/matter interaction times.
is the last light-matter interaction. The indexesν j are the conjugates to ν j and defined as follows: the conjugate of + is − and vice versa. Equation (218) implies that is the excitation in the material are caused by fluctuations in the optical field and vice versa. Equation (218) also holds in the L, R representation. Here we have ν m+1 = L and ν j ∈ {L, R}, j = m, . . . , 1.ν j . In the L, R representation the conjugate of "left" is "left" and the conjugate of "right" is "right": L = L,R = R. The material SNGF V  Three wave processes  Heterodyne  Homodyne  Incoherent  Coherent  Technique SFG  DFG  TPIF  TPEF  SFG  PDC  Modes c/c/c c/c/c c/c/q c/q/q c/c/q c/q/q ω3 ω1 + ω2 ω1 − ω2 The SNGF χ

B. Connection to nonlinear fluctuation-dissipation relations
Spontaneous fluctuations and response are uniquely related in the linear regime (Hashitsume et al., 1991) by the fluctuation-dissipation theorem, but not when they are nonlinear. Some nonlinear fluctuation-dissipation relations have been proposed for specific models under limited conditions (Bertini et al., 2001;Bochkov and Kuzovlev, 1981;Lippiello et al., 2008) but there is no universal relation of this type (Kryvohuz and Mukamel, 2012).
A different approach for calculating and analyzing quantum spectroscopy signals is based on the Glauber Sudarshan Prepresentation which expresses the field density matrix as an integral over coherent state density matrices |β β| weighted by a quasi-probability distribution P (β) It has been suggested Kira et al., 2011;Mootz et al., 2014) that, since the response of a material system to a field initially prepared in a coherent state |β is given by the classical response function CRF, the quantum response R QM may be recast as an average of the classical response R |β with respect to this quasi-probability By transforming to the response of a nonclassical state given by P (β), in which the fluctuations along one quadrature are squeezed below the classical limit, this form of data processing can uncover otherwise hidden features in the signal. For instance, it was put to use in (Almand-Hunter et al., 2014) to reveal the "dropleton", as new kind of quasiparticle excitation in semiconductors. While this analysis has proven successful, it misses the multimode nature of the entangled states, which lies at the heart of time-energy entanglement presented here. Based on our analysis of the entangled light state (33), we have shown that optical signals induced by such entangled states cannot be reduced to a simple sum over the signal of each quantum mode: The quantum correlations depicted in Fig. 3 are necessarily missed in a single-mode approach. For instance, the transition amplitude (99) may be written as a sum over the Schmidt modes, where R t shall denote the material response. Hence, the TPIF signal ∼ |T f g | 2 cannot be reduced to the signal of the indivdual Schmidt modes, |T f g | 2 ∼ k,k · · · = k · · · , and the above transformation fails to capture such a signal.
The idea behind the Glauber-Sudarshan quasi-probability in relation to quantum and classical light outlined in Kira et al., 2011;Mootz et al., 2014) makes an equality between classical light and coherent states. This is a potential source of confusion: For instance, the effect of revival of Rabi oscillations demonstrated by Rempe (Rempe et al., 1987) which is observed when a coherent state is treated quantum mechanically shows clearly that a coherent state is generally a quantum state of light. Therefore, Eq. (220) merely represents a transformation between two different basis sets for the quantum state. Of course, in the nonlinear response case, if the operators are normally ordered as in Eq. (74) the result of the coherent state is equivalent to the classical case. In Rempe's example, the nonclassical contributions arise due to commutator terms in Eq. (87) and the higher the order in perturbative expansion -the larger the contribution of the commutator terms. Therefore, in the strong field limit one can observe the revival of Rabi oscillations.
We now use the superoperator notation introduced in section I.A to show more broadly why the quantum response is different from the classical one so that Eq. (221) is violated. In Liouville space, the time-dependent density matrix is given by Eq. (69). We now make use of the algebraic relation of superoperators (Marx et al., 2008): Let us first assume that the electric field operators commute and set E − = 0. We then calculate the expectation value of a system A operator O A : Since the trace of a commutator vanishes, and since there are only V B − operators for system B, all correlation functions of The nonlinear response function is thus additive, and contains no cooperative terms. The time evolution of two coupled quantum systems and the field is generally given by a sum over Feynman paths in their joint phase space. Order by order in the coupling, dynamical observables can be factorized into products of correlation functions defined in the individual spaces of the subsystems. These correlation functions represent both causal response and noncausal spontaneous fluctuations (Cohen and Mukamel, 2003;Roslyak and Mukamel, 2010).
The linear response contains several possible combinations of field superoperators V ± V ± . V − V − represents a commutator, and thus vanishes since its trace is zero. Therefore V + V + (two anticommutators) and V + V − (commutator followed by anticommutator) are the only two quantities that contribute to the linear response. These two quantities are related by the universal fluctuation-dissipation relation (Hashitsume et al., 1991), Here is the response function, whereas denotes spontaneous fluctuations.
The classical response function C +− (ω) thus carries all relevant information about linear radiation matter coupling, including the quantum response. In the nonlinear regime the CRF is a specific causal combination of matter correlation functions given by one "+" and several "-" operators. e.g V + V − V − V − for the third order response. However, the quantum response may also depend on the other combinations. To n th order in the external field the CRF is one member of a larger family of 2 n quantities V + (ω n+1 )V ± (ω n )...V ± (ω 2 )V ± (ω 1 ) representing various combinations of spontaneous fluctuations (represented by V + ) and impulsive excitations (represented by V − ). For example, an "all +" quantity such as V + V + V + V + represents purely spontaneous fluctuations. The CRF does not carry enough information to reproduce all 2 n possible quantities which are accessible by quantum spectroscopy. The reason why the CRF and QRF are not simply related in is the lack of a fluctuationdissipation relation in the nonlinear regime (Bertini et al., 2001;Bochkov and Kuzovlev, 1981;Kryvohuz and Mukamel, 2012;Lippiello et al., 2008).
The field commutator E − is intrinsically related to vacuum modes of the field which may induce coupling between noninteracting parts of the system. One example where such an effect arising from E − is combined with the appearance of collective resonances, which occur for E + , has been recently investigated in the context of harmonic systems (Glenn et al., 2015). The response of classical or quantum harmonic oscillators coupled linearly to a classical field is strictly linear; all nonlinear response functions vanish identically. We have recently shown that quantum modes of the radiation field that mediate interactions between harmonic oscillator resulted in nonlinear susceptibilities. A third-order nonlinear transmission of the optical field yields collective resonances that involve pairs of oscillators, and are missed by the conventional quantum master equation treatment .

C. Heterodyne-detected sum and difference frequency generation with classical light
We first compare two experiments involving three classical modes. The third mode is singled out by the heterodyne detection which measures its time-averaged photon flux (photons 3 . The initial state of the field is given by a direct product of coherent states: |t = −∞ = |β 1 |β 2 |β 3 , where |β α are eigenfunctions of the mode α annihilation operator: a α |β α = β α |β α . Coherent states are the most classical states of quantum light, hence we shall refer to them as classical optical fields (Glauber, 2007a).
Only one type of optical SNGF contribute to the classical third order signalsignal: The conjugate material SNGF's become: where we have assumed that initially the material system is in the ground state, which implies that V (2) LRR (t 3 , t 2 , t 1 ) ≡ 0. We assume the same three-level system used earlier.
The L, R representation plus RWA allows the material SNGFs (228) to be represented by the loop diagrams shown in Fig.28(C), Fig.29(C). The rules for constructing these partially-time-ordered diagrams are summarized in Appendix A (Roslyak et al., 2009a). The signal is given by the causal χ (2) +−− (−ω 3 , ±ω 2 , ±ω 1 ) response function. This result is not new and can be obtained by using various combinations of a classical field from the outset. To calculate equation (228) we have to specify the phase matching conditions (frequencies and wave vectors of the optical modes). This will be done below.

D. Difference frequency generation (DFG)
In DFG the first mode k 1 promotes the system from its ground state |g into state |f . The second mode k 2 induces stimulated emission from |f to an intermediate |e ; and the third mode k 3 stimulates the emission from |e to |g , as sketched in Fig.29 (B). The signal is measured in the phase matching direction: k 3 = k 1 − k 2 (See Fig.28 (A)).
The corresponding loop diagrams corresponding are shown in Fig.28(C). The optical field SNGF yields: The material SNGF's are: The DFG signal in the frequency domain can be written as: Utilizing the rules given in Appendix A and the diagrams shown in Fig.28(C) we obtain: . A sum-over-states (SOS) expansion then yields Equation (233) indicates that the signal induced by classical optical fields is given by the second order classical response function (CRF).

E. Sum Frequency Generation (SFG)
In SFG the first two modes promote the molecule from its ground state |g through intermediate state |e into the state |f . The third mode induces stimulated emission from |f to the ground state |g as sketched in Fig.29 (B). The signal is generated in the direction: k 3 = k 1 + k 2 (See Fig.29 (A)).
The heterodyne-detected SFG signal can be obtained in an analogous manner to the DFG by utilizing the diagrams shown in Fig.29(C): where the CRF is given by In the coming sections equations (235) and (233) will be compared with other techniques involving various combinations of quantum and classical optical fields. These include-homodyne detected SFG, DFG and PDC where one or more optical modes are spontaneously generated and must be treated quantum mechanically.

F. Two-photon-induced fluorescence (TPIF) vs. homodyne-detected SFG
We now turn to techniques involving two classical and one quantum mode where the initial state of the optical field is: |t = −∞ = |β 1 1 |β 2 2 |0 3 . The modes interact with the three level material system (Fig.30(B)). We assume that ω eg = ω ef . This allows to focus on the resonant SNGF's and reduces the number of diagrams.
The two classical modes k 1 and k 2 promote the molecule from its ground state |g into the intermediate state |e and to the final state |f . The system then spontaneously decays back into one of the ground state manifold |g emitting a photon into the third mode k 3 which is initially in the vacuum state (See Fig.30 (B)). The phase-matching condition k 1 − k 1 + k 2 − k 2 + k 3 − k 3 = 0 is automatically satisfied for any k 3 . Therefore the spontaneous photons are emitted into a cone (See Fig.30 (A)).
We shall calculate the photon flux in the k 3 mode. Since the process involves three different modes and six light/matter interactions. The signal (218) must by expanded to fifth order in H int . The fifth order signal S (5) CCQ (CCQ implies two classical and one quantum mode) may be written as sum of 2 5 terms each given by a product of molecule/field six point SNGF's.
Only the corresponding diagrams shown in Fig.30(C). These satisfy the following conditions: 1. The creation operator of the quantum mode a † 3 must be accompanied by the corresponding annihilation operator a 3 .
2. The quantum modes de-excite the molecule, which implies that the annihilation operators must act on the bra and the creation operators act on the ket.
3. The coherent optical fields are tuned off resonance with respect to the ω f g transition. Hence, the signal is not masked by stimulated emission.
To address the collective properties of three wave mixing. We consider a collection of N noninteracting molecules positioned at r i so that V = i V i δ(r − r i ). The optical field can FIG. 30 Three wave process with two classical and one quantum mode: a) phase-matching condition; b) molecular level scheme; c) loop for the incoherent TPIF (C1) and coherent homodyne SFG (C2) interact with different systems at different times. Fig.30(C) shows processes involving three different molecules and six light/matter interactions. These modes can either interact with the same system (incoherent process, Fig.30(C1)) or two different systems (coherent process, Fig.30(C2)) (Marx et al., 2008).

Two-photon-induced fluorescence (TPIF).
This is an incoherent three wave process. Using the identity The relevant material SNGF is: Utilizing equations (236,237) the frequency domain signal can be written as: LR−−−− (−ω 3 ; ω 3 , −ω 2 , ω 2 , −ω 1 , ω 1 ) Note that the above expression is given in the mixed (L/R, +/−) representation. It can be recast into the +, − representation using χ ++−−−− )/2. The second term (two "plus" four "minus" indices) arises since one of the modes is nonclassical. The frequency domain material SNGF can be calculated from the diagram of Fig.30(D1): (239) in molecular states gives: Provided the energy splitting within the ground state manifold is small comparing to the optical transitions the signal can be recast in the form: is the transition amplitude. This is similar to the Kramers-Heisenberg form of ordinary (single-photon) fluorescence (Marx et al., 2008). As in single photon fluorescence, for a correct description of the TPIF the ground state must not be degenerate. Otherwise γ gg = 0 (the degenerate ground state the system has infinite life time) and the signal vanishes.

Homodyne-detected SFG
Here the optical field SNGF is given by equation (236), but instead of factor N it contains the factor: factor: where ∆k ≈ k 1 + k 2 − k 3 . ≈ reflects phase uncertainty given by the reciprocal of the molecular collection length which effectively narrows the optical cone. For large N the coherent part ∝ N (N − 1) dominates over the incoherent ∝ N response. For a small sample size, exact calculation of the optical field part of the SNGF is rather lengthy , but it can be performed in the same fashion as done by Hong and Mandel(Hong and Mandel, 1985) for the probability of photon detection.
To calculate the matter SNGF, we must work in the joint space of two molecules | 1,2 = | 1 | 2 interacting with the same field mode. The matter SNGF of the joint system can be factorized into product of each molecule SNGF's: where we have used the fact that the last interaction must be a "plus".
Since the two molecules are identical the following identity holds: where v g ∆t is the optical path length connecting molecules situated at r i and r j . Using this identity and equation (243) The matter SNGF in the frequency domain can be factorized into the square of an ordinary (causal) response function (one "+" and several "-"): The unique factorization of the coherent matter SNGF's can be also obtained in the L, R representation using the diagrammatic technique as shown in Fig.30(C2). The classical modes have to excite the molecules and the quantum mode has to de-excite them. Hence, the interaction with the first molecule ket (L) are accompanied by the conjugate interactions with the second molecule bra (R). We obtain that the ordinary, causal response function is give by equation (235) as χ (2) +−− (−(ω 1 + ω 2 ); ω 2 , ω 1 ) = χ (2) LLL (−(ω 1 + ω 2 ); ω 2 , ω 1 ). The homodyne-detected SFG signal is finally given by: Both homodyne and heterodyne SFG are given by the same causal response function χ (2) +−− . The main difference is that the latter satisfies perfect phase matching, while for the former this condition is only approximate. For sufficiently large samples the two techniques are identical.
To conclude, we present the total signal for the three wave process involving two classical and one quantum field (CCQ)) which includes both an incoherent and a coherent components: G. Two-photon-emitted fluorescence (TPEF) vs. type-I parametric down conversion (PDC).

TPEF
We now turn to three wave processes involving one classical and two quantum modes. We start with the incoherent response of N identical molecules initially in their ground state. The initial state of the optical field is |t = −∞ = |β 1 1 |0 2 |0 3 . The classical field k 1 pumps the molecule from its ground state |g into the excited state |f . The system then spontaneously emits two photons into modes k 2 , k 3 (See Fig.31 (B)) which are initially in the vacuum state. Such incoherent process which involves one classical and two quantum modes will be denoted two-photon-emitted fluorescence (TPEF). To our knowledge there is neither theoretical nor experimental work concerning this process.
This process is not phase sensitive since the phase matching condition k 1 −k 1 +k 2 −k 2 +k 3 −k 3 = 0 is automatically satisfied for any k 3 . Therefore spontaneously generated modes are emitted into two spatial cones (See Fig.31 (A)). For a single molecule the cones are collinear, as in Type I parametric down conversion (PDC). Due to this similarity we chose beams polarizations as is usually done for PDC of this type: the spontaneously generated photons have the same polarization along the x axis, and orthogonal to that of the classical mode polarized along the y axis.
The time-averaged photon flux in the k 3 mode is our TPEF signal. We again make use of the loop diagrams to identify the relevant SNGF contributing to the signal. Using the initial state of the field the diagrams must satisfy the following conditions: 1. The creation operators of a spontaneously generated modes a † 3 , a † 2 acting on the "ket" must be accompanied by the corresponding annihilation operators a 3 , a 2 acting on the "bra".
2. The first mode ω 1 (k 1 ) is off-resonance with both ω eg and ω f e transitions to avoid stimulated emission contributions.
The matter SNGF assumes the form: LLRR−− (t 6 , t 5 , . . . , t 1 ) = (248) Using equations (247), (248) the incoherent part of the frequency domain signal is given by: The corresponding SNGF can be calculated from the diagram in Fig.31(C1): Expansion in the molecular eigenstates brings the response function into the Kramers-Heisenberg form: Note that unlike TPIF, the TPEF signal depends on the SNGF other than causal response function χ We now turn to the coherent response from a collection of identical molecules which interact with one classical pumping mode and two spontaneously generated quantum modes (See Fig.31 (A,B)). This is known as Type I parametric down conversion (PDC), which is widely used for producing entangled photon pairs. Hereafter we assume perfect phase matching ∆k = k 1 − k 2 − k 3 which is the case for a sufficiently large sample (Gerry and Knight, 2005).
The initial condition for PDC are the same as for TPEF and most PDC experiments are well described by the causal response function χ (2) +−− . Therefore one would expect a connection between TPEF and PDC, similar to that of TPIF and homodyne-detected SFG. However, as we are about to demonstrate, for a complete description of the PDC process the causal second order response function is not enough.
We shall establish the corrections to the second order CRF caused by the quantum origin of the spontaneous modes. To do so we again resort to the CTPL diagrams (See Fig.31 (C2)). For the coherent response the optical field SNGF is given by equation (247) with the factor N replaced by N (N − 1). The material SNGF (248) can be factorized as: Note that this factorization is unique as T V 3 L (t)V 2, † R (t 4 )V 1, † L (t 1 ) = 0 due to the material initially being in its ground state. The coherent PDC signal in the frequency domain can then be written as: LL− (−(ω 1 − ω 2 ); ω 2 , ω 1 ) | 2 Here the generalized response function expanded in the eigenstates has form: This mixed representation can be recast in +, − and L, R representations: Comparing the CRF for heterodyne detected DFG (233) and the SNGF's for Type I PDC (254) we see that the latter is described not only by causal response function χ An interesting effects arise in type-I PDC if the detection process is included explicitly. The details of these effects are discussed in Section IV.I.
In summary of this section we give the signal for the three wave process involving one classical and two quantum fields (CQQ). This contains both an incoherent and a coherent component:

H. Type II PDC; polarization entanglement.
In Type II parametric down-conversion, the two spontaneously-generated signals have orthogonal polarizations. Because of birefringence, the generated photons are emitted along two non-collinear spatial cones known as ordinary and extraordinary beams (See Fig.32 (A)). Polarization-entangled light (Gerry and Knight, 2005) is generated at the intersections of the cones. An x polarization filter and a detector are placed at one of the cones intersections. The detector cannot tell from which beam a photon is obtained. To describe the process we need five optical modes: one classical |1 |y and four quantum modes {|2 |x , |2 |y , |3 |x , |3 |y }.

I. Time-and-frequency resolved type-I PDC
Previously we briefly summarized various wave-mixing signals in a very simple toy model without addressing the detection process in detail. In the following we include the photon counting detection described in Section III.D and demonstrate how does the actual generation process in type-I PDC relate to nonlinear response and how different it is compared to semiclassical theory.
The standard calculation of nonlinear wave mixing assumes that all relevant field frequencies are off resonant with matter. It is then possible to adiabatically eliminate all matter degrees of freedom and describe the process by an effective Hamiltonian for the field that contains a nonlinear cubic coupling of three radiation modes (Mandel and Wolf, 1995). For SFG this reads: (260) and for DFG and PDC: All matter information is embedded into a coefficient that is proportional to χ (2) (Gerry and Knight, 2005) that is defined by the semiclassical theory of radiation-matter coupling. Langevin quantum noise is added to represent vacuum fluctuations caused by other field modes (Avenhaus et al., 2008;Glauber, 2007b;Scully and Zubairy, 1997) and account for photon statistics. The microscopic theory of type I PDC presented below holds if cascade of two photons in a three level system is generated in both on and off resonance. The resonant case is especially important for potential spectroscopic applications (Mukamel, 1995), where unique information about entangled matter (Lettow et al., 2010) can be revealed. Other examples are molecular aggregates and photosynthetic complexes or biological imaging . Second, it properly takes into account the quantum nature of the generated modes through a generalized susceptibility that has a very different behavior near resonance than the semiclassical χ (2) . χ (2) is derived for two classical fields and a single quantum field. While this is true for the reverse process (sum frequency generation) it does not apply for PDC, which couples a single classical and two quantum modes. Third, macroscopic propagation effects are not required for the basic generation of the signal. Fourth, gated detection (Dorfman and Mukamel, 2012a) yields the finite temporal and spectral resolution of the coincident photons limited by a Wigner spectrogram. For either time or frequency resolved measurement of the generated field, the signal can be expressed as a modulus square of the transition amplitude that depends on three field modes. This is not the case for photon counting. Shwartz et al. recently reported PDC in diamond, where 18keV pump field generates two X-ray photons (Shwartz et al., 2012). Calculation in (Dorfman and Mukamel, 2012b) reproduces experimental data in the entire frequency range without adding Langevin noise.
The nature of entangled light can be revealed by photon correlation measurements that are governed by energy, momentum and/or angular momentum conservation. In PDC, a nonlinear medium is pumped by electromagnetic field of frequency ω p and some of the pump photons are converted into pairs of (signal and idler) photons with frequencies ω s and ω i , respectively (see Fig. 33a) satisfying ω p = ω s + ω i .
To get more detailed description of wave-mixing process in PDC we recast an effective semiclassical Hamiltonian in Eq. (261) as follows whereâ †(s) andâ (i) † are creation operators for signal and idler modes, β (p) is the expectation value of the classical pump field, ∆k = k p − k s − k i , j runs over molecules, χ The " + −− indices in Eq. (263) signify two commutators followed by an anti commutator. The bottom line of the semiclassical approach is that PDC is represented by 3point matter-field interaction via the second order susceptibility χ (2) +−− that couples the signal, idler and pump modes. However, it has been realized, that other field modes are needed to yield the correct photon statistics. Electromagnetic field fluctuations are then added as quantum noise (Langevin forces) (Scully and Zubairy, 1997).
Below we summarize a microscopic calculation of the photon coincidence counting (PCC) rate in type I PDC (Roslyak et al., 2009a). We show that PDC is governed by a quantity that resembles but is different from (263). In contrast with the semiclassical approach, PDC emerges as a 6-mode twomolecule rather than 3-mode matter-field interaction process and is represented by convolution of two quantum susceptibilities χ (2) LL− (ω s , ω i , ω p ) and χ (2) * LL− (ω s , ω i , ω p ) that represent a pair of molecules in the sample interacting with many vacuum modes of the signal (s, s ) and the idler (i, i ). Field fluctuations are included self consistently at the microscopic level. Furthermore the relevant nonlinear susceptibility is different from the semiclassical one χ (2) +−− and is given by Eq. (264) has a single commutator and is symmetric to a permutation of ω i and ω s = ω p − ω i , as it should. Eq. (263) has two commutators and lacks this symmetry.

The bare PCC rate
The measured PCC rate of signal and idler photons starts with the definition (175). The basis for the bare signal in Eq. (174) is given by the time-ordered product of Green's functions of superoperators in the interaction picture (see (Dorfman and Mukamel, 2012b)).
which represents four spectral modes arriving at the detectors, where modes s, s , i, i are defined by their frequencies ω s , ω i , In type I PDC the sample is composed of N identical molecules initially in their ground state. They interact with one classical pump mode and emit two spontaneously generated quantum modes with the same polarization into two collinear cones. The initial state of the optical field is given by |0 s |0 i |β p . A classical pump field promotes the molecule from its ground state |g to the doubly excited state |f (see Fig 33b).
Due to the quantum nature of the signal and the idler fields, the interaction of each of these fields with matter must be at least second order to yield a non vanishing signal. The leading contribution to Eq. (265) comes from the four diagrams shown in Fig. 33c (for rules see Appendix A (Roslyak et al., 2009a)). The coherent part of the signal represented by interaction of two spontaneously generated quantum and one classical modes is proportional to the number of pairs of sites in the sample ∼ N (N − 1), which dominates the other, incoherent, ∼ N response for large N . Details of the calculation of the correlation function (265) are presented in Ref. (Dorfman and Mukamel, 2012b). We obtain for the "bare" frequency domain PCC rate R where E (p) (ω) ≡ E (p) (ω)β (p) is a classical field amplitude, E (p) (ω) is the pump pulse envelope and D(ω) = ωD(ω) whereD(ω) = V ω 2 /π 2 c 3 is the density of radiation modes. For our level scheme (Fig. 33b) the nonlinear susceptibility χ Eq. (266) represents a 6-mode (ω p , ω i , ω s , ω p , ω i , ω s ) field-matter correlation function factorized into two generalized susceptibilities each representing the interaction of two quantum and one classical mode with a different molecule. Because of two constraints ω p = ω s + ω i , ω p = ω s + ω i that originate from time translation invariance on each of the two molecules that generate the nonlinear response, Eq. (266) only depends on four field modes. Each molecule creates a coherence in the field between states with zero and one photon. By combining the susceptibilities from a pair of molecules we obtain a photon occupation number that can be detected. Thus, the detection process must be described in the joint space of the two molecules and involves the interference of four quantum pathways (two with bra and two with ket) with different time orderings. Note that this pathway information is not explicit in the Langevin approach.
For comparison, if all three fields (signal, idler and pump) are classical, the number of material-field interactions is reduced to three -one for each field. Then the leading contribution to the field correlation function yields the semiclassical  nonlinear susceptibility χ where G(ω) = 1/(ω − H 0− / + iγ) is the retarded Liouville Green's function, and γ is lifetime broadening. χ (2) LL− possesses a permutation symmetry with respect to s ↔ i (both have L index). In contrast the semiclassical calculation via χ (2) ++− is non symmetric with respect to s ↔ i ( one + and one − indexes), which results in PCC rate that depends upon whether the signal or idler detector clicks first (Shwartz et al., 2012).

Simulations of typical PDC signals
Fig. 33d-k compares both susceptibilities calculated for a typical KTP crystal (PPKTP) represented by a degenerate three-level system with parameters taken from Ref. (D'auria et al., 2008;Wu et al., 1986). Fig. 33f,h,j show that far from resonances (ω i = ω eg , ω p = ω gf ) the semiclassical and quantum susceptibilities coincide and depend weakly on the frequencies ω p and ω i . This is the regime covered by the semiclassical theory, where the susceptibility is assumed to be a constant. Similar agreement between classical and quantum susceptibilities can be observed if the pump is resonant with two-photon transition ω p ω gf but the idler is off resonance Fig. 33g,h,k. However, close to resonance - Fig.  33d,e the two susceptibilities are very different. The semiclassical susceptibility χ (2) +−− vanishes at resonance, where the quantum susceptibility χ (2) LL− reaches its maximum. To put the present ideas into more practical perspective we show in Fig. 34 the PCC rate for a monochromatic pump ω p and mean signal detector frequencyω s . The quantum theory yields one strong resonant peak atω s = ω p − ω eg and two weak peaks at ω p = ω gf andω s = ω eg if the idler detector is resonant with the intermediate state |e :ω i = ω eg - Fig. 34a. However, if we tune the idler detector to a different frequency, for instance ω eg −ω i = 2GHz there is an additional peak at ω s = ω p −ω i - Fig. 34b. Similarly, when the pump consists of two monochromatic beams ω p = ω p (panels c,d) the number of peaks are doubled compare to single monochromatic pump. Clearly, one can reproduce the exact same peaks for ω p as for ω p .

Spectral diffusion
Spectral diffusion (SD) which results from the stochastic modulation of frequencies can manifest itself either as discrete random jumps of the emission frequency (Santori et al., 2010;Siyushev et al., 2009;Walden-Newman et al., 2012) or as a broadening of a hole burnt in the spectrum by a narrowband pulse (Wagie and Geissinger, 2012;Xie and Trautman, 1998). We focus on the latter case assuming that the electronic states of molecule α = a, b are coupled to a harmonic bath described by the HamiltonianĤ α B = k ω k (â †α kâ α k + 1/2) (see Fig. 36b). The bath perturbs the energy of state ν. This is represented by the Hamiltonian whereq ν is a collective bath coordinatê d mn,k represents bath-induced fluctuations of the transition energies (m = n) and the intermolecular coupling (m = n).
We define the line-shape function where the bath spectral density is given by β = k B T with k B being the Boltzmann constant and T is the ambient temperature. We shall use the overdamped Brownian oscillator model for the spectral density, assuming a single nuclear coordinate (ν α = ν α ) where Λ α is the fluctuation relaxation rate and λ α is the system-bath coupling strength. The corresponding lineshape function then depends on two parameters: the reorganization energy λ α denoting the strength of the coupling to the bath and the fluctuation relaxation rate Λ α = 2λ α k B T / . In the high-temperature limit k B T Λ α we have For a given magnitude of fluctuations ∆ α , α = a, b the FWHM of the absorption linewidth (Mukamel, 1995) The absorption and emission lineshape functions for a pair of molecules obtained in the slow nuclear dynamics limit: Λ α ∆ α are given by (Mukamel, 1995) σ where 2λ α is the Stokes shift and ∆ α = 2λ α k B T / is a linewidth parameter. Together with the relaxation rate Λ α (see Eq. (274)) these parameters completely describe the SD model and govern the evolution of the emission linewidth between the initial time given by Eq. (276) (see Fig. 35a) and long time given by Eq. (277) (see Fig. 35b).
The corresponding four-point matter correlation function may be obtained by the second order cumulant expansion (Mukamel, 1995) where ω α ≡ ω eα − ω gα is the absorption frequency, Φ α (t 1 , t 2 , t 3 , t 4 ) is the four-point lineshape function We focus on the SD in the "hole burning" limit (HBL) which holds under two conditions: First, the dephasing is much faster than the fluctuation timescale, i.e. t k Λ −1 α , k = 1, 2, 3, 4. Second, if excitation pulse duration σ −1 p and the inverse spectral (σ j ω ) −1 , and temporal (σ j T ) −1 , j = r, s gate bandwidths are much shorter than the fluctuation time scales, one may neglect the dynamics during the delay between population evolution and its detection. This parameter regime is relevant to crystals which store information in the form of reversible notches that are created in their optical absorption spectra at specific frequencies. Long storage times (Longdell et al., 2005), high efficiencies (Hedges et al., 2010), and many photon qubits in each crystal (Shahriar et al., 2002) can be achieved in this limit. The HBL limit is natural for long-term quantum memories where entanglement is achieved with telecom photons, proved the possibility of quantum internet (Clausen et al., 2011;Saglamyurek et al., 2011).
Time-and-frequency resolved fluorescence is the simplest way to observe SD. The molecular transition frequency is coupled linearly to an overdamped Brownian oscillator that represents the bath (see Fig. 36b). This fluorescence signal given by Eq. (284) is depicted as a series of the snapshot spectra at different times for molecule a in Fig. 35c. It shows a time dependent frequency redshiftω a (t) and time dependent spectral broadening given byσ a0 (t). Initiallyω a (0) = ω 0 a + λ a whereas at long timesω a (∞) = ω 0 a − λ a , where 2λ α is the Stokes shift. Same for molecule b is shown in Fig. 35d. Because of the different reorganization energies λ a , λ b and relaxation rates Λ a , Λ b the Stokes shift dynamics and dispersion are different. Even when the absorption frequencies are the same ω a = ω b , the fluorescence can show a different the profile due to SD. This affects the distinguishability of the emitted photons as will be demonstrated below.

J. Generation and entanglement control of photons produced by two independent molecules by time-and-frequency gated photon coincidence counting (PCC).
As discussed in section II, entangled photons can be produced by a large variety of χ (2) -processes. Alternatively, photons can be entangled by simultaneous excitation of two remote molecules and mixing the spontaneously emitted photons on the beam splitter. This entangled photon generation scheme is suitable e.g. for spectroscopic studies of single molecules by single photons discussed below. The manipulation of single photon interference by appropriate time-andfrequency gating discussed below is presented in details in Ref. (Dorfman and Mukamel, 2014a). The generated entangled photon pair can be used in logical operations based on optical measurements that utilize interference between indistinguishable photons.
We examine photon interference in the setup shown in Fig.  36a. A pair of photons is generated by two remote twolevel molecules a and b with ground g α and excited state e α , α = a, b. These photons then enter a beam splitter and are subsequently registered by time-and-frequency gated detectors s and r. There are three possible outcomes: two photons registered in detector s, two photons registered in r or coincidence where one photon is detected in each. The ratios of these outcomes reflects the photon Bose statistics and depends on their degree of indistinguishability. If the two photons incident on the beam splitter are indistinguishable the photon coincidence counting signal (PCC) vanishes. This causes the Hong-Ou-Mandel (HOM) dip when varying the position of the beam splitter which causes delay T between the two photons. The normalized PCC rate varies between 1 for completely distinguishable photons (large T ) and 0 when they are totally indistinguishable (T = 0). For classical fields and 50:50 beam splitter the PCC rate may not be lower than 1/2. We denote the photons as indistinguishable (distinguishable) if the PCC rate is smaller (larger) than 1/2. PCC is typically measured by time-resolved detection (Gerry and Knight, 2005;Glauber, 2007b). Originally performed with entangled photons generated by parametric down conversion (PDC) (Hong et al., 1987) the shape of the dip vs delay is usually related to the two-photon state envelope which is governed by an effective PDC Hamiltonian (Dorfman and Mukamel, 2012b). Bath effects can become important for remote emitters and have been introduced phenomenologically (Lettow et al., 2010) below. We present a microscopic description of PCC with bath fluctuations by formulating the signal in the joint field-matter space measured by simultaneous time-and-frequency resolved detection. This complex measurement can be achieved using high-speed photodiode which converts a fast optical signal into a fast electric current, fast oscilloscopes to observe the waveform, wide bandwidth spectrum analyzers and other elements. Short pulse characterization using time-frequency map such as frequency-resolved optical gating (FROG) (Trebino, 2000), spectral phase inter-ferometry for direct-field reconstruction (SPIDER) (Dorrer et al., 1999) are well established tools for ultrafast metrology (Walmsley and Dorrer, 2009;Wollenhaupt et al., 2007). Extending these techniques to a single photon time and frequency resolved detection is challenging and may be achieved if combined with on-chip tunable detectors (Gustavsson et al., 2007) or upconversion processes (Gu et al., 2010;Ma et al., 2011).
1. PCC of single photons generated by two remote emitters.
The time-and-frequency gated PCC signal is described by the two pairs of loop diagrams shown in Fig. 36c. Each loop represents one molecule (a or b) which undergoes four fieldmatter interactions and each detector interacts twice with the field. Fig. 36c shows that after interacting with the pump (with its ket) at time t 2 molecule a evolves in the coherence ρ eaga during time interval t 2 . The second interaction of the pump with the bra then brings the molecule into a population state ρ eaea which evolves during interval t 1 until the first interaction with spontaneous emission mode occurs with ket. The molecule then evolves into the coherence ρ gaea during t 1 until the second bra-interaction of spontaneous mode. During population and coherence periods, the characteristic timescale of the dynamics is governed by population relaxation and dephasing, respectively.
The relevant single particle information for molecule α is given by the four point dipole correlation function where V and V † are the lowering and raising dipole transition operators, respectively. Diagrams i in Fig. 36c represent noninterfering terms given by a product of two independent fluorescence contributions of the individual molecules. Diagrams ii represent an interference in the joint space of the two molecules and involves the interference of eight quantum pathways (four with the bra and four with the ket) with different time orderings. Each molecule creates a coherence in the field between states with zero and one photon |0 1| and |1 0|. By combining the contributions from a pair of molecules we obtain a photon population |1 1| that can be detected (Dorfman and Mukamel, 2012a,b). For a pair of identical molecules, the beam splitter destroys the pathway information making the molecules indistinguishable and giving rise to quantum interference.
The PCC signal (175) (see Ref. (Dorfman and Mukamel, 2014a)) has been calculated for the output fields E 3 and E 4 (i = 3, s = 4) of the beam splitter (see Fig. 36a). The fields in the output 3, 4 and input 1, 2 ports of the 50:50 beamsplitter are related by where ±cT is small difference of path length in the two arms. Taking into account that Eq. (175) should be modified to include the beam splitter position and absorb it into the gating spectrograms, the time-and-frequency resolved PCC in this case is given by Here Γ j = {t j , ω j } represents the set of parameters of the matter plus field incident on the detector j = r, s. Eq. (280) is given by the spectral and temporal overlap of the Wigner spectrograms of detectors W Eqs. (281) -(282) contain all relevant field matter interactions. The modified expressions for the gating spectrograms that include delay T are given in (Dorfman and Mukamel, 2014a).

Time-and-frequency gated PCC
Under SD and HBL conditions the PCC signal (175) is given by where expressions in the last line represent permutation of the molecules a and b, Γ j = {t j ,ω j } represents a set of gating parameters for the detector j = r, s. C α (Γ = {t, ω}) is the time-and-frequency resolved fluorescence of molecule α = a, b corresponding to diagram i in Fig. 36c: ω 0 α = ω α − λ α is the mean absorption and fluorescence frequency. I j α (Γ 1 , t 2 , τ ) and I j α (t 1 , Γ 2 , τ ) with t 1 < t 2 is the interference contribution α = a, b, j = r, s corresponding to diagram ii in Fig. 36c and all the remaining parameters are listed in Ref. (Dorfman and Mukamel, 2014a). The contribution of Eq. (284) is represented by an amplitude square coming from each molecule in the presence of fluctuations. The interference term (285) generally cannot be recast as a product of two amplitudes . Eqs. (283) -(285) are simulated below using the typical parameters of the two photon interference experiments (Ates et al., 2012;Bernien et al., 2012;Coolen et al., 2008;Lettow et al., 2010;Patel et al., 2010;Sanaka et al., 2012;Santori et al., 2002;Sipahigil et al., 2012;Trebbia et al., 2010;Wolters et al., 2013).
3. Signatures of gating and spectral diffusion in the HOM dip.
Photon indistinguishability depends on the molecular transition frequencies. Fig. 35e shows that for fixed time and frequency gate bandwidths σ j ω and σ j T , j = r, s the photons are distinguishabe as long as the transition energy offset ω 0 b −ω 0 a is larger than the gate bandwidth and are indistinguishable otherwise. The SD timescale is a key parameter affecting the degree of indistinguishability. Using Eq. (275) we fixed the absorption linewidth Γ α and varied Λ α and ∆ α . The PCC signal (283) depicted in Fig. 35f shows that if the molecules have nearly degenerate transition energy offset for slower fluctuations the photons are indistinguishable. Increasing the SD rate of one of the molecules increases the photon distinguishability when both time and frequency gates are broader than the difference in transition frequencies.
We further illustrate the effect of frequency and time gating on spectral diffusion. Fig. 35g shows that if two molecules have different SD timescales and the frequency gate bandwidth is narrow the photons are rendered distinguishable and HOM dip is 0.6. By increasing the σ ω the photons become indistinguishable and HOM dip is 0.48 for σ ω = 120 MHz and 0.35 for σ ω = 200 MHz. In all three cases we kept the time gate fixed. Alternatively we change the time gate bandwidth while keeping the frequency gate fixed. Fig. 35h shows that initially indistinguishable photons at σ T = 80 MHz with HOM dip 0.675 become indistinguishable at σ T = 110 MHz with HOM dip 0.45 and at σ T = 150 MHz with HOM dip 0.275. Thus, if the presence of the bath erodes the HOM dip the manipulation of the detection gating allows to preserve the quantum interference

V. SUMMARY AND OUTLOOK
The term quantum spectroscopy broadly refers to spectroscopy techniques that make use of the quantum nature of light. Photon counting studies obviously belong to this category. Studies that detect the signal field such as heterodyne detection or fluorescence are obtained by expanding the signals in powers of the field operators. These depend on multipoint correlation functions of the incoming fields. Classical spectroscopy is then recovered when these can be factorized into products of field amplitudes. Otherwise the technique is quantum. Spectroscopy is classical if all fields are in a coherent state and the observable is given by normally ordered product.
Another important aspect of quantum spectroscopy, which we had only touched briefly in section IV.B, concerns the nonclassical fluctuations of quantum light and their exploitation as spectroscopic tools (Benatti et al., 2010;Giovannetto et al., 2011). These may also lead to novel features in nonlinear optical signals (López Carreño et al., 2015).
We now classify and briefly survey the main features of quantum spectroscopy. First, the unusual time/frequency windows for homodyne, heterodyne and fluorescence detection arising due to the quantum nature of the light generation resulting in the enhanced resolution of the signals not accessible by classical light. Second, photon counting and interferometric detection schemes constitute a class of multidimensional signals that are based on detection and manipulation of a single photons and are parametrized by the emitted photons rather by the incoming fields. Third, the quantum nature of light manifest in collective effects in many-body systems by projecting entanglement of the field to matter. This allows to e.g. prepare and control of higher excited states in molecular aggregates, access dark multi particle states, etc. Fourth, due to the lack of fluctuation dissipation relations, quantum light can manifest new combinations of field and corresponding matter correlation functions not governed by semiclassical response functions such as in parametric down conversion, sum-or difference-frequency generation, two-photon-induced fluorescence, etc. Finally, elaborate pulse shaping techniques that have been recently scaled down to single photon level provide an additional tool for multidimensional measurements using delay scanning protocols not available for classical laser experiments.
The reason for the potential advantage of quantum spectroscopy was traced back to the strong time-frequency correlations inherent to quantum light, and the back-action of the interaction events onto the quantum field's state. The combination of the two effects leads to the excitation of distinct wavepackets, which can be designed to enhance or suppress selected features of the resulting optical signals.
We have described the nonlinear signals in terms of convolutions of multi-time correlation functions of the field and the matter. This approach naturally connects to the established framework of quantum optics, where field correlation functions are analyzed (Glauber, 1963), with nonlinear laser spectroscopy, which investigates the information content of matter correlation functions. As such, it provides a flexible platform to explore quantum light interaction with complex systems well beyond spectroscopic applications. This could include coherent control with quantum light (Schlawin and Buchleitner, 2015;Wu et al., 2015), or the manipulation of ultracold atoms with light (Mekhov and Ritsch, 2012). Entangled quantum states with higher photon numbers (Shalm et al., 2013) promise access to the χ (5) -susceptibility and its additional information content. Combination of quantum light with strong coupling to intense fields in optical cavities (Herrera et al., 2014;Hutchison et al., 2012;Schwartz et al., 2013)  blue lines indicate the evolution of the matter density matrix. An arrow pointing towards (away from) it corresponds to a matter excitation V † (τ ) (de-excitation V (τ )), which is accompanied by a photon annihilation E(τ ) [creation E † (τ )]. In the Liouville space formulation, interactions on the left side represent left-superoperators (1), and interaction on the right right-superoperators (2).

VI. ACKNOWLEDGEMENTS
We gratefully acknowledge the support of the National Science Foundation through Grant No. CHE-1361516, the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy. KD was supported by the DOE grant. We also thank the support and stimulating environment at FRIAS (Freiburg Institute for Advanced Studies). FS would like to thank the German National Academic Foundation for support.

Appendix A: Diagram construction
In leading (fourth-) order perturbation theory, the population of a two-excitation state f is given by where P f (t) = |f (t) f (t)| is the projector onto the final state, the interaction Hamiltonian H int as given in Eq. (7), and Γ denotes the set of control parameters in the light field. Eq. (A1) contains 4 4 = 256 terms, of which only very few contribute to the signal. These contributions can be conveniently found using a diagrammatic representation. Its basic building blocks are shown in Fig. 37: The evolution of both the bra and the ket side of the density matrix in Eq. (10) [ ψ| and |ψ in Eq. (17)] are represented by vertical blue lines, and sample excitations (de-excitations) are represented by arrows pointing towards (away from) the density matrix.

Loop diagrams
The following rules are used to construct the field and matter correlation function from the diagrams (Marx et al., 2008): 1. Time runs along the loop clockwise from bottom left to bottom right. 2. The left branch of the loop represents the "ket", the right represents the "bra".

= tr
3. Each interaction with a field mode is represented by an arrow line on either the right (R-operators) or the left (L-operators).
4. The field is marked by dressing the lines with arrows, where an arrow pointing to the right (left) represents the field annihilation (creation) operator E α (t) (E † α (t)).

5.
Within the RWA, each interaction with the field annihilates the photon E α (t) and is accompanied by applying the operator V † α (t), which leads to excitation of the state represented by ket and deexcitating of the state represented by the bra, respectively. Arrows pointing "inwards" (i.e. pointing to the right on the ket and to the left on the bra) consequently cause absorption of a photon by exciting the system, whereas arrows pointing "outwards" (i.e. pointing to the left on the bra and to the right on the ket) represent deexcitating the system by photon emission.
6. The observation time t, is fixed and is always the last.
As a convention, it is chosen to occur from the left. This can always be achieved by a reflection of all interactions through the center line between the ket and the bra, which corresponds to taking the complex conjugate of the original correlation function.
7. The loop translates into an alternating product of interactions (arrows) and periods of free evolutions (vertical solid lines) along the loop.
8. Since the loop time goes clockwise along the loop, periods of free evolution on the left branch amount to propagating forward in real time with the propagator give by the retarded Green's function G. Whereas evolution on the right branch corresponds to backward propagation (advanced Green's function G † ).
9. The frequency arguments of the various propagators are cumulative, i.e. they are given by the sum of all "earlier" interactions along the loop. Additionally, the ground state frequency is added to all arguments of the propagators.
10. The Fourier transform of the time-domain propagators adds an additional factor of i(−i) for each retarded (advanced) propagator.
11. The overall sign of the SNGF is given by (−1) N R , where N R stands for the number of R superoperators.

Ladder diagrams
The same rules may be applied to evaluate ladder diagrams, with the exception of rule 1, which reads 1. Time runs from bottom to top. where we defined the Hermite functions H n . Here, we have a = 1/(2σ 2 p ) + γT 2 1 , b = 1/(2σ 2 p ) + γT 1 T 2 , and c = 1/(2σ 2 p ) + γT 2 2 . The Hermite functions are then given by with the Hermite polynomials h n , and we further defined the parameters and r n = α 1 + µ 2 4acσ 2 p µ n , )µ e f G f e (ω a + ω b − ω a )µ e f G ee ;e e (ω a − ω a )µ ge G eg (ω a )µ ge , )µ e f G ee ;e e (ω a − ω a )µ ge G eg (ω a )µ ge . (C11) for photon absorption at time t can be calculated in first-order perturbation theory, which yields (Glauber, 2007b) ψ f |E(t)|ψ i · a|d|b , where d is the dipole moment of the atom and E(t) = E † (t) + E(t) is the electric field operator (we omit the spatial dependence). Clearly, only the annihilation part of the electric field contributes to the photon absorption process. The transition probability to find the field in state ψ f at time t is given by the modulus square of the transition amplitude Since the initial state of the field ψ i is rarely known with certainty, we must trace over all possible initial states as determined by a density operator ρ. Thus, the output of the idealized detector is more generally given by tr ρE † (t)E(t) .
Simultaneous time-and-frequency resolved measurement must involve a frequency (spectral) gate combined with time gate -a shutter that opens up for very short interval of time. The combined detector with input located at r G is represented by a time gate F t centered att followed by a frequency gate F f centered atω (Stolz, 1994). First, the time gate transforms the electric field E(r G , t) = sÊ s (r G , t) withÊ s (r G , t) = E(r G , ω s )e −iωst as follows: Then, the frequency gate is applied E (tf ) (t,ω; r G , ω) = F f (ω,ω)E (t) (t; r G , ω) to obtain the time-and-frequencygated field. We assume that the time gate is applied first. Therefore, the combined detected field at the position r D can be written as where E(t) is the electric field operator (8) in the Heisenberg picture. Similarly, one can apply the frequency gate first and obtain frequency-and-time-gated field E (f t) .
The following discussion will be based on Eq. (E4). Eq. (E5) can be handled similarly. For gaussian gates whereσ 2 ω = σ 2 T + σ 2 ω ,σ 2 Note that σ T and σ ω can be controlled independently, but the actual time and frequency resolution is controlled byσ T andσ ω , respectively, which always satisfy Fourier uncertaintỹ σ ω /σ T > 1. For lorentzian gates the detector time-domain and Wigner spectrograms are given by The gated signal is given by (E13) where the angular brackets denote ... ≡ Tr[ρ(t)...]. The density operator ρ(t) is defined in the joint field-matter space of the entire system. Note, that Eq. (E13) represents the observable signal, which is always positive since it can be recast as a modulus square of an amplitude (Eq. (E2)).For clarity we hereafter omit the position dependence in the fields assuming that propagation between r G and r D is included in the spectral gate function.

The bare signal
The bare signal assumes infinite spectral and temporal resolution. It is unphysical but carries all necessary information for calculating photon counting measurement. It is given by the closed path time-loop diagram shown in Fig. 21 . We assume an arbitrary field-matter evolution starting from the matter ground state g that promotes the system up to some excited state. The system then emits a photon with frequency ω s that leaves the matter in the state e. This photon is later absorbed by the detector.
In the absence of dissipation (unitary evolution) the matter correlation function can be further factorized into a product of two amplitudes that correspond to unitary evolution of braand ket-. This transition amplitude can be recast in Hilbert space without using superoperators and is given by This gives for the bare Wigner spectrogram n(t , ω ) = e ∞ 0 dτ e −iω τ ×T eg (t − τ /2)T * eg (t + τ /2). (E15)

Spectrogram-overlap representation for detected signal
In the standard theory (Stolz, 1994), the detected signal is given by a convolution of the spectrograms of the detector and bare signal. The detector spectrogram is an ordinary function of the gating parameters whereas the bare signal is related to the field prior to detection. We now show that when the process is described in the joint matter plus field space the signal can be brought to the same form, except that now the bare signal is given by a correlation function of matter that further includes a sum over the detected modes. We denote this the spectrogram-overlap (SO) representation of the signal. Alternatively one can introduce a spectrogram-superoperatoroverlap (SSO) representation where field modes that interact with detector are included in the detector spectrogram, which becomes a superoperator in the field space. Details of this representation are presented in Ref. (Dorfman and Mukamel, 2012a). Below we present the signals in the time domain, which can be directly read of the diagram (Fig. 21b). We then recast them using Wigner spectrograms, which depict simultaneously temporal and spectral profiles of the signal. We now define the detector Wigner spectrogram if the spectral gate applied first, using Eq. (E5). The detector spectrogram is alternatively given by Combining Eqs. (E4) -(E16) we obtain that the gated signal is given by the temporal overlap of the bare signal and detector Wigner spectrogram W D (t,ω; t , ω )n(t , ω ).
This is the conventional form (Stolz, 1994) introduced originally for the field space alone. Eq. (E15) contains explicitly the multiple pairs of radiation modes s and s that can be revealed only in the joint field plus matter space. Eventually this takes into account all the field matter interactions that lead to the emission of the detected field modes. All parameters of F f and F t can be freely varied. The spectrogram will always satisfy the Fourier uncertainty ∆t∆ω > 1.
Together with the gated spectrogram (E16) the bare signal (E15) represents the time and frequency resolved gated signal. Note, that in the presence of a bath, the signal (E15) is no longer given by a product of two amplitudes.T eg (t) is then an operator in the space of the bath degrees of freedom. Therefore, one has to replace product of amplitudes in Eq. (E15) by T eg (t − τ /2)T * eg (t + τ /2) , where ... corresponds to averaging over the bath degrees of freedom. The convolution of two operatorsT eg reveals the multiple pathways between these initial and final states of matter that allows to observe them through the simultaneous time and frequency resolution.
We now consider two limiting cases. In the absence of a frequency gate, then F f (ω,ω) = 1 we get W D (ω,t; t, τ ) = δ(τ )F * t (t + τ /2,t)F t (t − τ /2,t). For the narrow time gate |F t (t,t)| 2 = δ(t −t) we then obtain the time resolved measurement nt = e |T eg (t)| 2 . (E22) In the opposite limit, where there is no time gate, i.e. F t (t,t) = 1, and the frequency gate is very narrow, such that F f (t,ω) = √ γ π e −iωt−γt θ(t) at γ → 0, then W D (ω,t; t, τ ) = e −iωτ . In this case we obtain the frequency resolved measurement nω = e |T eg (ω)| 2 , where T eg (ω) = ∞ −∞ dte iωt T eg (t). Eqs. (E22) and (E23) indicate that if the measurement is either purely time or purely frequency resolved, the signal can be expressed in terms of the modulus square of a transition amplitude. Interference can then occur only within T eg in Hilbert space but not between the two amplitudes. Simultaneous time and frequency gating also involves interference between the two amplitudes; the pathway is in the joint ket plus bra density matrix space.
In the presence of a bath, the signal can be written as a correlation function in the space of bath coordinates T * eg (t)T eg (t) for Eq. (E22) and T * eg (ω)T eg (ω) for Eq. (E23).

Multiple detections
The present formalism is modular and may be easily extended to any number of detection events. To that end it is more convenient to use the time domain, rather than Wigner representation. For coincidence counting of two photons measured by first detector with parametersω i ,t i followed by second detector characterized byω s ,t s the time-and-frequency resolved measurement in SO representation is given by where the detector spectrogram for mode ν = i, s reads D(t ν ,ω ν ; t ν , τ ν ) = ∞ −∞ dω ν 2π e −iων τν × |F f (ω ν ,ω ν )| 2 F * t (t ν + τ ν /2,t ν )F t (t ν − τ ν /2,t ν ).
The bare signal is given by the loop diagram in Fig. 21c which reads B(t s , τ s ; t i , τ i ) = − e T eg (t s − τ s /2, t i − τ i /2)T * eg (t s + τ s /2, t i + τ i /2).

(E26)
The transition amplitude for the ket reads and for the bra