Two-photon spectra of quantum emitters

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, so as to obtain a continuous distribution that describes the probability of emission of each photon at a given energy. An example of such a one-photon spectrum (1PS) S(ω) is given in figure 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.

Zoom In Zoom Out Reset image size

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 paper, we will focus on the simplest case where N = 2. There is ample literature on two-photon correlations [5]. Typically, photons are detected as a function of time (and/or other commuting variables such as space) in photon counting configurations such as the celebrated Hanbury 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 ω₁ and another at energy ω₂. 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 ω₁ detected at time t₁ and another of frequency ω₂ detected at time t₂, although we will limit here to simultaneous detection for brevity. This is a generalization of the conventional (colour 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–11], although its interpretation relies chiefly on physical grounds or numerical simulations. Theoretically, almost all the effort has been devoted to resonance fluorescence [12–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 paper, 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⁽²⁾_Γ(ω₁,ω₂). 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 of N-photon correlations and/or arbitrary time delays is laid down in [1]. In section 3, we address the simplest case of a single mode: a 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 and 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 (the Jaynes–Cummings model) shows distinct and peculiar features in the 2PS, betraying 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 semi-classical 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.

We have recently shown how to compute efficiently and easily two-photon correlations between photons of frequencies ω₁ and ω₂ 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 whose photon 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 normalized correlation function, or the two-photon spectrum (2PS), is therefore obtained simply in terms of the population operators of the sensors n₁ and n₂ as

$$\begin{equation} g_{\Gamma_1,\Gamma_2}^{(2)}(\omega_1,\omega_2)=\lim_{\varepsilon \rightarrow 0}\frac{\langle{n_1 n_2}\rangle}{\langle{n_1}\rangle\langle{n_2}\rangle}, \end{equation} \tag{ 1 }$$

where ε represents the coupling to the sensors, small enough not to disturb the system, and Γ₁, Γ₂ are the sensor decay rates that provide the detector linewidths or, in other words, 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 linewidth Γ. This corresponds to applying a single filter of Lorentzian shape (a so-called interference filter) or to the effect of the detectors resolution in a colour-blind measurement, and is thus of prime interest.

In practical terms, the calculation can be performed in two different ways, as discussed in [1]. The first method relies on the actual inclusion of the sensors in the dynamics of the system, supplementing the Hamiltonian and the master equation with corresponding terms. The main drawback is that the size of the Hilbert space is increased by a factor of 2² (and the density matrix by a factor of 4²). 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 [1], consists in using the formal semi-analytical expressions for the populations and cross correlations between the sensors to leading order in the couplings, ε² and ε⁴, respectively. Then, the relevant quantities 〈n₁〉, 〈n₂〉 and 〈n₁n₂〉 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 averages 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 $\mathbf {v}=(1,O,{O}^{\dagger },{O}^{\dagger }O,\ldots )$ (with subsequent terms depending on the system and its dynamics). The observables of interest for the sensors cross-correlations are specified by the mean values of this vector: $\langle \mathbf {v}\rangle =(1, \langle {O}\rangle , \langle {{O}^{\dagger }}\rangle ,\langle {{O}^{\dagger }O}\rangle ,\ldots )$. 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 the normal order in each case. That is, $\mathbf {T}_+\langle \mathbf {v}\rangle =(\langle {{O}^{\dagger }}\rangle ,\langle {{O}^{\dagger }O}\rangle ,\langle {({O}^{\dagger })^2}\rangle ,\langle {({O}^{\dagger })^2O}\rangle ,\ldots )$ and $\mathbf {T}_-\langle \mathbf {v}\rangle =(\langle {O}\rangle ,\langle {O^2}\rangle ,\langle {{O}^{\dagger }O}\rangle ,\langle {{O}^{\dagger }O^2}\rangle ,\ldots )$, respectively. With these matrices, one obtains the populations of the sensors to leading order in ε, as the first element [··· ]₁ of the vector:

$$\begin{equation} \langle{n_j}\rangle=\frac{2 \varepsilon^2}{\Gamma_j}\Re\left[\mathbf{T}_+ \frac{-1}{\mathbf{M}+[-\mathrm{i}\omega_j-\frac{\Gamma_j}{2}]\mathbf{1}}\mathbf{T}_-\langle\mathbf{v}^{\mathrm{ss}}\rangle\right]_1,\quad j=1,\,2. \end{equation} \tag{ 2 }$$

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 are given by

$$\begin{eqnarray} &&\fl \langle{n_1 n_2}\rangle=\frac{2 \varepsilon^4}{\Gamma_1+\Gamma_2}\Re \left[ \mathbf{T}_+ \frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_2-\Gamma_1-\frac{\Gamma_2}2)\mathbf{1}}\right.\nonumber\\ &&\fl\quad\times \left\{ \mathbf{T}_- \frac{-1}{\mathbf{M}-\Gamma_1\mathbf{1}} \left( \mathbf{T}_+ \frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_1-\frac{\Gamma_1}2)\mathbf{1}}\mathbf{T}_- + \mathbf{T}_- \frac{-1}{\mathbf{M}+(\mathrm{i}\omega_1-\frac{\Gamma_1}2)\mathbf{1}} \mathbf{T}_+ \right)\right.\nonumber\\ &&\fl\quad + \mathbf{T}_- \frac{-1}{\mathbf{M}+(\mathrm{i}\omega_1-\mathrm{i}\omega_2-\frac{\Gamma_1+\Gamma_2}2)\mathbf{1}} \left( \mathbf{T}_- \frac{-1}{\mathbf{M}+(\mathrm{i}\omega_1-\frac{\Gamma_1}2)\mathbf{1}} \mathbf{T}_+ + \mathbf{T}_+ \frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_2-\frac{\Gamma_2}2)\mathbf{1}}\mathbf{T}_- \right)\nonumber\\ &&\fl\quad \left.\left.+\mathbf{T}_+\frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_1-\mathrm{i}\omega_2-\frac{\Gamma_1+\Gamma_2}2)\mathbf{1}} \mathbf{T}_- \left(\frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_2-\frac{\Gamma_2}2)\mathbf{1}} +\frac{-1}{\mathbf{M}+(-\mathrm{i}\omega_1-\frac{\Gamma_1}2)\mathbf{1}} \right)\mathbf{T}_- \right\}\mathbf{v}^{\mathrm{ss}}\right]_1\nonumber\\ &&\fl\quad +[1\leftrightarrow 2], \end{eqnarray} \tag{ 3 }$$

where $[1\leftrightarrow 2]$ means the interchange of sensors 1 and 2, that is, permuting $\omega _1\leftrightarrow \omega _2$ and $\Gamma _1\leftrightarrow \Gamma _2$ everywhere.

The advantage of this second method is twofold. Firstly, 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. Secondly, 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 (11 different matrix inversions and numerous multiplications for 〈n₁n₂〉). 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 indiscriminantly to obtain the results presented in the rest of this paper.

From now on, the linewidths of the detectors will be taken equal, Γ₁ = Γ₂ = Γ, for the sake of simplicity. We will thus simply write g⁽²⁾_Γ(ω₁,ω₂). 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 (colour-blind) second-order coherence function at zero delay is recovered: lim_Γ→∞g⁽²⁾_Γ(ω₁,ω₂) = g⁽²⁾. 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_Γ→0g⁽²⁾_Γ(ω₁,ω₂) = 1 provided that ω₁ ≠ ω₂. The limiting value becomes lim_Γ→0g⁽²⁾_Γ(ω,ω) = 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 shall refer to this correlation tendency to bunch when ω₁ ≈ ω₂ 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].

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 quantum optics terminology, 'single mode') whose Hamiltonian simply reads as (ℏ = 1)

$$\begin{equation} H_{O}=\omega_O{O}^{\dagger}O. \end{equation} \tag{ 4 }$$

There are two fundamental possibilities for the operator O: it can be the annihilation operator of a 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 made) and a continuous incoherent pump (to populate the system)—is given by the master equation

$$\begin{equation} \partial_t\rho=\mathrm{i}[\rho,H_O]+\frac{\gamma_O}2\mathcal{L}_{O}(\rho)+\frac{P_O}2\mathcal{L}_{{O}^{\dagger}}(\rho) \end{equation} \tag{ 5 }$$

with the Lindblad super-operator $\mathcal {L}_{O}(\rho )=(2 O\rho {O}^{\dagger }-{O}^{\dagger }O\rho -\rho {O}^{\dagger }O)$.

The 1PS, provided by equation (2), is the same for both systems O = a,σ. It is plotted in figure 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⁽²⁾[a] = 2, while the two-level system displays the exact opposite behaviour, antibunching, g⁽²⁾[σ] = 0 (the mode whose statistics is measured appears in square brackets to lift the ambiguity). The 2PS also exhibits structures in the frequency correlations that are qualitatively different. For such simple cases, equation (3) can be solved exactly and analytical formulae for the 2PS obtained:

$$\begin{equation} g_{\Gamma}^{(2)}[O](\omega_1,\omega_2)=\frac{\Gamma_O}{(\Gamma+\Gamma_O)^2}\left[\Gamma_O+\frac{\Gamma^2 (2\Gamma +\Gamma_O)}{\Gamma^2+(\omega_1-\omega_2)^2}+\tilde g_{\Gamma}^{(2)}[O](\omega_1,\omega_2)\right] \end{equation}\noindent \tag{ 6 }$$

with a common expression regardless of the mode O = a,σ spelt out above, and a term $\tilde g_{\Gamma }^{(2)}[O]$ specific to each case, given by

$$\begin{equation} \tilde g_{\Gamma}^{(2)}[a]=2\Gamma\left(1+\frac{\Gamma}{\Gamma_a}\right) \end{equation}\noindent \tag{ 7 }$$

for the harmonic oscillator (note that it is frequency independent), and by

$$\begin{equation} \tilde g_{\Gamma}^{(2)}[\sigma](\omega_1,\omega_2)=4\Gamma-2\Gamma(2\Gamma+\Gamma_\sigma)\left( \frac{\frac{ 3\Gamma+\Gamma_\sigma}{2}}{(\frac{3\Gamma+\Gamma_\sigma}{2})^2+\omega_1^2} +\frac{ \frac{3\Gamma+\Gamma_\sigma}{2}}{(\frac{3\Gamma+\Gamma_\sigma}{2})^2+\omega_2^2} \right) \end{equation} \tag{ 8 }$$

for the two-level system. The single filter case g⁽²⁾_Γ(ω,ω) shows that for the two-level 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: lim_Γ→∞g⁽²⁾_Γ[σ](ω₁,ω₂) = g⁽²⁾[σ] = 0). The photon-counting correlation of the filtered peak indeed reads g⁽²⁾_Γ[σ](0,0) = 2(Γ_σ/Γ)/(3 + Γ_σ/Γ), e.g. the filter linewidth must be over 66 times that of the emitter to reach the 1% level of accuracy, an awkward requirement in practice. This is an effect of the indeterminacy in the photon arrival time due to a better knowledge of its frequency. For this reason, post-processing of the coincidences through deconvolution of the spectral response is important to characterize 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 microscopic or macroscopic nature. In contrast, the two-level system has an energy scale and the information is localized 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 figure 1, we plot equation (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 colour code where red corresponds to bunching (g⁽²⁾ ≳ 2), blue to antibunching (g⁽²⁾ ≈ 0) and white to non-correlated emission (g⁽²⁾ = 1). We are more interested, in this paper, in the qualitative patterns that emerge in two-photon spectroscopy than the quantitative results (optimization of given correlations is straightforward but would bring us to a too voluminous discussion).

The common characteristic of the harmonic oscillator and the two-level system is the indistinguishability bunching line, for frequencies that are indeed indistinguishable within the detector linewidth |ω₁ − ω₂| < Γ. This accounts for the diagonal feature on all plots. Apart from that, the harmonic oscillator correlations lack any structure and are always above 1, corresponding to the expected bunching, whereas the two-level system correlations mainly assume values below 1. The strongest antibunching is obtained in a non-trivial configuration of detection  $\mathcal {C}$ with different frequencies that attempt to maximize the overlap with the 1PS without entering the indistinguishability bunching region

$$\begin{equation} \mathcal{C}=((|\omega_1-\omega_\sigma|<\Gamma_\sigma)\vee(|\omega_2-\omega_\sigma|<\Gamma_\sigma))\wedge(|\omega_1-\omega_2|>\Gamma). \end{equation} \tag{ 9 }$$

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 centre, rather than between photons both coming from the central peak. One can easily compute where exactly to place the filters by evaluating g⁽²⁾_Γ[σ](ω,−ω) with equations (6) and (8). For instance, in the case Γ = Γ_σ, the minimum goes down to  $2(\sqrt {2}-1)/5\approx0.16$, for $\omega _1=-\omega _2=\sqrt {(15\sqrt {2}-4)/31}\Gamma _\sigma\approx 0.75\,\Gamma _\sigma$. This is the closest that the frequencies can get to the central one without suffering from 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 loses 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 characterizing 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 $\frac {\gamma _\phi }{2}\mathcal {L}_{{\sigma }^{\dagger }\sigma }(\rho )$ [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 pure dephasing (c) and with pure dephasing (d) are compared in figure 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].

The bare emitters already have a 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 single-mode emitters account for, or are at the heart of, a large gamut of quantum optical systems. We will consider three paradigmatic cases, the subjects of intense theoretical and experimental studies: two coupled harmonic oscillators [28, 29], two coupled two-level systems [30–33] and the Jaynes–Cummings model [8, 34–38] where one harmonic oscillator is coupled to a two-level system. All of them are described by the following Hamiltonian:

$$\begin{equation} H_C=\omega_{O_1} {O}^{\dagger}_1 O_1 +\omega_{O_2} {O}^{\dagger}_2 O_2 +g ({O}^{\dagger}_1 O_2 +{O}^{\dagger}_2 O_1), \end{equation} \tag{ 10 }$$

where, as in the previous section, O_1,2 describes two modes that can be either a 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

$$\begin{equation} \partial_t\rho=\mathrm{i}[\rho,H_C]+\left[\frac{\gamma_{O_1}}2\mathcal{L}_{O_1}+\frac{\gamma_{O_2}}2\mathcal{L}_{O_2}+\frac{P_{O_1}}{2}\mathcal{L}_{{O}^{\dagger}_1}+\frac{P_{O_2}}{2}\mathcal{L}_{{O}^{\dagger}_2}\right](\rho). \end{equation} \tag{ 11 }$$

In the rest of the paper, for the sake of brevity, we will limit to the situation where the first mode O₁ is detected and the other mode O₂ is pumped (P_O1 = 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 figure 2, we show the 2PS of the O₁-mode for different values of γ_O1, from 4g (weak coupling) to 0.1g (strong coupling), fixing the lifetime of the second mode to a small value γ_O2 = 0.001g, also on the basis of typical experimental values. In the linear regime, when P_O2 → 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₁. In the case of the coupled harmonic oscillators and the coupled two-level systems, the 2PS has the form of the single mode, cf figures 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 single-photon sources: the (classical) cavity enhances the emission and collects the light while still retaining the valuable quantum features of the (quantum) emitter.

Zoom In Zoom Out Reset image size

As the quality of the coupling increases (γ_O1 → 0.1g), the system enters the strong coupling with new dressed states, or polaritons, emerging. We call p_± the annihilation operators for the upper (+) and the lower polariton (−). The 1PS, still identical for the three coupling models, turns into the so-called 'Rabi doublet', with underlying splitting $2R=2\sqrt {g^2-[(\gamma _{O_1}-\gamma _{O_2})/4]^2}$ between the polaritons (the expression for the splitting observed in PL is more complicated and is given in [39]). The 2PS, however, present, again, different features that we discuss in turn.

4.1. Coupled harmonic oscillators

The coupled harmonic oscillators, which we denote as a and b, are characterized 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 colour-blind statistics always remains thermal, for both the bare and dressed modes: g⁽²⁾[a] = g⁽²⁾[b] = g⁽²⁾[p₊] = g⁽²⁾[p₋] = 2. The cross correlations, however, depend on the system parameters

$$\begin{equation} g^{(2)}[a;b]\equiv\frac{\langle{{a}^{\dagger}a{b}^{\dagger}b}\rangle}{\langle{n_a}\rangle\langle{n_b}\rangle}=2-\frac{P_a/\langle{n_a}\rangle+P_b/\langle{n_b}\rangle}{\Gamma_a+\Gamma_b} \end{equation}\noindent \tag{ 12 }$$

(at resonance) and is always between 1 and 2. In the strong-coupling regime, g⁽²⁾[a;b] → 1, i.e. the two modes become uncorrelated. The correlations between polaritons, g⁽²⁾[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 (second row in figure 2) presents, in the strong coupling, two hyperbolic dips where g⁽²⁾_Γ(ω₁,ω₂) drops to its minimum value of 1. The hyperbolas, given by the equation

$$\begin{equation} \omega_1 \omega_2=-g^2, \end{equation}\noindent \tag{ 13 }$$

correspond to the detected frequencies (ω₁,ω₂) hitting the polariton branches which, as a function of detuning Δ, read $\omega _\pm (\Delta )=\Delta /2\pm \sqrt {g^2+(\Delta /2)^2}$. When one detects a photon at ω₁, 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 a polariton that is 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 evidences the existence of 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\sqrt {2}g$ instead of 2g. The renormalization by $\sqrt {2}$ is due to the fact that the information is carried at the two-photon level instead of one in the 1PS.

When the modes are detuned, the polariton correlation interferences still form an hyperbola,

$$\begin{equation} (\omega_1-\omega_{\mathrm{hyp}}) (\omega_2-\omega_{\mathrm{hyp}})=-V_{\mathrm{hyp}}^2 \end{equation}\noindent \tag{ 14 }$$

but its centre, ω_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 figure 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 non-correlation between the polaritons, g⁽²⁾[p₊;p₋] ≈ 1. In a single-particle picture, one can define a frequency-dependent operator for mode a, $\alpha (\omega )=\sqrt {L^a_{-}(\omega )} p_{-}+\sqrt {L^a_{+}(\omega )} p_{+}$, using the relative spectral weights

$$\begin{equation} L^a_{\pm}(\omega)= \frac{l_{\pm}^a}{\pi}\frac{\gamma_{\pm}/2}{(\gamma_{\pm}/2)^2+(\omega-\omega_{\pm})^2} \end{equation}\noindent \tag{ 15 }$$

that take into account the quantum state of the system through the coefficients l^a_± [40]. The parameters γ_± and ω_± are the full-width at half-maximum and positions of the polariton modes, p_±. The 1PS is well approximated by 〈α^†(ω)α(ω)〉. In a similar way, 〈α^†(ω₁)α^†(ω₂)α(ω₂)α(ω₁)〉 captures some of the features of the exact results in figure 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, equation (14), specified by

$$\begin{equation} \omega_{\mathrm{hyp}}\approx\frac{l_{-}^a}{l_{-}^a+l_{+}^a}\omega_{+}+ \frac{l_{+}^a}{l_{-}^a+l_{+}^a}\omega_{-}, \end{equation}\noindent \tag{ 16 }$$

$$\begin{equation} V_{\mathrm{hyp}}\approx \Re\sqrt{g^2-\left(\frac{4l_{-}^a l_{+}^a}{(l_{-}^a+l_{+}^a)^2}\right)^2\left(\frac{\Gamma_a-\Gamma_b}{4}+\mathrm{i}\frac{\omega_a-\omega_b}{2}\right)^2}. \end{equation}\noindent \tag{ 17 }$$

In figure 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, figure 3(a), the normalized weight is l^a₊ ≈ 0 and, in accordance with equation (16), the centre of the anticrossing ω_hyp ≈ ω_a with a vertex V_hyp ≈ g. In the symmetric situation where γ_a = γ_b (figure 3(b)), both weights are equal, l^a₊ ≈ l^a₋, yielding a balanced position of the anticrossing between the two modes, $\omega _{\mathrm {hyp}}\approx \frac {\omega _a+\omega _b}{2}$, and a larger splitting given by the same Rabi splitting as the 1PS, $V_{\mathrm {hyp}}\approx \sqrt {g^2+(\omega _a-\omega _b)^2/4}$. Finally, in figure 3(c) where γ_a = 0, the quantum state with l^a₋ ≈ 0 yields ω_hyp ≈ ω_b and the same vertex as in the case of figure 3(a), V_hyp ≈ g.

Zoom In Zoom Out Reset image size

4.2. Coupled two-level systems

4.3. The Jaynes–Cummings model

In mixing the types of emitters, one harmonic oscillator (cavity) and one two-level system (atom, quantum dot, superconducting qubit, etc), a much richer dynamics is unravelled by the 2PS (fourth row in figure 2). In the strong coupling regime, a ladder-type level structure (shown in figure 4) arises thanks to the nonlinear splitting $2{\rm E}_n=2\sqrt {(\sqrt {n}g)^2-[(\gamma _a-\gamma _\sigma )/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 an extended set of two-photon relaxations, which we now discuss in detail.

Zoom In Zoom Out Reset image size

The simplest and most straightforward type of two-photon correlation is the one that corresponds to a sequence 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 the 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 figure 4(ii)) leads to bunching. Such processes are 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 a very strong coupling, which 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 ω₁ or ω₂ = ± R, corresponding to the situation depicted in figure 4(i).

As the system gets deeper into strong coupling, more such correlations from transitions between different rungs are revealed. They are given by ω₁ = E_n ± E_n−1 and ω₂ free (vertical lines) and their counterpart $\omega _1\leftrightarrow \omega _2$ (horizontal lines). All together, they constitute the grid of horizontal and vertical lines shown in black in figure 4(a). This pattern is neatly reproduced by the actual 2PS (see figure 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].

Zoom In Zoom Out Reset image size

There is also the now familiar diagonal of indistinguishability bunching, in red in figure 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 $\omega _1+\omega _2=2\omega _a\pm E_2\approx \pm \sqrt {2}g$ (dashed green lines in figure 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 whose unspecified energy allows this process to be detected at any point on a line. This process is depicted in figure 4(iii). The correlation for such a direct two-photon emission is, naturally, bunched.

Increasing the excitation power, the higher rungs get populated and an additional leapfrog process $\left | {3\pm }\right \rangle \rightrightarrows |{1\pm }\rangle$ can be observed. The four satellite antidiagonal lines that result, satisfying ω₁ + ω₂ = 2ω_a ± (E₃ − E₁) or ±(E₃ + E₁), are plotted in figure 4(b) and are materialized in the 2PS in figure 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 ω₁ − ω₂ = ± 2R, which meet with the antibunching dips at the first rung resonances. They are plotted as dotted-dashed blue lines in figure 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 with 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 figure 4(iv). 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 figure 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 ω₁ − ω₂ = ± 2E_n. Such diagonals from the second rungs are also visible in figure 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.

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 involves 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⁽²⁾[a] → 1. The cavity 1PS reduces to a single line that narrows with increasing occupancy of the mode as γ_L ≈ g²/(2γ_a〈n_a〉²) [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 ladder operator can be decomposed into four possible transitions, σ = |0〉〈1| = c²| + 〉〈 − | − cs| + 〉〈 + | + cs| − 〉〈 − | − s²| − 〉〈 + |. The first and last terms give rise to sidebands in the 1PS at positions $\Omega _\pm {\approx } \pm 2g\sqrt {\langle {n_a}\rangle }$ (at resonance), as shown schematically in figure 6(c).

Zoom In Zoom Out Reset image size

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〉 = Ω²_L/g² and Ω_± = ± 2Ω_L. This is a successful description of the experimental observations that report the theoretical lineshape shown in figure 6(a). A fully quantized picture of the one-emitter laser, such as the Jaynes–Cummings model under incoherent pumping, provides the counterpart lineshape in figure 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 light field is established by the two-level emitter itself. They do share, however, essentially the same phenomenology, with the fundamental difference of quantization 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 intensive experimental [11, 48, 49] and theoretical [12–16, 50–53] efforts (note that the Mollow triplet of a cavity QED system in strong coupling under incoherent pumping has not yet been reported experimentally at the time of writing). Recently, the full 2PS (more specifically a closely related generalized 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 figure 6. They are smoother than in the quantum regime (figure 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:

$$\begin{eqnarray} &&\fl \sigma\sigma=s^3c\left| {-}\right\rangle(\left\langle {+}\right|+\rangle-\left\langle {-}\right|-\rangle)\left\langle {+}\right|+c^2s^2\left| {+}\right\rangle(\left\langle {+}\right|+\rangle-\left\langle {-}\right|-\rangle)\left\langle {+}\right|\nonumber\\ &&-c^3s\left| {+}\right\rangle(\left\langle {+}\right|+\rangle-\left\langle {-}\right|-\rangle)\left\langle {-}\right|-c^2s^2\left| {-}\right\rangle(\left\langle {+}\right|+\rangle-\left\langle {-}\right|-\rangle)\left\langle {-}\right|. \end{eqnarray} \tag{ 18 }$$

Each of the 2³ 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 figure 6(f), which start from | + 〉:

• (Ω₊,0): one would expect this combination to produce bunched or cascaded photons but, as it corresponds to two different paths with amplitudes of probability with opposite signs (first two terms in equation (18)), destructive interference leads to g⁽²⁾_Γ[σ](Ω₊,0) = 0. Such an 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⁽²⁾_Γ[σ](Ω₊,Ω₋) = 1.
• (0, 0): a similar situation of bunching and interference leads to g⁽²⁾_Γ[σ](0,0) = 1.
• (Ω₊,Ω₊): photons from the same sideband do not form a cascade so their simultaneous emission is suppressed (despite the indistinguishability bunching): g⁽²⁾_Γ[σ](Ω₊,Ω₊) = 0.

These ideal correlations are shown in figure 7(a) in the full peak detection region (in yellow), where Ω_L ≫ Γ_σ. With a less intense laser such as that in figure 6, the features are still visible but slightly tempered due to the overlap between peaks, as is the case in experiments [11].

Zoom In Zoom Out Reset image size

Thanks to the ease of use of our general solution, we are able to compare the semiclassical model for the laser figure 6(d) with the Jaynes–Cummings model in the lasing regime figure 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 clearer in the strongly driven situation depicted in figure 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 centre 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, equations (22). With a narrower detector, leapfrogs are enhanced due to the longer uncertainty in the time of two-photon emission, as shown in figures 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

$$\begin{equation} \omega_1+\omega_2 \approx 0, \quad \omega_1+\omega_2 \approx \Omega_\pm \end{equation} \tag{ 19 }$$

as schematically depicted in figure 6(i). They are more pronounced in figure 6(g), as the perfect laser approximation, H_L, generates dressed states homogeneously split for all intensities Ω_L.

Zoom In Zoom Out Reset image size

If the detectors are narrower than the peaks, Γ < Γ_σ, the interference effect described above disappears (see figure 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 that lim_Γ→0g⁽²⁾_Γ[σ](ω₁,ω₂) = 1, while for equal frequencies, lim_Γ→0g⁽²⁾_Γ[σ](ω,ω) = 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 figure 7(b)). In contrast, we have a totally different result under coherent excitation of a perfect laser whenever ω₁ + ω₂ = 2ω_L = 0: lim_Γ→0g⁽²⁾_Γ[σ](0,0) = 1 and lim_Γ→0g⁽²⁾_Γ[σ](Ω₊,Ω₋) = 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 completely isolates 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

$$\begin{equation} S^{\mathrm{coh}}_{\Gamma_{\mathrm{min}} }[\sigma](0)=S^{\mathrm{incoh}}_{\Gamma_{\mathrm{min}} }[\sigma] (0), \end{equation}\noindent \tag{ 20 }$$

where S_Γ[σ](ω) = S^coh_Γ[σ](ω) + S^incoh_Γ[σ](ω). In the limit of intense lasing, Γ_min ≈ Γ³_σ/(4Ω²_L). In the intermediate region, Γ_min < Γ < Γ_σ, correlations from the central peak are bunched (and are 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 [55, 57] 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_\sigma {\gg } \sqrt {2}g$.

Like for any other system, very broad filters Γ ≫ Ω₊ recover the full antibunching of the two-level system, therefore lim_Γ→∞g⁽²⁾_Γ[σ](0,0) = g⁽²⁾_σ[σ] = 0. This is independent of the type of excitation. Under very low coherent excitation (dotted red line in figure 7(b)), where the 1PS is dominated by the coherent part, correlations monotonically go from 1 to 0 as

$$\begin{equation} \lim_{\Omega_{\mathrm{L}}\rightarrow 0} g_{\Gamma}^{(2)}[\sigma] (0,0)=\left(\frac{\gamma_\sigma}{\gamma_\sigma+\Gamma}\right)^2. \end{equation}\noindent \tag{ 21 }$$

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–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 providing photons one by one.

Finally, the impact of detuning from the laser (in the semiclassical approximation) is shown in figure 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 interference closer 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 corresponds 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].

One can understand the asymmetry in the 2PS with detuning, g⁽²⁾_Γ[σ](Ω₊,Ω₊) ≠ g⁽²⁾_Γ[σ](Ω₋,Ω₋), as a manifestation of the dissimilar two-photon dynamics of the dressed states, namely, the τ-dynamics of the colour-blind correlations corresponding to each transition: g⁽²⁾[| − 〉〈 + |](τ) and g⁽²⁾[| + 〉〈 − |](τ). These two functions are identical in the limiting τ-values, equal to 0 at τ = 0 and to 1 at τ → ∞, but evolve differently at intermediate τ. In the limit of high excitation and full filtering of each peak, such dressed state two-photon dynamics is overall well mapped by the frequency-resolved correlations, g⁽²⁾_Γ[σ](Ω₊,Ω₊,τ) and g⁽²⁾_Γ[σ](Ω₋,Ω₋,τ), as long as τ > 1/Γ. However, as we have shown, the τ = 0 value is not well mapped, being larger than zero when the detection is included (spoiling perfect antibunching). The frequency-resolved functions 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).

In conclusion, we have performed the first systematic investigation of the 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 generalized 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 summarized 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 [64] for an example within a four-level 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 us 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-virtual-state anticorrelations, where the emission does not correspond to dressed state transitions but takes 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, although the higher the N, the weaker the process. In the Jaynes–Cummings system, they are evidenced by higher-order correlations, $g_{\Gamma }^{(N)}(\omega _1,\ldots ,\omega _N)$, at the frequencies
  $$\begin{equation} \sum_{i=1}^N \omega_i= N \omega_a \pm (E_n - E_{n-N})\quad \mathrm{or} \end{equation} \tag{ 22 }$$
  $$\begin{equation} = N \omega_a \pm (E_n + E_{n-N}) \end{equation} \tag{ 23 }$$
  but their further characterization is beyond the scope of this paper.
• The intricate pattern of N-photon correlations that is formed by the combination of all the above processes in the Jaynes–Cummings model, which 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, which features new types of patterns such as circles rather than simply straight lines. For high rungs of excitation, n ≫ 1, the lines in equation (22) converge to the central antidiagonal, $\sum _i \omega _i\approx 0$, while the lines in equation (23) agglomerate at the outer positions, $\sum _i \omega _i\approx \pm 2\sqrt {n}g$. 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 us 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 artefacts 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 are 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 resolution. Such correlations can be further taken advantage of by considering finite time delays and/or optimizing the frequency windows. Further interesting systems to apply two-photon spectroscopy are ultra-strong coupling systems [65], closely spaced atoms in an optical lattice [66], the biexciton two-photon emission in a quantum dot [42, 64] or the dynamics of Bose–Einstein condensation [67].

AGT acknowledges support from the FPU programme AP2008-00101 (MICINN); FPL from the RyC programme; CT from MAT2011-22997 (MINECO) and S-2009/ESP-1503 (CAM); MJH from the Emmy Noether project HA 5593/1-1 and from CRC 631 (DFG); and EdV from the Alexander von Humboldt Foundation.