Two-photon spectra of quantum emitters

We apply our recently developed theory of frequency-filtered and time-resolved N-photon correlations to study the two-photon spectra of a variety of systems of increasing complexity: single mode emitters with two limiting statistics (one harmonic oscillator or a two-level system) and the various combinations that arise from their coupling. We consider both the linear and nonlinear regimes under incoherent excitation. We find that even the simplest systems display a rich dynamics of emission, not accessible by simple single photon spectroscopy. In the strong coupling regime, novel two-photon emission processes involving virtual states are revealed. Furthermore, two general results are unraveled by two-photon correlations with narrow linewidth detectors: i) filtering induced bunching and ii) breakdown of the semi-classical theory. We show how to overcome this shortcoming in a fully-quantized picture.


Introduction
At the heart of the modern theory of light lies the concept of the photon, the quantum of the electromagnetic field [2]. Regardless of the wave behaviour of light in a variety of contexts, ultimately, every measurement is accounted for by clicks on a detector. When the number of quanta is so large that the granularity of the field disappears due to coarse-graining, the limit of classical electrodynamics is recovered [3]. It is one of the wonders of quantum mechanics that the light field still exhibits wave-like properties when tracked at the single-particle level. Most popularly, the interference patterns in a double slit experiment is reconstructed by the accumulation of a large number of single clicks, each observed in isolation [4].
Usually, even in the regime of a small number of quanta in the field, the experimental observables are obtained by averaging over a large enough repetition of the experiment. For instance, a luminescence spectrum is obtained by integrating the signal over long times, such as to obtain a continuous distribution that describes the probability of emission of each photon at a given energy. An example for such a one-photon spectrum (1PS) S(ω) is given in Fig. 1(a) for the simplest conceivable case, that of a single mode emitter. The lineshape is simply a Lorentzian with a linewidth that is given by the inverse lifetime, corresponding, in the time domain, to an exponential decay of the excitation: the emission is memoryless and has the same probability to occur at any given time, as befits a quantum particle. This results in a spread of energy of the emitted photon. Arguably, however, this lineshape is also typical of a classical system, in the form of a damped resonance.
To access the dynamics of photons at a quantum level, with no counterpart from a classical theory involving a continuous field, it is necessary to invoke N -photon correlations. An average is still taken over a large number of single-detection events, but the latter now involves N correlated clicks in each shot. In this text, we will focus on the simplest case where N = 2. There is an ample literature on two-photon correlations [5]. Typically, photons are detected as a function of time (and/or other commuting variables like space) in photon counting configurations such as the celebrated Hanburry Brown-Twiss setup [6]. The specificity of our approach is that we study correlations in the frequency domain, i.e., we compute correlations between a photon detected at an energy ω 1 and another at energy ω 2 . The conjugate time information is still present through the linewidth of the frequency windows. Real time can also be included, computing correlations between a photon of frequency ω 1 detected at time t 1 and another of frequency ω 2 detected at time t 2 , although we will limit here to simultaneous detection for brevity. This is a generalisation of the conventional (color-blind) theory. The two-photon spectroscopy of a quantum emitter is obtained when such correlations are computed (theoretically) or measured (experimentally) for all the energies spanning the range of emission of the system. Experimentally, the implementation is straightforward, as it essentially consists in inserting filters between the system and the conventional photon-counting setup. The technique is increasingly popular in the laboratory [7,8,9,10,11] although its interpretation relies chiefly on physical grounds or numerical simulations. Theoretically, almost all the effort has been dedicated to resonance fluorescence [12,13,14,15,16], owing to the difficulty of actual computation for other systems, prior to our input [1] that yields exact results for arbitrary systems.
In this text, we present a comprehensive analysis of a variety of fundamental quantum optical systems and show the extremely rich dynamics and the novel insights that such frequency-resolved two-photon correlations unravel. In Section 2, we give a brief overview of our theory to compute two-photon spectra (2PS), denoted g (2) Γ (ω 1 , ω 2 ). The notation refers to the Glauber precept that places photon correlations as the pillar of quantum optics [5]. The mode which photons are collected from will be appended in square brackets for clarity but only when necessary to avoid confusion. The full theory for N -photon correlations and/or arbitrary time delays is laid down in Ref. [1]. In Section 3, we address the simplest case of a single mode: an harmonic oscillator or a two-level system. These serve as the first illustration of our theme that the 2PS reveals a considerably richer picture of the dynamics of emission for an otherwise identical 1PS lineshape, depending on both the mode statistics or the nature of its broadening (homogeneous or inhomogeneous). In Section 4, we analyse the more complex dynamics of coupled modes. Here again, even though the low-pumping regime results in identical 1PS lineshapes, the coupling of two harmonic oscillators, two two-level systems, or one harmonic oscillator to a two-level system (Jaynes-Cummings model), show distinct and peculiar features in the 2PS, unravelling the different nature of the three couplings. In the latter cases, of great importance in cavity-QED, we thereby identify an entire class of two-photon emission processes, to which one is completely oblivious through one-photon spectroscopy. In Section 5, we address the high excitation regime of a two-level system, considering two ways of driving it coherently: with the standard semiclassical model [17], that is, through the coupling to a classical term in the Hamiltonian, describing an ideal laser, and with a fully-quantized model [18] where a very large number of photons builds up a coherent field in a cavity. Both provide the expected Mollow triplet lineshapes in the 1PS but only the full quantum model provides a physical 2PS. In Section 6, we conclude and discuss further directions.

Theory of one-and two-photon spectra
We have recently shown how to compute efficiently and easily two-photon correlations between photons of frequencies ω 1 and ω 2 through the coupling of two "sensors" to the open quantum system of interest [1]. In the limit of vanishing coupling of the sensors to the mode which photons correlations are to be measured, the sensors intensity-intensity cross-correlations recover exactly the frequency-resolved two-photon correlations defined in the literature through integrals that are, however, too cumbersome to compute directly in most cases [19]. In the case of zero delay between the photons, the normalised correlation function, or two-photon spectrum (2PS), is therefore obtained simply in terms of the population operators of the sensors n 1 and n 2 , as: where ε represents the coupling to the sensors, small enough not to disturb the system, and Γ 1 , Γ 2 are the sensors decay rates that provide the detectors linewidths or, in other terms, the respective widths of the frequency windows. The case of equal frequencies g Γ (ω, ω) measures the photon-statistics of photons detected in a frequency window of Lorentzian shape of linewidth Γ. This corresponds to applying a single filter or to the effect of the detectors resolution in a colour-blind measurement, and is thus also of prime interest.
In practical terms, the calculation can be performed in two alternative ways, as discussed in Ref. [1]. The first method relies in the actual inclusion of the sensors in the dynamics of the system, supplementing the Hamiltonian and the master equation with the corresponding terms. The main drawback is that the size of the Hilbert space is increased by a factor 2 2 (and the density matrix by a factor 4 2 ). However, the computational operations are simple. In one of the possible computational approaches to the problem, one merely needs to solve the steady state of the full system (a set of homogeneous linear equations), which only requires the inversion of a matrix [20]. An alternative method, presented in detail in the supplemental material of Ref. [1], consists in using the formal semi-analytical expressions for the populations and cross correlations between the sensors to leading order in the couplings, ε 2 and ε 4 respectively. Then, the relevant quantities n 1 , n 2 and n 1 n 2 can be obtained through these expressions (that are provided below) and the master equation of the system only with no increase in the numerical complexity. In order to implement this method, one first needs to determine a vector of the system operators whose average are needed to compute the coupling to the sensors. Calling O the annihilation operator of the mode of interest, which correlations one wishes to compute, such a vector can be written in the form v = (1, O, O † , O † O . . .) (with subsequent terms depending on the systems and its dynamics). The observables of interest for the sensors cross-correlations are specified by the mean values of this vector . A regression matrix M can be obtained from the master equation of the system [20] which rules the dynamics of these correlators according to ∂ t v = M v . The steady state (if it exists) is obtained in the limit of infinite times: v ss = lim t→∞ e Mt v(0) , for any initial state v(0). Next, one builds two re-ordering matrices, T ± , which, when acting on the vector v , introduce in all correlators an extra operator O † for T + and an extra operator O for T − , keeping normal order in each case. That is, . .), respectively. With these matrices, one obtains the populations of the sensors to leading order in ε, as the first element [. . .] 1 of the vector: Regarded as a function of the frequency ω j , this produces the 1PS as measured by a detector of linewidth Γ j . Similarly, the cross-correlations that provide the 2PS is given by: where [1 ↔ 2] means the interchange of sensors 1 and 2, that is, permuting ω 1 ↔ ω 2 and Γ 1 ↔ Γ 2 everywhere. The advantage of this second method is twofold. First, it is very useful when the Hilbert space is small as it can lead to closed-form analytical expressions, as we show in the next Section. Second, it may also be numerically advantageous as the matrices involved correspond to the original (smaller) Hilbert space of the bare system, in the absence of sensors. The computational price to pay is in a higher number of matrix operations (eleven different matrix inversions and numerous multiplications for n 1 n 2 ). In some situations, the advantage of a smaller Hilbert space dominates. In others, it is more convenient and straightforward to explicitly include the sensors and solve for the enlarged system, with a similar overall numerical efficiency. We used both methods indistinctively to obtain the results presented in the rest of this text.
From now on, the linewidths of the detectors will be taken equal, Γ 1 = Γ 2 = Γ, for the sake of simplicity. We will thus simply write g (2) Γ (ω 1 , ω 2 ). The width Γ of the frequency windows is a fundamental parameter that cannot be dispensed with, unlike the 1PS case where an ideal detector is usually assumed in theoretical works. In the 2PS, Γ provides the uncertainty in the frequency (Γ) and time of detection (1/Γ) [21], as required by the Heisenberg principle. In the limit of very broad filters (Γ → ∞) the standard (color-blind) second order coherence function at zero delay is recovered: lim Γ→∞ g (2) Γ (ω 1 , ω 2 ) = g (2) . The opposite limit of very narrow filters (Γ → 0), results in a systematic tendency regardless of the underlying system, namely photons are uncorrelated, and those with the same frequency are thermally bunched. This can be understood on physical grounds: in order for the detectors to provide high precision in frequency (Γ → 0), their interaction time with the system has to be long (1/Γ → ∞), so that the collected photons correspond to all the possible times in the dynamics. This leads to an apparent uncorrelated statistics: lim Γ→0 g (2) Γ (ω 1 , ω 2 ) = 1 provided that ω 1 = ω 2 . The limiting value becomes lim Γ→0 g (2) Γ (ω, ω) = 2 in the case of detection at equal frequencies as there are 2! ways that two indistinguishable photons can be collected by two different detectors. We refer to this correlation tendency to bunch Harmonic Oscillator Two-level system Two-level system Figure 1. (a) The one-photon spectrum (1PS) common to a single harmonic oscillator and a two-level system, and the corresponding two-photon spectra (2PS), g Γ O [O](ω 1 , ω 2 ), for (b) an harmonic oscillator and (c)-(d) a two-level system, in units of the total broadening Γ O . In (d), the emitter is subject to pure dephasing, γ φ = 0.75Γ σ , with Γ σ such that the 1PS remains the one shown in (a). Beyond the bunching on the diagonal, common to all 2PS, clear deviations are observed depending on the type of emitter: the harmonic oscillator is structureless while the two-level system concentrates its antibunching in a butterfly-shaped region of non-overlapping frequencies. Dephasing favours antibunching. The color code is: 0 blue; 1 white; 2 red.
when ω 1 ≈ ω 2 and Γ is small (as compared to the relevant linewidths in the system), as indistinguishability bunching. This effect has been observed experimentally in the filtering of a single mode laser [22].

Single mode emitters
Two-photon frequency-resolved spectroscopy is such a nascent field of optical characterization that even the trivial systems need to be investigated. Namely, our starting point is the free mode (or, in the quantum optics terminology, "single mode"), which Hamiltonian simply reads: There are two fundamental possibilities for the operator O, namely, it can be the annihilation operator of an harmonic oscillator, in which case we denote it O = a, or it can be that of a two-level system, in which case we denote it O = σ. Their respective quantum algebra is given by the bosonic commutation rule (aa † − a † a = 1) and the fermionic anticommutation rule (σσ † + σ † σ = 1). The dynamics of the density matrix ρ of the emitter-including decay (required to bridge with the external world where the measurement is performed) and a continuous incoherent pump (to populate the system)-is given by the master equation with the Lindblad super-operator The 1PS, provided by Eq. (2), is the same for both systems O = a, σ. It is plotted in Fig. 1(a). It is a Lorentzian lineshape with linewidth Γ O (inverse mode lifetime) given by Γ a = γ a − P a for the harmonic oscillator and Γ σ = γ σ + P σ for the two-level system [23].
The two-photon correlations, however, are dramatically different, since the harmonic oscillator displays bunching, g (2) [a] = 2, while the two-level system displays the exact opposite behaviour, antibunching, g (2) [σ] = 0 (the mode which statistics is measured is here denoted explicitly in square brackets to lift the ambiguity). The 2PS also exhibits structures in the frequency correlations that are qualitatively different. For such simple cases, Eq. (3) can be solved exactly and analytical formulas for the 2PS obtained: with a common expression regardless of the mode O = a, σ spelt out above, and a term g (2) Γ [O] specific to each case, given by: for the harmonic oscillator (note that it is frequency independent), and bỹ for the two-level system. The single filter case g Γ (ω, ω) shows that for the twolevel system, correlations are maximum when filtering the peak itself, ω = 0, and that a large overlap of the peak is needed to recover a good antibunching of the two-level system (the ideal result is recovered in the limit of broadband detection: (2) [σ] = 0). The photon-counting correlation of the filtered peak indeed reads g (2) Γ [σ](0, 0) = 2(Γ σ /Γ)/(3 + Γ σ /Γ), e.g., the filter linewidth must be over sixty-six times that of the emitter to reach the 1% level of accuracy, an awkward requirement in practice. For this reason, post-processing of the coincidences through deconvolution of the spectral response is important to characterise fairly a quantum emitter.
For the harmonic oscillator, on the other hand, g Γ [a](ω, ω) = 2 regardless of the filtering window Γ and of which part of the peak is detected (maximum or any point in the tail). This is a manifestation of the classical character of the harmonic oscillator: it has no local information, any part behaves like the whole and one cannot tell apart from a filtered window whether the information is of a microscopic or macroscopic nature. In contrast, the two-level system has an energy scale and the information is localised in the energy window proper to the dynamics of the emitter.
The full 2PS of these systems further reveals such fundamental aspects of these emitters. In Fig. 1, we plot Eq. (6) in the case where the detector linewidth matches that of the 1PS peak, Γ = Γ O . In this and all subsequent density plots of the 2PS, we use a color code where red corresponds to bunching (g (2) 2), blue to antibunching (g (2) ≈ 0) and white to non-correlated emission (g (2) = 1). We are more concerned in this text with the qualitative patterns that emerge in two-photon spectroscopy than quantitative results (the optimisation of given correlations is straightforward but would bring us to a too voluminous discussion).
The common characteristic between the harmonic oscillator and the two-level system is the indistinguishability bunching line, for frequencies that are indeed indistinguishable within the detector linewidth |ω 1 − ω 2 | < Γ. This accounts for the diagonal feature on all plots. Apart from that, the harmonic oscillator correlations lack any structure and are always above one, corresponding to the expected bunching, whereas the two-level system correlations mainly assume values below one, the expected antibunching, furthermore in a nontrivial configuration of detection C with different frequencies that attempt to maximise the overlap with the 1PS without entering the indistinguishability bunching region: The largest antibunching is thus obtained on the antidiagonal. In practical terms, for a given detector linewidth, it is therefore better to perform photon coincidences between the left and right elbows of the 1PS, with a slight shift of both windows away from the center, rather than between photons both coming from the central peak. One can easily compute where exactly to place the filters by evaluating g Γ [σ](ω, −ω) with Eqs. (6) and (8). For instance, in the case Γ = Γ σ , the minimum goes down to 2( This is the closest that the frequencies can get to the central one without suffering from the indistinguishability bunching. It is a considerable improvement, at the small cost of a reduced signal, on measuring the correlations from the central peak with a detector of the same linewidth, that provides an antibunching of 0.5 only. When both detection windows are far from the peak, the emitter looses its quantum character and behaves like its classical counterpart, exhibiting bunching of its photons on the diagonal and uncorrelated emission otherwise. We complete the study of the correlations of a single two-level system emitter by characterising the effect of inhomogeneous broadening. This is relevant, for instance, in semiconductor systems (quantum dots) where the solid state environment induces fluctuations and pure dephasing on the levels. In this context, filtered two-photon correlations have already shed light on the timescales and origin of fluctuations in quantum dots [10]. Theoretically, inhomogeneous broadening is introduced in the dynamics through another Lindblad term of the form [24,25]. Pure dephasing occurs at the rate γ φ , diminishing coherence in the system without affecting directly the populations. Its effect on the 1PS is merely to increase the linewidth of the peak, Γ σ = γ σ + P σ + γ φ . For this reason, it is typically difficult to measure a radiative lifetime from spectroscopy and time-resolved measurements are usually invoked for that purpose. For a given 1PS, however, the 2PS changes quantitatively, which allows radiative and non-radiative contributions to be measured directly from spectroscopy measurements.
The 2PS spectra without (c), and with (d), pure dephasing are compared in Fig. 1. While the statistics of the two-level system as a whole is not affected by γ φ , its 2PS structure is. Namely, pure dephasing extends the condition of antibunching (9) to a wider range of frequencies and enhances the anticorrelations. This counter-intuitive result is a manifestation in 2PS of the profitable effect of dephasing for some quantum correlations [26,27].

Coupled single mode emitters
The bare emitters already display a rich structure of two-photon correlations. We now show that the next simple system, that of two emitters coupled linearly, also exhibits two-photon correlations of unsuspected complexity. Two coupled singlemode emitters account for, or are at the heart of, a large gamut of quantum optical systems. We will consider three paradigmatic cases, subject of intense theoretical and experimental studies: two coupled harmonic oscillators [28,29], two coupled two-level systems [30,31,32,33] and the Jaynes-Cummings model [34,35,8,36,37,38] where one harmonic oscillator is coupled to a two-level system. All of them are described by the following Hamiltonian: where, as in the previous section, O 1,2 describes two modes that can be either an harmonic oscillator or a two-level system. The dynamics of the system is completed with the inclusion of decay and pumping through the master equation: In the rest of the text, for the sake of brevity, we will limit to the situation where the first mode O 1 is detected and the other mode O 2 is pumped (P O 1 = 0). This reduces the number of parameters and remains close to the experimental situation in cavity-QED where typically the emitter is incoherently pumped and the system is observed through the cavity emission. In Fig. 2, we show the 2PS of the O 1 -mode for different values of γ O 1 , from 4g (weak coupling) to 0.1g (strong coupling), fixing the lifetime of the second mode to a small value γ O 2 = 0.001g, also on the basis of typical experimental values. The detectors linewidths are set equal to that of the emitting mode, Γ = γ O 1 , as previously. In the linear regime, when P O 2 → 0, all the models converge to the same 1PS (first row), whereas their 2PS exhibit very different structures depending on the nature of the emitters. When the losses overcome the coupling, i.e., in the weak-coupling regime, the 1PS converges to that of the bare mode O 1 . In the case of the coupled harmonic oscillators and the coupled two-level systems, the 2PS has the form of the single mode, cf. Fig. 1(a) and (b), respectively. However, in the case of the Jaynes-Cummings model, although the cavity mode is detected, a statistics reminiscent of the two-level system is observed, even in weak-coupling. This hints at the value of such systems for quantum applications, such as a single-photon sources: the (classical) cavity enhances the emission and collects the light while still retaining the valuable quantum features of the (quantum) emitter. Note that while juxtaposing the two filters on both sides of the central peak was only an improvement in the case of the bare two-level system, it is here mandatory, as the central peak is strongly bunched.
As the quality of the coupling increases (γ O 1 → 0.1g), the system enters in the strong-coupling with new dressed states, or polaritons, emerging. We call p ± the annihilation operators for the upper (+) and lower polariton (−). The 1PS, still identical for the three coupling models, turns into the so-called "Rabi doublet", with splitting 2 between the polaritons (the splitting observed in PL is more complicated and is given in Ref. [39]). The 2PS, however, present, again, different features that we discuss in turns.

Coupled harmonic oscillators
The coupled harmonic oscillators, that we denote a and b, are characterised by bunched correlations between any two frequencies, as expected due to the bosonic character of both modes. Similarly to the case of a single harmonic oscillator, the color-blind statistics always remains thermal, for both the bare and dressed modes: g (2) [a] = g (2) [b] = g (2) [p + ] = g (2) [p − ] = 2. The cross correlations, however, depend on the system parameters: (at resonance) and is always between 1 and 2. In the strong-coupling regime, g (2) [a; b] → 1, i.e., the two modes become uncorrelated. The correlations between polaritons, g (2) [p + ; p − ], are also close to 1 as long as polaritons are well defined (well separated spectrally), due to the formation of their independent polaritonic dynamics.
Coming back to frequency-resolved correlations, the 2PS of the mode a (2nd row in Fig. 2) presents, in strong-coupling, two hyperbolic dips where g (2) Γ (ω 1 , ω 2 ) drops to its minimum value of 1. The hyperbolas, given by the equation corresponds to the detected frequencies (ω 1 , ω 2 ) hitting the polariton branches which, as a function of detuning ∆, read ω ± (∆) = ∆/2 ± g 2 + (∆/2) 2 . When one detects a photon at ω 1 , the classical system does not provide information on whether this photon comes from the main peak itself (where the dressed state is located and the mode emits predominantly) or from any other point on its tail. This photon is detected as if emitted by the polariton being at precisely this energy. It has thus the fitting correlation of the coupled oscillators with another photon detected at the energy of the other polariton with the same detuning as the would-be polariton accounting for the first photon. The 2PS thus evidence the existence of the full polariton dispersion via a change in the statistics, even though the modes are at resonance in the 1PS. This shows again the richer discriminative power of two-photon spectroscopy. The 2PS anticrossing is wider than the 1PS one given that the distance between vertices (±g, ∓g) is 2 √ 2g instead of 2g. The renormalisation by √ 2 is due to the fact that the information is carried at the two-photon level instead of one in the photoluminescence.
When the modes are detuned, the polariton correlation interferences still form an hyperbola, but its center, ω hyp , and vertex, V hyp , depend on the set of relevant parameters of the system, namely {P a , P b , γ a , γ b , ω a , ω b }, as shown in Fig. 3. In order to gain further insight into this dependence, let us assume that polaritons are well defined and, therefore, the 1PS is essentially the sum of two Lorentzian peaks associated with each polariton, S[a](ω) ≈ L a − (ω) + L a + (ω). We can neglect the dispersive contributions to the spectrum [40] thanks to the negligible overlap between these peaks and assume complete uncorrelation between the polaritons, g (2) [p + ; p − ] ≈ 1. In a single particle picture, one can define a frequency dependent operator for mode a, α(ω) = L a − (ω)p − + L a + (ω)p + , using the relative spectral weights, that take into account the quantum state of the system through the coefficients l a ± [40]. The parameters γ ± and ω ± are the full-width half-maximum and positions of the polariton modes, p ± . The 1PS is well approximated by α † (ω)α(ω) . In a similar way, α † (ω 1 )α † (ω 2 )α(ω 2 )α(ω 1 ) , captures some of the features of the exact results in Fig. 3. The correlations computed in this way are less well reproduced than the 1PS as this approach neglects the multi-photon dynamics. However, due to the linear nature of the system, they produce interferences between the well defined polaritons with the same hyperbola, Eq. (14), specified by: In Fig. 3, we exemplify how the quantum state affects the anticrossing of the 2PS for a detuned situation, where ω − ≈ ω a = 0 and ω + ≈ ω b = 5g. For simplicity, only one mode is pumped in the linear regime (P a = 0, P b γ b ). The total decay of the system is fixed (γ a + γ b = 0.1g) whereas the relative value of the two rates is varied to alter the quantum state of the system. In the asymmetric case where γ b = 0, Fig.  3(a), the normalized weight is l a + ≈ 0 and, in accordance with Eq. (16), the center of the anticrossing ω hyp ≈ ω a with a vertex V hyp ≈ g. In the symmetric situation where Fig. 3(b), both weights are equal, l a + ≈ l a − , yielding a balanced position of the anticrossing between the two modes, ω hyp ≈ ωa+ω b 2 , and a larger splitting given by the same Rabi splitting as the 1PS, V hyp ≈ g 2 + (ω a − ω b ) 2 /4. Finally, in Fig. 3(c) where γ a = 0, the quantum state with l a − ≈ 0, yields ω hyp ≈ ω b and the same vertex as in the case of Fig. 3(a), V hyp ≈ g.

Coupled two-level systems
The coupled two-level systems (3rd row in Fig. 2) differs from the coupled harmonic oscillators in much the same respect than the two-level system does from the harmonic oscillator. There are now, unsurprisingly, regions of antibunching and a richer pattern of correlations.
As the system enters strong-coupling, a grid of two vertical and two horizontal lines of antibunching appear, corresponding to one frequency fixed at one of the peaks at ±R and the other frequency scanning through the 1PS. The stronger antibunching is obtained at their crossing. There is also the diagonal of indistinguishability bunching. Finally, a neat feature revealed by the 2PS is the strong superpoissonian antidiagonal, when ω 1 + ω 2 = 0. This is the first occurrence of an important and recurrent pattern in two-photon spectroscopy: these correlations correspond to direct two-photon relaxation from the doubly occupied state |1, 1 ⇒ |0, 0 , through a virtual (offresonant) intermediate state. As the energy of this auxiliary step is not constrained, the feature appears on an entire line. Only thanks to the coupling can one two-level system effectively emit two excitations at the same time (within the detector time window 1/Γ). We use the denomination of "leapfrog" to label this type of relaxation, as the intermediate rung is not populated and the system effectively jumps over it by emitting two photons simultaneously. It becomes more prominent as the detector width, Γ, decreases and the time uncertainty of the measurement increases. This is a fundamental process that is common to all quantum nonlinear systems. It supports the idea that the emission of two identical photons can be enhanced (and selected) by placing such systems in a high-Q cavity [41,42,43]. In the simple system that is the two-coupled two-level systems, it appears in isolation. We discuss it in more extent in the next paragraph where it appears along with similar processes, making sharper its own features.

Jaynes-Cummings model
In mixing the types of emitters, one harmonic oscillator (cavity) and one two-level system (atom, quantum dots, superconducting qubit. . . ), a much richer dynamics is unravelled by the 2PS (4th row in Fig. 2). In strong-coupling, a ladder-type level structure (shown in Fig. 4) arises thanks to the nonlinear splitting 2 ( √ ng) 2 − [(γ a − γ σ )/4] 2 between the states with n excitations [20]. The one-photon relaxations between these states give rise to a rich multiplet spectrum [44], and to a considerably richer set of two-photon relaxations, that we now discuss in details. The simplest and most straightforward type of two-photon correlations is the one that precisely follows from sequences of two one-photon relaxations, with each photon emitted from one state of the Jaynes-Cummings ladder to another. Depending on the type of relaxation, a given type of correlation is expected. For instance, at low pumping, detecting emission from the two types of polaritons in any given rung of the ladder leads to antibunching, as the de-excitation went one way or the other. Similarly, detecting emission in cascades from the ladder in consecutive events (as depicted in Fig. 4ii) leads to bunching. Such processes are the extensions to all possible rungs of the ladder of the correlations observed with the coupled two-level systems. In fact, when the system is not in very strong-coupling, and thus makes it difficult to resolve the dynamics of the higher rungs, the 2PS of the Jaynes-Cummings resembles that of the coupled two-level systems, with four antibunching dips appearing at ω 1 or ω 2 = ±R, corresponding to the situation depicted in Fig. 4i.
As the system gets deeper into strong-coupling, more such correlations from transitions between different rungs are revealed. They are given by ω 1 = E n ± E n−1 and ω 2 free (vertical lines) and their counterpart ω 1 ↔ ω 2 (horizontal lines). All together they constitute the grid of horizontal and vertical lines shown in black in Fig. 4(a). This pattern is neatly reproduced by the actual 2PS (see Fig. 5, for the best system) with enhancement/suppression of the correlations at every crossing point depending on whether the combination of frequencies corresponds to consecutive/non-consecutive processes [1].
There is also the now familiar diagonal of indistinguishability bunching, in red in Fig. 4(b), and other diagonal and antidiagonal features, better resolved in the best systems. These correspond to emission through virtual states during the short interaction time with the detector, ∼ 1/Γ, and have, to the best of our knowledge, not been reported up to now. The two antidiagonal lines at ω 1 + ω 2 = 2ω a ± E 2 ≈ ± √ 2g (dashed green lines in Fig. 4(b)) are the leapfrog processes already encountered with the coupled two-level systems: a two-photon transition from the second rung directly to the ground state, transiting through a virtual intermediate state which unspecified energy allows this process to be detected at any point on a line. This process is depicted in Fig. 4iii. The correlation for such a direct two-photon emission is, naturally, bunched.
Increasing the excitation power, the higher rungs get populated and the additional leapfrog process |3± ⇒ |1± can be observed. The four satellite antidiagonal lines that result, satisfying ω 1 + ω 2 = 2ω a ± (E 3 − E 1 ) or ±(E 3 + E 1 ), are plotted in Fig. 4(b) and materialised in the 2PS in Fig. 5(d). These lines are more visible in some regions (of strong bunching) but upon closer inspection, a deviation of the statistics due to such processes is realized everywhere.
A feature not observed in the previous cases is the secondary diagonal lines at ω 1 −ω 2 = ±2R, that meets with the antibunching dips at the first rung resonances. They are plotted as dotted-dashed blue lines in Fig. 4(b). In the actual 2PS, these are globally weaker and more clearly distinguished in some regions only, like the leapfrog lines from the higher rungs, and are also realized everywhere if looked at closely. Unlike leapfrog processes or indistinguishability bunching, however, they do exhibit an antibunching tendency as compared to the surrounding region (even if the overall correlation remains bunched). They correspond to the anticorrelated emission from both polaritons of the first rung into the same virtual final state, as depicted in Fig. 4iv. These lines become more pronounced when the mode under observation is also the one that is excited (as stated previously, we consider the situation where one mode is excited and the other one is observed). In such a pumping configuration, strong polariton-to-virtual-state correlations also appear in the coupled two-level systems (not shown) while, as seen in Fig. 2, third row, they are absent in the other pumping-detection scheme. This suggests that the interference between relaxations from |+ and |− plays a role. They are thus not peculiar to the Jaynes-Cummings model but to any strongly coupled, nonlinear system. Higher excitation, makes higher order diagonal lines appear, corresponding to higher rung dressed states, such as ω 1 − ω 2 = ±2E n . Such diagonals from the second rungs are also visible in Fig. 5(d). In this figure, the effect of increasing pumping is shown to "melt" the pattern of correlations as more real states are populated and more features distributed over the 2PS.
All these results show that two-photon spectroscopy unfolds the mechanisms of relaxation that are not apparent in photoluminescence. This can be used to characterize a system or exploit quantum correlations.

Lasing and the Mollow triplet
We have just discussed the effect of increasing the excitation power in the case of the Jaynes-Cummings model: higher rungs get populated and give rise to additional relaxation processes, either from real states that involve photon by photon de-excitations or through two-photon emission that involve intermediate virtual states. Further increasing P σ may bring the Jaynes-Cummings system into the lasing regime [45,46] where the optical field becomes coherent and g (2) [a] → 1. The cavity 1PS reduces to a single line that narrows with increasing occupancy of the mode as γ L ≈ g 2 /(2γ a n a 2 ) [47], where n a ≈ P σ /(2γ a ). The inverse quantity, 1/γ L , corresponds to the coherence time of the laser. More interesting is the 1PS of the two-level system (or quantum emitter), which converges to a Mollow triplet [17], with some specificities of its own in the cavity QED version as compared to resonance fluorescence due to the incoherent nature of the excitation [18]. Therefore, in this discussion we analyse the emission directly from the two-level system; this is the mode to which the sensors are now connected.
A Mollow triplet arises when a two-level system is driven by intense laser light. The splitting from very high rungs of the Jaynes-Cummings ladder becomes homogeneous giving rise to similar dressed states, |+ = c |0 − s |1 and |− = s |0 + c |1 , with |0 the ground state and |1 the excited state of the two-level system. The

Coherent excitation Incoherent excitation
Narrow filters Full-peak filters Figure 6. Normalised 1PS (first row) and 2PS (second and third rows) for the Mollow triplet under coherent (left) and incoherent excitation (right). Parameters are chosen so that the sidebands are at the same positions Ω ± = ±6.16g and have the same broadening with P σ = 2g and n a = 9.5. Vertical/horizontal gridlines mark the positions of the upper (green), lower (orange) and central black (peaks), arising from the transitions between dressed states depicted in (c). In (d) and (e) the detector width, Γ = 3P σ /2, is wide enough to filter full spectral peaks. Clear correlations appear for the pairs of transitions in (f). In (g) and (h) the detector width, Γ = P σ /2, is smaller than the spectral peaks. In this case, the dominant feature is the leapfrog triplet of antidiagonal lines, corresponding the to two-photon de-excitation processes sketched in (i).
A standard way to describe the laser excitation is through the semiclassical Hamiltonian term H L = Ω L (σ + σ † ) that couples the classical term Ω L (describing a macroscopic laser field) and the quantum term (σ + σ † ) (describing the quantum emitter) [17]. In this approximation, the laser has an infinite coherence time (γ L = 0) with intensity n a = Ω 2 L /g 2 and Ω ± = ±2Ω L . This is a successful description of the experimental observations, that report the theoretical lineshape, shown in Fig. 6(a). A fully quantised picture for the one-emitter laser, such as the Jaynes-Cummings model under incoherent pumping, provides the counterpart lineshape in Fig. 6(b). Parameters were chosen so that n a coincide in both cases (to 9.5) and the three peaks match in position and broadening, with P σ = 2g for the Jaynes-Cummings. They are not expected to match exactly since in the former case the coherence of the laser is perfect and independent of the emitter, while in the latter case, the coherence of the lightfield is established by the two-level emitter itself. They do share, however, essentially the same phenomenology, with the fundamental difference of quantisation or not of the light-field between the two. In the lineshapes, this causes the important difference of a finite linewidth γ L for the Rayleigh scattering peak on top of the two-level system triplet in the Jaynes-Cummings model, while it is a Dirac δ function in the conventional Mollow triplet theory due to the infinite coherence time of the laser. We will see that, although this is a rather innocent departure in the 1PS, it has dramatic consequences for the photon statistics.
We now describe the 2PS of these two cases. Correlations between full, well separated peaks of the semiclassical Mollow triplet have been the subject of intense experimental [48,49,11] and theoretical [12,13,50,51,52,14,15,53,16] efforts (note that the Mollow triplet of a cavity QED system in strong-coupling under incoherent pumping has still not be reported experimentally at the time of writing). Recently, the full 2PS (more specifically a closely related generalised Mandel Q parameter) of resonance fluorescence was obtained [16] within the generating function formalism, via single-molecule photon counting statistics with spectral resolution. The 2PS for both cases are compared in Fig. 6. They are smoother than in the quantum regime, Fig. 2, and with patterns more dominantly aligned along diagonals and antidiagonals rather than horizontal and vertical lines, implying that the system is close to a macroscopic, classical state as opposed to being dominated by correlations between well defined real states. One can understand the main features of these plots again through a decomposition into two-photon dressed state transitions: Each of the 2 3 terms represents a possible deexcitation path and all together, they interfere destructively to give rise to the expected total antibunching σσ = 0. The first line corresponds to the four cascade transitions depicted in Fig. 6(f), that start from |+ : • (Ω + , 0): One would expect this combination to produce bunched or cascaded Γ [σ](0, 0) ≈ 3 at small Γ when the coherent scattering is a significant fraction of the filtered light. Vanishing Γ yields the statistics of the classical laser since the semiclassical zero linewidth cannot display indinstinguishability bunching. (b) Comparison of the central peak correlations, (0, 0) (solid red), in the case of an incoherently pumped two-level system in free space (dotted black) and in a cavity (solid black). The laser intensity here is lower: Ω ± = ±6.16g as in Fig. 6. The linear regime, Ω L γ σ , is plotted with a dotted red line as a reference in units of γ σ = g. When the photon field is quantized, the correct limit g photons but, as it corresponds to two different paths with amplitudes of probability with opposite signs (first two terms in Eq. (18)), destructive interference leads to g (2) Γ [σ](Ω + , 0) = 0. Such debunching effect occurs within the detector timescale 1/Γ [49,54].
• (Ω + , Ω − ): Their cascade configuration produces a strong bunching, with a well defined time order that depends on the detuning with the laser [13]. However, at resonance, destructive interference debunches again the statistics to g Γ [σ](Ω + , Ω + ) = 0. These ideal correlations are shown in Fig. 7(a) in the full peak detection region (in yellow), where Ω L Γ σ . With a less intense laser such as than in Fig. 6, the features are still visible but slightly tempered due to the overlap between peaks, as is the case in experiments [11].
Thanks to the ease of use of our general solution, we are able to compare the semiclassical model for the laser 6(d) with the Jaynes-Cummings model in the lasing regime 6(e). In these qualitatively similar plots, the interference effect described above is even more evident as it extends beyond the four points (Ω ± , 0), forming two blue rings of antibunched correlations. These interference rings are more clear in the strongly driven situation depicted in Fig. 8. Their origin, as for the leapfrogs, stems from the uncertainty introduced by including the detector physics, which extends the interference appearing at the points (Ω ± , 0) and (Ω ± , Ω ± ) to two circles. At the center of these interference rings, (Ω ± /2, Ω ± /2), two red spots of enhanced emission appear from the concentration of leapfrog outer lines in the Jaynes-Cummings model, Eqs. (22). With a narrower detector, leapfrogs are enhanced due to the longer uncertainty in the time of two-photon emission, as shown in Fig. 6(g) and (h). In the same way that the Jaynes-Cummings multiplets from high rungs of excitation converge to the Mollow triplet in the 1PS, the antidiagonal leapfrog lines in the 2PS converge to a leapfrog triplet at: as schematically depicted in Fig. 6(i). They are more pronounced in Fig. 6(g), as the perfect laser approximation, H L , generates dressed states homogeneously split for all intensities Ω L . If the detectors are narrower than the peaks, Γ < Γ σ , the interference effect described above disappears (see Fig. 7). Furthermore, as has been previously discussed, when Γ is smaller than any peak width, photons are detected from any point of the dynamics and appear as uncorrelated. For different frequencies this means lim Γ→0 g (2) Γ [σ](ω 1 , ω 2 ) = 1 while for equal frequencies, lim Γ→0 g (2) Γ [σ](ω, ω) = 2 due to the indistinguishability of the photons. In this limit, the statistics observed can no longer be attributed to the system, but to the detection process. The two-level system excited incoherently (with or without the cavity) always recovers this limit correctly as, in the Jaynes-Cummings model, the smallest width, γ L is still finite and sets a clear lower boundary (see Fig. 7(b)). In contrast, we have a totally different result under coherent excitation of a perfect laser whenever ω 1 + ω 2 = 2ω L = 0: lim Γ→0 g (2) Γ [σ](0, 0) = 1 and lim Γ→0 g (2) Γ [σ](Ω + , Ω − ) = 2. Due to the implicit assumption of zero-linewidth for the laser field, it is not possible in these cases for the detector linewidth to be thinner and reach the uncorrelated limit. Instead, the limit Γ → 0 isolates completely the Rayleigh peak so correlations become exactly 1, i.e., those of the laser. This occurs with a detector linewidth below Γ min , defined as that at which the coherent part of the filtered 1PS is as large as the incoherent part, where [σ](ω). In the limit of intense lasing, Γ min ≈ Γ 3 σ /(4Ω 2 L ). In the intermediate region, Γ min < Γ < Γ σ , correlations from the central peak are bunched (and equal to 3 in the ideal case), due to the mixture of scattered and emitted light.
In real lasers, γ L = 0 due to phase fluctuations. Non-monochromatic theoretical models, where the phase varies stochastically [55,56], take into account both the amplitude and speed of fluctuations, recovering not only a finite γ L for the laser mode [57,55] but also physical correlations at all limits [58]. The Jaynes-Cummings model or one-atom laser is free from this pathology as fluctuations are intrinsic to the dynamics, provided by the interplay between coherent exchange, incoherent pump and decay. However, not being fully coherent unless very far above threshold, one requires a very good system (γ a g) in a very intense and well defined lasing regime, where γ L γ a , Γ min , that is P σ √ 2g. Like for any other system, very broad filters Γ Ω + recover the full antibunching of the two-level system, therefore lim Γ→∞ g (2) Γ [σ](0, 0) = g (2) σ [σ] = 0. This is independent of the type of excitation. Under very low coherent excitation (dotted red line in Fig. 7(b)), where the 1PS is dominated by the coherent part, correlations monotonously go from 1 to 0 as This regime has been exploited to create an ultra-coherent single photon source, i.e., indistinguishable in frequency thanks to the inherited long laser coherence time [59,60,61,62,63]. It is interesting to note that filtering the coherent peak, Γ < γ σ , yields the coherent statistics of a laser and to the destruction of the antibunching even if the coherent part strongly dominates the 1PS. Only the full emission has the property of interfering destructively and provide photons one by one. Finally, the impact of detuning from the laser (in the semiclassical approximation) is shown in Fig. 8. In the absence of pure dephasing, the 1PS remains symmetric around the laser frequency (here ω L = 0). Interestingly, the 2PS does not, even if detuning is small as compared to the driving, as shown in (b). The asymmetry in the 2PS consists in the disappearance of the ring of interferences nearer to the two-level system frequency (in this case ω σ > ω L ). When detuning dominates and Ω + ≈ ω σ , as in (c), the region around this peak becomes more antibunched, as it correspond to a real transition in a two-level system. On the other hand, the region around the other sideband, (Ω − , Ω − ), becomes more bunched, as it corresponds to a virtual laser transition [12,13]. Two red lines cross at (Ω − , Ω − ), only interrupted by the blue interference ring. A strong bunching point emerges at (Ω + , Ω − ), accompanied by the typical cascade τ -dynamics, making the system suitable for heralded single photon emission [11].

Resonance
Large detuning (c) Small detuning (b) (a) Figure 8. 2PS for a two-level system under coherent excitation deep in the Mollow triplet regime (Ω + = (2Ω L ) 2 + (ω σ − ω L ) 2 = 300γ σ ) in the cases where (a) the laser is on resonance with the emitter (Ω L = 150γ σ , as in Fig. 7(a)) and (b) the laser is slightly detuned (Ω L = 100γ σ ) or (c) largely detuned (Ω L = 10γ σ ). In contrast with the 1PS, the 2PS is asymmetric with detuning. Also, in this regime where a highly coherent field is involved, a new type of pattern emerges in the correlations, namely, circles.
at τ = 0, include an average over the short-time dressed-state dynamics, within the detector time-scale, 1/Γ. Therefore, the interaction with the filter/detector can be considered a dephasing mechanism of the dressed-state dynamics, which breaks the symmetry of the 2PS even at τ = 0. This view is confirmed by the behaviour of the physical time-dependent 1PS [21], in the absence of pure dephasing, which is also asymmetric in general before reaching its symmetric steady state lineshape (not shown).

Conclusions and further directions
In conclusion, we have performed the first systematic investigation of the two-photon spectra (2PS) of a variety of quantum optical systems, using our recently developed formalism to compute frequency and time resolved N -photon correlations [1]. We have focused here on the case N = 2 of two-photon correlations at zero time delay. Thanks to our formalism, such investigations can be generalised to higher photon number and/or arbitrary time delays. We studied systems of increasing complexity, starting from the simplest possible case that is the quantum harmonic oscillator and proceeding with the two-level system and then all their possible combinations, namely, two coupled harmonic oscillators, two two-level systems and a mixture of both that amounts to the celebrated Jaynes-Cummings model. The latter, in contrast with our starting point, is an extremely rich and complicated system that stands as the pillar of cavity quantum electrodynamics. We have outlined what are the common features shared by all these systems and what are those specific to their underlying structure and dynamics. These results constitute the backbone for two-photon spectroscopy. The main finding can be summarised as follows: • There is a universal bunching of photons when filtering an emission line below its linewidth, regardless of its inherent statistics, due to the indistinguishable and uncorrelated arrival of photons (indistinguishability bunching). This manifests as a diagonal on all 2PS.
• An antibunching arises for a two-level system with a characteristic butterfly-shape that is reproduced by any two-level transition that can be resolved in isolation from a more complex level structure (see Ref. [64] for an example within a fourlevel system).
• The 2PS of various types of coupled systems differ greatly from each other even when their 1PS are identical, and within the same system, some symmetries of the 1PS are lifted in 2PS. This shows the much higher degree of characterization accessible via two-photon spectroscopy.
• The 2PS of coupled harmonic oscillators features an anticrossing (two hyperbolas) of uncorrelated emission corresponding to the independent and classical polaritons dynamics. Coupled two-level systems or the Jaynes-Cummings model on the other hand give rise to rich patterns of correlations permitting to investigate the underlying physical picture in terms of relaxation between quantum states.
• Such correlations are due to fundamental processes that can be identified through simple equations that locate them in the 2PS. They are i) leapfrog processes, i.e., two-photon emission jumping over an intermediate rung of the level structure through a virtual state of indeterminate frequency and ii) polariton-to-virtualstate anticorrelations, where the emission does not correspond to dressed states transitions but take into account the dynamical nature of the system, with the final state provided by a virtual state or fluctuation of the system.
• Leapfrog processes are particularly noteworthy for applications in quantum information processing. They are present in any nonlinear quantum system, as simple as an anharmonic oscillator. N -photon leapfrogs with N > 2 also take place, though the higher the N , the weaker the process. In the Jaynes-Cummings system, they are evidenced by higher order correlations, g Γ (ω 1 , . . . , ω N ), at the frequencies N i=1 ω i = N ω a ± (E n − E n−N ) or (22) = N ω a ± (E n + E n−N ) , but their further characterization is out of scope of this text.
• The intricate pattern of N -photon correlations that is formed by the combination of all the above processes in the Jaynes-Cummings model, that is neatly resolved in a system well into the strong-coupling regime, evolves when brought in the lasing regime into the qualitatively different 2PS of the Mollow triplet, that features new types of patterns such as circles rather than simply straight lines. For high rungs of excitation, n 1, the lines in Eq. (22) converge to the central antidiagonal, i ω i ≈ 0, while the lines in Eq. (23) agglomerate at the outer positions, i ω i ≈ ±2 √ ng. This formation of the Mollow triplet in the 1PS gives rise to leapfrog triplet in the 2PS, with the same origin as in the linear regime in terms of two-photon transitions between dressed states.
Beyond spelling out the dynamics of emission at the quantum level, the theory of two-photon spectroscopy also allows to address fundamental theoretical issues of the quantum formalism and illustrates the deep link between the quantum dynamics and the detection process. For instance, in the limit of an ideal detector, the semiclassical description of the Mollow triplet fails to recover the fundamental indistinguishability bunching. This is due to the artificial δ line of the scattering peak produced by the semiclassical approximation. This shortcoming can be solved by turning to a fully quantized theory or upgrading the semiclassical theory to get rid of the artifacts caused by a fixed phase.
Another conclusion from our study is that the dynamics of "real processes", involving photon by photon de-excitations, are related to the system parameters and less influenced by detection. The corresponding correlations are best resolved when the related peaks are separated and fully filtered. On the other hand, "virtual processes", involving virtual states such as the leapfrogs and polariton-to-virtual state interferences, happen within the time of interaction with the detector, 1/Γ, and conserving energy within its linewidth, Γ. Therefore, their correlations become more prominent using narrow filters.
We have thus amply demonstrated that two-photon spectroscopy unravels a rich two-photon dynamics and interference effects in open quantum systems by a precise disposition of filters of given resolutions. Such correlations can be further taken advantage of by considering finite time delays and/or optimising the frequency windows. Further interesting systems to apply two-photon spectroscopy are ultra strong coupling systems [65], closely spaced atoms in optical lattice [66], the biexciton two-photon emission in a quantum dot [42,64] or the dynamics of Bose-Einstein condensation [67].