Frequency-resolved Monte Carlo

TL;DR

This paper introduces a frequency-resolved Monte Carlo method adapted for quantum optics, enabling detailed simulation of photon emissions and correlations, including complex phenomena like the Mollow triplet and leapfrog processes.

Contribution

It extends quantum Monte Carlo simulations to frequency-resolved photon correlations, allowing for accurate modeling of photon emissions with known energies in quantum optical systems.

Findings

Simulation results agree with theoretical frequency-resolved photon correlation predictions.

Demonstrates leapfrog processes increase two-photon emission rates in the Mollow triplet regime.

Applicable to photon-counting experiments and quantum optics applications.

Abstract

We adapt the Quantum Monte Carlo method to the cascaded formalism of quantum optics, allowing us to simulate the emission of photons of known energy. Statistical processing of the photon clicks thus collected agrees with the theory of frequency-resolved photon correlations, extending the range of applications based on correlations of photons of prescribed energy, in particular those of a photon-counting character. We apply the technique to autocorrelations of photon streams from a two-level system under coherent and incoherent pumping, including the Mollow triplet regime where we demonstrate the direct manifestation of leapfrog processes in producing an increased rate of two-photon emission events.

Equations

$${\partial }_t\rho =i[\rho ,H_{\sigma }+{\omega }_{\xi }{\xi }^{†}\xi +\epsilon {\sigma }^{†}\xi +{\epsilon }^{\ast }\sigma {\xi }^{†}]+\frac{{\gamma }_{\sigma }}{2}{\mathcal{L}}_{\sigma }\rho +\frac{{\gamma }_{\xi }}{2}{\mathcal{L}}_{\xi }\rho ,$$

$${\partial }_t\vec{w}[\mu ,\nu ]=\{M+[(\mu -\nu )i{\omega }_{\xi }-(\mu +\nu )\frac{{\gamma }_{\xi }}{2}]1\}\vec{w}[\mu ,\nu ]+i\epsilon \mu T_{+}\vec{w}[\mu -1,\nu ]-i{\epsilon }^{\ast }\nu T_{-}\vec{w}[\mu ,\nu -1]+O({\epsilon }^2),$$

$$\vec{w}[\mu ,\nu ]=i|\epsilon |{\left\{M+[(\mu -\nu )i{\omega }_{\xi }-(\mu +\nu )\frac{{\gamma }_{\xi }}{2}]1\right\}}^{-1}\left(-e^{i\theta }\mu T_{+}\vec{w}[\mu -1,\nu ]+e^{-i\theta }\nu T_{-}\vec{w}[\mu ,\nu -1]\right)+O({\epsilon }^2),$$

$${\partial }_t\rho =i[\rho ,H_{\sigma }+{\omega }_{\xi }{\xi }^{†}\xi ]+\frac{{\gamma }_{\sigma }}{2}{\mathcal{L}}_{\sigma }\rho +\frac{{\gamma }_{\xi }}{2}{\mathcal{L}}_{\xi }\rho +\sqrt{\alpha {\gamma }_{\sigma }{\gamma }_{\xi }}\{[\sigma \rho ,{\xi }^{†}]+[\xi ,\rho {\sigma }^{†}]\}.$$

$$\frac{{\gamma }_{\sigma }}{2}{\mathcal{L}}_{\sigma }\rho +\frac{{\gamma }_{\xi }}{2}{\mathcal{L}}_{\xi }\rho +\sqrt{\alpha {\gamma }_{\sigma }{\gamma }_{\xi }}\{[\sigma \rho ,{\xi }^{†}]+[\xi ,\rho {\sigma }^{†}]\}=\frac{1}{2}{\mathcal{L}}_{\hat{o}}\rho +\frac{{\chi }_1{\gamma }_{\sigma }}{2}{\mathcal{L}}_{\sigma }\rho +\frac{{\chi }_2{\gamma }_{\xi }}{2}{\mathcal{L}}_{\xi }\rho +\frac{\sqrt{\alpha {\gamma }_{\sigma }{\gamma }_{\xi }}}{2}[\rho ,{\xi }^{†}\sigma -{\sigma }^{†}\xi ],$$

$${\partial }_t\vec{w}[\mu ,\nu ]=\{M+[(\mu -\nu )i{\omega }_{\xi }-(\mu +\nu )\frac{{\gamma }_{\xi }}{2}]1\}\vec{w}[\mu ,\nu ]-\sqrt{\alpha {\gamma }_{\sigma }{\gamma }_{\xi }}\{\mu T_{+}\vec{w}[\mu -1,\nu ]+\nu T_{-}\vec{w}[\mu ,\nu -1]\}.$$

$$|\psi (t+\delta t)⟩=\frac{c_k|\psi (t)⟩}{⟨\psi (t)|c_k^{†}{c}_k|\psi (t)⟩}\mathrm{or}|\psi (t+\delta t)⟩=\frac{\exp (i\tilde{H}\delta t)|\psi (t)⟩}{⟨\psi (t)|\exp (i{\tilde{H}}^{†}\delta t)\exp (i\tilde{H}\delta t)|\psi (t)⟩},$$

$$c_1=\sqrt{(1-{\chi }_1){\gamma }_{\sigma }}\sigma +\sqrt{(1-{\chi }_2){\gamma }_{\xi }}\xi ,c_2=\sqrt{{\chi }_1{\gamma }_{\sigma }}\sigma \mathrm{and}{c}_3=\sqrt{{\chi }_2{\gamma }_{\xi }}\xi ,$$

$$\tilde{H}$$

$$=H_{\sigma }+H_{\xi }-i\sqrt{\alpha {\gamma }_{\sigma }{\gamma }_{\xi }}{\xi }^{†}\sigma -\frac{i}{2}({\gamma }_{\sigma }{\sigma }^{†}\sigma +{\gamma }_{\xi }{\xi }^{†}\xi ).$$

$$H_{\sigma }={\omega }_{\sigma }{\sigma }^{†}\sigma$$

$$g^{(2)}(\tau )=1-\exp [-({\gamma }_{\sigma }+P_{\sigma })\tau ],$$

$$g^{(n)}=\frac{⟨s^{†n}{s}^n⟩}{⟨s^{†}s{⟩}^n}=\frac{⟨{\xi }^{†n}{\xi }^n⟩}{⟨{\xi }^{†}\xi {⟩}^n},$$

$${\gamma }_{\sigma }⟨s^{†}s⟩=\Gamma ⟨{\xi }^{†}\xi ⟩,$$

$$⟨s^{†n}{s}^n⟩={\left(\frac{\Gamma }{{\gamma }_{\sigma }}\right)}^n⟨{\xi }^{†n}{\xi }^n⟩,$$

$$⟨{\xi }^{†}\xi ⟩=\frac{P_{\sigma }(P_{\sigma }+\Gamma +{\gamma }_{\sigma })}{(P_{\sigma }+{\gamma }_{\sigma })((P_{\sigma }+\Gamma +{\gamma }_{\sigma }{)}^2+{\Delta }^2)}$$

$$\langle\xi^{\dagger}\xi\rangle=\Big{[}2\gamma_{01}\Omega^{2}(\gamma_{10}(\gamma_{11}^{2}+4\Delta^{2})^{2}(\gamma_{12}^{2}+4\Delta^{2})+4\Omega^{2}(\gamma_{10}\gamma_{11}^{2}\gamma_{12}\gamma_{32}+4(2\gamma_{10}^{3}+16\gamma_{10}^{2}\gamma_{01}+23\gamma_{10}\gamma_{01}^{2}+8\gamma_{01}^{3})\Delta^{2}\\ -16(\gamma_{10}-2\gamma_{01})\Delta^{4})+32\gamma_{11}(\gamma_{10}^{2}+4\Delta^{2})\Omega^{4})\Big{]}\bigg{/}\Big{[}(\gamma_{10}^{2}+4\Delta^{2})(\gamma_{11}^{2}+4\Delta^{2})(\gamma_{01}^{2}+4\Omega^{2})\times\\ \{\gamma_{10}^{4}+6\gamma_{10}^{3}\gamma_{01}+12\gamma_{10}\gamma_{01}(\gamma_{01}^{2}+2\Delta^{2}+4\Omega^{2})+\gamma_{10}^{2}[13\gamma_{01}^{2}+8(\Delta^{2}+2\Omega^{2})]+4[\gamma_{01}^{4}+4(\Delta^{2}-2\Omega^{2})^{2}+\gamma_{01}^{2}(5\Delta^{2}+8\Omega^{2})]\}\Big{]}$$

$${\omega }_1+{\omega }_2=\Delta ,\mathrm{with}\Delta =-{\Omega }_{+},0,{\Omega }_{+}.$$

$${\partial }_t\rho =i[\rho ,H_{\sigma }+\sum _{i=1,2}{\omega }_{{\xi }_i}{\xi }_i^{†}{\xi }_i]+\frac{{\gamma }_{\sigma }}{2}{\mathcal{L}}_{\sigma }\rho +\sum _{i=1,2}\left(\frac{{\gamma }_{{\xi }_i}}{2}{\mathcal{L}}_{{\xi }_i}\rho +\sqrt{{\alpha }_i{\gamma }_{\sigma }{\gamma }_{{\xi }_i}}\{[\sigma \rho ,{\xi }_i^{†}]+[{\xi }_i,\rho {\sigma }^{†}]\}\right)$$

$$c_1=\sqrt{{\chi }_0{\gamma }_{\sigma }}\sigma +\sqrt{(1-{\chi }_2){\gamma }_{{\xi }_1}}{\xi }_1,c_2=$$

$$c_4=\sqrt{{\chi }_2{\gamma }_{{\xi }_1}}$$

$$\tilde{H}=H_{\sigma }+H_{{\xi }_1}+H_{{\xi }_2}-i(\sqrt{{\alpha }_1{\gamma }_{\sigma }{\gamma }_{{\xi }_1}}{\xi }_1^{†}\sigma +\sqrt{{\alpha }_2{\gamma }_{\sigma }{\gamma }_{{\xi }_2}}{\xi }_2^{†}\sigma )-\frac{i}{2}({\gamma }_{\sigma }{\sigma }^{†}\sigma +{\gamma }_{{\xi }_1}{\xi }_1^{†}{\xi }_1+{\gamma }_{{\xi }_2}{\xi }_2^{†}{\xi }_2).$$

Full text

1 Introduction

The Monte Carlo method, of estimation by random sampling, was invented by Ulam as a practical way to estimate the chance of winning a game of solitaire, that a direct combinatorial approach had proved too challenging for the then convalescing scientist who was playing cards on a sick leave. [1] Since this had direct applications for the more “serious” problems tackled at Los Alamos (this was in 1946), such as the computation of neutron diffusion, the activity had to receive a codename, which came after the Monaco ward, famous for its casino. The technique indeed relies on chance by sampling randomly to get hold of a small but representative enough sample to describe a system. This is a surprisingly powerful technique that combines efficiency and accuracy, with applications in virtually all fields of human endeavours. [2]

In quantum physics, the Quantum Monte Carlo technique finds many ramifications in several fields. [3] In the context of interest in this work, that of quantum optics, several methods have been developed in the early 90s (see Ref. [4] for a review). Of these, the quantum jump approach [5, 6, 7, 8] (see Ref. [9] for an introduction) is particularly appealing as it links the wavefunction collapse to the emission of a photon. Assuming an ideal detector covering the full $4\pi$ solid angle surrounding the emitter, this allows to perform a computer experiment of photo-detections. From such “clicks” (as we will call a detected photon), one can for instance compute the Glauber correlation functions $g^{(n)}$ that measure the deviations of intensity correlations at the $n$th order from uncorrelated light, but one can also compute less easily accessible quantities such as exclusive probability densities, e.g., detecting the next photon at a time $\tau$ after one detection, with no other photon in between ($g^{(2)}$ assumes any photon rather than the next one), or distributions of time delays between nearest neighbours, probabilities to detect any given number or even configuration of photons in a time window, or any other type of binning “experiment”. Such Quantum Monte Carlo-generated photons have also been used to support the introduction of the $N$-photon “bundle”, [10] to distinguish the case of $N$-photon sources from strongly-correlated emission at the $N$ photon level.

In this text, we apply the quantum-jump Monte Carlo technique to the case of filtered emission, that is to say, as applied to a stream of photons going through an interference (i.e., Lorentzian) filter. Mean values for the correlators can be conveniently obtained with the theory of frequency-resolved photon-correlations. [11] This theory predicts strong correlations in frequency windows that had been neglected by both theorists and experimentalists until recently, [12, 13] as they lie far away from the luminescence peaks. Such correlations can clearly be turned into a resource and when technology will be mature to exploit photons as qubits, this aspect will certainly become compelling. As it should be useful to go beyond mean values and get access to time series for a variety of purposes, one could turn to the Quantum Monte Carlo technique applied to color-resolved photons. While the quantum Monte Carlo method has been used to compute the power spectrum as well as time-series of photon emission, [14, 15, 16, 17, 18, 19, 20] its combined use for both time and energy-resolved photons has, to the best of our knowledge, not been provided before in both a practical and exact form, as the few attempts in this direction [21, 22] have involved a cavity in the weak-coupling limit, which comes at the price of vanishing signal or approximate correlations if coupling of the filter is not weak enough. Monte Carlo methods have been used instead for their computational advantage or to access particular configurations such as the resonance fluorescence spectrum of an effective-three level system in a bright period exclusively (such systems are noted for their intermittent emission). In contrast, our approach allows to extract streams of photons from any frequency windows of a quantum source, using all the signal theoretically available and taking into account self-consistently the effect of its filtering, with an exact treatment of its effect on the correlations. This allows to revisit photon counting experiments with the added energy degree of freedom, that are already challenging without the frequency constrains.

In this text, we first prove that the technique is exact and, subsequently, we cover two cases for illustration. Namely, we show the effect of filtering a two-level system, and describe how this spoils the antibunching and quantum character of such sources in a practical context, although from a theoretical point of view, the saturated emitter is the brightest single-photon emitter, but of photons with wildly fluctuating frequencies (so not indinstinguishable). We illustrate how the loss of antibunching in time deviates notably from the single-exponential approximation used in the literature. [23]. We compare both the cases of coherent and incoherent excitations at low pumping. Then, at high pumping in the coherent case, thus bringing the problem in the Mollow triplet regime, we show how this turns a simple system into a versatile, tunable quantum source, with applications such as quantum spectroscopy [24, 25, 26] or photon sources with tunable statistics. [27] We will also address arguments [28] that imply that the strongly correlated emission is an artifact of normalization, which we will rebute by explicitly exhibiting these strongly correlated photons thanks to the Monte Carlo simulations. Before that, however, we briefly summarize the theory of frequency-resolved photon correlations and its main conclusions, which are to be found in greater details elsewhere, [11, 12, 24, 29, 22, 30, 31, 32, 33] and prove that the proposed Monte Carlo method scheme is Mathematically equivalent to this exact theory.

2 Theory

Although the process of filtering the light emitted from an optical source has a clear interpretation—the emitted photons are detected in a certain frequency window—its theoretical description used to be far from trivial. [34] With the introduction of the sensor method, [11] this became a straightforward task no more complicated than any old problem of computing quantum correlators, getting rid of all the complicated tasks of normal-ordering and time-integrals in spaces of many dimensions. The technique relies on coupling the system to sensors with strength $\varepsilon$ and taking the limit of vanishing coupling. In such a limit, it is enough to consider only two levels of the sensors as their populations remain $\ll 1$, without affecting the system’s dynamics. The computation of normalized correlators results in quantities that are $\varepsilon$ independent to first order and exact in the limit $|\varepsilon|\rightarrow 0$ (taking, in general, $\varepsilon\in\mathbb{C}$). These results are also absolute in the sense that they do not depend on detection efficiency or other details of the measurement, but characterize the source’s emission in given frequency windows. Alternatively, [35] such a coupling can also be made through the so-called “cascaded formalism”, [36, 37] that describes the dynamics of “detectors”, which are physical objects with a sizable coupling to the source (unlike sensors that have vanishing coupling) but that also do not alter the dynamics of the source, regardless of how strongly they are affected by it. Each method presents some advantages: the sensor method is straightforward to implement while the cascading formalism allows to characterize the detector’s dynamics beyond normalized correlations. An important difference with respect to computational cost is that while the sensor is typically described by a two-level system, the detector must be described by an harmonic oscillator, since it does get populated, whereas the sensor only acts as a probe in the limit of vanishing coupling. As a result, the sensor always provides exact results as a two-level system while the detector must be truncated high enough to provide a close-enough approximation, which depends on the dynamics of the system, and is therefore a tricky question. Furthermore, when considering cross-correlations, instead of $N$ two-level systems, one is dealing with $N$ harmonic oscillators and the problem becomes numerically forbidding, while the sensors’ Hilbert space scales as $2^{N}$ which is still tractable for $N$ much larger than anything that has been considered so far experimentally. For autocorrelations of the $N$th order with the sensor method, one can also use an harmonic oscillator truncated to $N$ excitations instead of $N$ two-level systems degenerate in frequency, which is also exact and with no need of checking convergenve for higher truncations, as is the case of the cascaded formalism. So in cases where correlations are requested rather than an actual signal, we believe the sensor method to be preferable, as it is both more efficient and more robust. In the present text, however, we specifically require a signal and will therefore turn to the cascaded formalism.

The mathematical equivalence of the two approaches for normalized autocorrelations can be established as follows. On the one hand, the sensor method “plugs” sensors to the dynamics. Formally, calling $\sigma$ the annihilation operator of a source and $\xi$ that of a sensor probing it, we can describe their joint dynamics by a Liouvillian equation

[TABLE]

where $H_{\sigma}$ is the Hamiltonian that describes the internal dynamics of the source, $\gamma_{\sigma}$ is the decay rate of the source, $\gamma_{\xi}$ is the decay rate of the sensor and $\mathcal{L}_{c}\rho\equiv 2c\rho{c^{\dagger}}-{c^{\dagger}}c\rho-\rho{c^{\dagger}}c$. The dynamics of an arbitrary operator $\xi^{\dagger\mu}\xi^{\nu}\sigma^{\dagger m}\sigma^{n}$ under the action of this Liouvillian is described with the notations of Ref. [11], by the equation (we set $\hbar=1$ along the paper)

[TABLE]

where $\mathbf{1}$ is the unit matrix, $T_{\pm}$ are normal-ordering superoperators for the $\sigma$ operators and $\vec{w}[\mu,\nu]$ is a vector of correlators for the $\mu$th and $\nu$th powers of the sensor operators $\xi$, $\xi^{\dagger}$ and spanning in normal order all powers of the $\sigma$, $\sigma^{\dagger}$ operators, i.e., $\vec{w}[\mu,\nu]\equiv(\langle\xi^{\dagger\mu}\xi^{\nu}\rangle,\langle\xi^{\dagger\mu}\xi^{\nu}\sigma\rangle,\langle\xi^{\dagger\mu}\xi^{\nu}{\sigma^{\dagger}}\rangle,\cdots,\langle\xi^{\dagger\mu}\xi^{\nu}\sigma^{\dagger m}\sigma^{n}\rangle,\cdots)^{T}$. The $O(\epsilon^{2})$ notation means that all other terms are of higher order than $\varepsilon$. The matrix $M$ provides the dynamics for the source, $\partial_{t}\vec{w}[0,0]=M\vec{w}[0,0]+O(\epsilon^{2})$, and is independent of the sensor at the lowest order in $\varepsilon$. At this stage, we do not assume any property of $\sigma$ or $\xi$, which could be bosonic (in which case $\mu$, $\nu$, $m$ and $n$ are unbounded) or fermionic (in which case $\mu$, $\nu$, $m$ and $n$ are 0 or 1). Equation (2) can be integrated, which yields

[TABLE]

where $\varepsilon=|\varepsilon|e^{i\theta}$. This in turn can be solved recursively, down to $w[0,0]$ where the equation self-truncates. This provides $\mu\nu+\mu+\nu$ terms, all with the same leading order of $\varepsilon^{\mu+\nu}$. More details on this derivation can be found in the Supplementary Material of Ref. [11]. The important point is that normalised correlators are ratios of components of $\vec{w}[\mu,\nu]$ with the $(\mu+\nu)$th power of components of $\vec{w}[1,1]$, themselves of order $\varepsilon$, so that such ratio cancel out $\varepsilon$ to leading order. While the higher-order terms do not have to cancel exactly, as is the case of the leading-order terms, they become negligible as the sensor coupling is made smaller, which does not preclude, however, a finite result for the normalized correlators, obtained from the $\varepsilon$ independent term, in contrast to unnormalised correlators that vanish with $\varepsilon$. Therefore, in the limit $\varepsilon\rightarrow 0$, the result becomes exact.

On the other hand, the cascaded formalism, which aims at exciting a target without affecting the source, provides a similar type of cancellation, although not restricted to vanishing coupling. From a causality point of view, it is clear that such a source/detector scenario where only one affects the other can be realized. The source that emitted a photon towards a detector may not even exist anymore by the time the detector is excited. This is achieved formally through interferences that cancel the back-action from the detector to the source. The master equation describing this asymmetric coupling reads

[TABLE]

The last three-terms of Eq. (4) can be re-written in the Lindblad form as

[TABLE]

where $\hat{o}=\sqrt{(1-\chi_{1})\gamma_{\sigma}}\,\sigma+\sqrt{(1-\chi_{2})\gamma_{\xi}}\,\xi$ is the joined decay operator of the whole system, source and detector, and the interpretation of the factor $\chi_{k}$ becomes that of factors that quantify the amount of signal that each part, source and detector, generates on its own and that the joined system generates as a whole. The detector, which must have a finite lifetime to couple to the source, thus also has an intrinsic frequency window with effect of filtering the emission it detects, whence the connection to the sensors formalism. The factor $\alpha=(1-\chi_{1})(1-\chi_{2})$, for $0\leq\chi_{1},\chi_{2}\leq 1$, takes into account that the source can have several decay channels. This is required for instance when only fluorescence is wanted without contamination from another source, e.g., a laser exciting it (experimentally this is typically achieved by detecting at right angle from the exciting beam).

Our proof proceeds by showing that $\xi^{\dagger\mu}\xi^{\nu}\sigma^{\dagger m}\sigma^{n}$ has the same equation as in the sensor formalism, by computing explicitly the equation for $\partial_{t}\vec{w}[\mu,\nu]$ in the cascaded formalism, Eq. (5). This reads, to all orders in the coupling in this case:

[TABLE]

Remarkably, this equation has the same form as Eq. (2) with $\varepsilon\rightarrow i\sqrt{\alpha\gamma_{\sigma}\gamma_{\xi}}$. Even though $\varepsilon$ is complex and a vanishing quantity in Eq. (2), with higher order corrections, and $\sqrt{\alpha\gamma_{\sigma}\gamma_{\xi}}$ it real and finite in Eq. (6), both methods provide exactly the same normalised correlators, as these coupling parameters enter in both the numerators and denominators with the same power and cancel out. The result becomes exact for vanishing coupling in the case of sensors and is exact in all cases with the cascaded formalism, regardless of their normalisation. Note as well that $\theta$, the phase of the coupling $\varepsilon$, has an effect on the dynamics only if the Lindblad equation features products of different operators in its dissipative terms, which is the case for the cascaded formalism with $\mathcal{L}_{\hat{o}}$ that brings cross terms of $\sigma$ and $\xi$. The sensor formalism, however, has no such joint decay emission and the phase of $\varepsilon$ does not play any role, so that $\varepsilon$ could have beeen set real. This achieves to prove the mathematical equivalence of the sensor method with the cascaded formalism for the computation of normalized correlators.

Since the sensor formalism has been shown[11] to be equivalent to normalized photon correlations according to photo-detection theory, [34] the above equivalence of the sensor and cascaded formalisms shows that applying the quantum Monte Carlo method to the detector[38] realises a sampling of the emission in the corresponding frequency windows, from which one can reconstruct the frequency-resolved photon correlations. That is to say, this allows us to simulate the photon emission with both time and energy information, which is what we are going to illustrate in the following. Note that with both Eqs. (2) and (6), any given correlator $\langle\xi^{\dagger\mu}\xi^{\nu}\rangle$ can be computed exactly (by recurrence) in terms of lower order ones $\langle\xi^{\dagger\mu^{\prime}}\xi^{\nu^{\prime}}\rangle$ only, with $\mu^{\prime}\leq\mu$ and $\nu^{\prime}\leq\nu$. This means that both methods can be applied using $N$ two-level systems as detectors at different frequencies, in order to compute cross correlations, or with a single harmonic oscillator truncated at $N$ excitations, to compute the $N$th order monocromatic autocorrelation function $g^{(N)}$. This is however not sufficient with the cascaded formalism for computing the density matrix (full state) of the detectors or for doing Monte Carlo simulations of the emission. In such cases, one must model the $N$ detectors as harmonic oscillators with a high enough truncation to provide converged results. The simulation is conveniently implemented through the quantum-jump approach. The dynamics of the system is thus described by a wavefunction $\left|\psi(t)\right\rangle$ that occasionally undergoes a process of “collapsing”, attributed to the emission of a photon, that one records in the simulation as a detector would register a click in an experiment. The collapse is decided in each infinitesimal time interval $\delta t\rightarrow 0$, where the evolution of the wavefunction is governed by two elements: a non-Hermitian Hamiltonian and random quantum jumps. In a system described by the master equation $\partial_{t}\rho=i[\rho,H]+(1/2)\sum_{k}\mathcal{L}_{c_{k}}\rho$, the non-Hermitian Hamiltonian is constructed as $\tilde{H}=H-(i/2)\sum_{k}{c_{k}^{\dagger}}c_{k}$, and the evolution of the wavefunction is given either by

[TABLE]

depending on whether the system undertook or not a quantum jump, respectively. The probability that this happened in the interval $\delta t$ due to the operator $c_{k}$ is proportional to the mean value of this operator, namely $p_{k}=\left\langle\psi\right|{c_{k}^{\dagger}}c_{k}\left|\psi\right\rangle\delta t$. For the system described by the master equation (4), the collapse operators are

[TABLE]

whereas the non-Hermitian Hamiltonian is given by

[TABLE]

The Monte Carlo approach within the cascaded formalism context has already been considered [37] but without connecting it to frequency filtering. Now that we have established this correspondence mathematically, we illustrate in the following how the clicks collected through Eqs. (4-5), in the frequency window determined by the detector, match indeed with the correlations predicted by the theory of frequency-resolved photon correlations. [11] Here, we apply this technique to the driven two-level system, under both coherent and incoherent emission, at low and large pumping. While some of the underlying physics has been published elsewhere, this will allow us to revisit it from another angle and provide additional results. We will consider both the cases of autocorrelations and cross-correlations. Although the same principles could be extended to even more than two detectors, we postpone such discussions and their further applications to future works.

2.1 The two-level system

The simplest system in quantum physics is the two-level system. We find it inspiring that it is still a Research problem, constantly offering new questions. Manifestly, in a quantum universe, the two-level system is as complicated as anything else. We first discuss a very basic problem, namely, the effect of filtering on photon emission from the Monte Carlo point of view. Starting with incoherent pumping, at rate $P_{\sigma}$, with Hamiltonian

[TABLE]

for a two-level system with free energy $\omega_{\sigma}$, and Liouvillian $\partial_{t}\rho=i[\rho,H_{\sigma}]+(\gamma_{\sigma}/2)\mathcal{L}_{\sigma}\rho+(P_{\sigma}/2)\mathcal{L}_{{\sigma^{\dagger}}}\rho$, with $\gamma_{\sigma}$ the inverse lifetime, one finds a simple enough dynamics of Glauber’s second order correlator:

[TABLE]

with, in particular, $g^{(2)}(0)=0$, that is, perfect antibunching. A (conventional) Monte Carlo simulation using the technique explained above, is shown in Fig. 1. The upper panel shows the fluctuations in the detection times of a million photons from such a source. As such, this realizes a random walk, similar to a random (Poissonian) process, and at large timescales there is nothing noticeable. On the short timescale, however, one can observe clear correlations of antibunching, as shown in the series of clicks indicated by blue ticks in Fig. 1(c). Namely, photons tend to repel each other and appear more orderly than if they would be uncorrelated, as is the case of the second series of photon detections, shown for comparison with black ticks in Fig. 1(d). The uncorrelated series exhibits the counter-intuitive “Poisson clumping” or “Poisson burst” effect, [39] made famous by von Bortkiewicz’s horse kicking casualties in the Prussian army and still of recurrent appearance in the medias as intuition repels the notion that a burst of accidents in, say, a hospital, is a natural random process rather than negligence. The strongly-correlated character of the two-level system emission becomes clear and compelling when computing intensity correlations $g^{(2)}(\tau)$ from the clicks, defined as the density of probability of finding two photons with a time difference $\tau$. Specifically, from the times of detection $t_{i}$, we compute $t_{i}-t_{j}$ for all $1\leq i\leq N$ with $N$ the total number of detected photons (here $N=10^{6}$) and compare the density of time differences to that from uncorrelated clicks with the same intensity. Note that in a typical experiment, a first photon starts a timer and a second stops it, and a distribution of the time difference between successive photons is used as a good approximation. In our case, we compute the exact correlations by collecting all the time differences within the correlation window of interest. This is shown for $|\tau|\leq 50/\gamma_{\sigma}$ in Fig. 1(e), left. One sees an overall plateau, indicating that photons have the same distribution for long-time separations as if they were emitted by a Poisson process (randomly). But one also observes a clear dip at $\tau\approx 0$, indicating that at such close distances, photons behave very differently than uncorrelated ones, namely, the occurrence of small time delays is strongly suppressed. This is better resolved in Fig. 1(e), right. Such a behaviour defines antibunching, $g^{(2)}(0)<g^{(2)}(\tau)$, with coincidences, i.e., simultaneous detection of two photons, less likely to occur than other closely spaced detections, with perfect suppression of coincidences when $g^{(2)}(0)=0$. Since these correlations wash out at long times, one has $\lim_{\tau\rightarrow\infty}g^{(2)}(\tau)=1$. The time it takes to reach this plateau is the second-order coherence time. We do not overlap these results of the Monte Carlo signal with the theory curve, Eq. (10), since, with one million points, it is exact to within the plot accuracy. Beside the statistical noise, that starts to be apparent for $\tau>1/(2\gamma_{\sigma})$, the Monte Carlo data provides a smooth curve in the window of strong correlations. In our simulation, the $\Delta t$ was $0.01/\gamma_{\sigma}$ and the binning size was taken twice as large, corresponding to the two closely-spaced vertical lines on the right panel of Fig. 1(e), bounding $g^{(2)}(0)$ from below due to this small uncertainty. With a binning size equal to the Monte Carlo timestep, one recovers the perfect antibunching at the origin, although on two grid points, so also producing a small error (the result would be perfect only in the limit of vanishing timesteps).

These results provide the background for our approach in the filtered signal. The general question is, what happens to the emitted photons if a filter is interposed on their way to the detector? This does not simply subtract a fraction, it also redistributes those that make it through, to provide them with possibly very different statistical properties, as we now discuss in more details.

3 Emission of a filtered two-level system

3.1 Incoherent excitation

The effect of a Lorentzian filter on the statistics of emission of an incoherently excited two-level system is shown in Fig. 2. The theory predicts thermalization and loss of antibunching with narrowing filtering. The exact way how this happens is discussed elsewhere, [23] and the theoretical result is shown on the density plot in Fig. 2 along with eight Monte-Carlo simulations of roughly 10 000 clicks each (25 000 for the narrowest filter in case $viii$). Extracts of the recorded clicks are shown, comparing them 1) in the same time window (black ticks), with effect of having much less clicks for narrower filters, and also 2) when rescaling the unit of time so that the intensities are the same (green ticks). In the latter case, one can compare the statistical distributions, and observe the transition from antibunched clicks ($i$) to thermal ones ($viii$) passing by auxiliary distributions. In the former case, one observes the characteristic antibunching, equally-spaced like distribution of a two-level emitter. In the latter case, one finds the wildly fluctuating thermal (or chaotic) light, with pronounced bunching in the form of long gaps of no emission followed by gusts of emission. This can be differentiated even with the naked eye from the Poisson distribution, whose tendency for “clumping” does not get as dramatic as the thermal case. One can follow the transition neatly from these various sets of clicks, passing by the case of uncorrelated light.

Since the isoline $g^{(2)}_{\Gamma}(\tau)=1$ is not straight  [23] (it is shown as a dotted line in the density plot of Fig. 2), the passage from antibunching to bunching does not transit through exactly uncorrelated (or coherent light), although the deviation is too small to be appreciated on a small sample. To observe such fine variations, one needs to acquire a large statistical ensemble and condense the correlations in a single object, such as $g^{(2)}_{\Gamma}$, as is shown in the eight panels at the bottom of Fig. 2. The case $v$ of close-to-uncorrelated light is also shown separately from the Monte Carlo data to reveal its fine structure. The other cases have a simpler shape of a dip that turns into a hump. The correlation time also changes dramatically, as is observable both from the density plot and the Monte Carlo histograms. As the emission thermalizes, its fluctuations occur on longer timescales. This is the reason for the increased noise in panels $vi$–$viii$. There, one should increase the binning and consider larger time windows, as shown in green for case $viii$ that assumes a binning of $\Delta t\gamma_{\sigma}=1$ instead of $0.1$ for the other cases, and plot the correlations in a time window $|t\gamma_{\sigma}|\leq 100$ instead of 10, as indicated on the respective axes, recovering the excellent agreement with theory displayed by the antibunched cases.

Now in possession of the statistical data, and with the insurance of its accuracy given its agreement with the theory, it is possible to undertake various types of analysis that would not be so straightforward theoretically, as has been described in the introduction. We will not go in this direction now and leave for future works and/or other colleagues (the statistical data is available on a public repository). Instead, we now turn to the case of coherent excitation, that also presents features of interest.

3.2 Coherent excitation

The case of filtering coherent is shown in Fig. 3. Here as well, there is thermalization, although this occurs in the case of strong-driving, [23] it is interesting to consider the effect of filtering and approach it from the Monte Carlo perspective. Taking one slice featuring these oscillations, we collect $10^{5}$ clicks, a small portion of which is shown as ticks at the bottom of Fig. 3. Computing the autocorrelations, we find indeed strong oscillations from a very good antibunching with steep bunching elbows, in agreement with the theory. This produces even more pronounced correlations in the photon-detection events, where the spacing appears more regular and between clumps of photons. As far as continuous streams are concerned, this suggests that such strongly-oscillating $g^{(2)}$ do in fact provide more ordered time series than the conventional antibunching of the type of Eq. (10). Such questions are however beyond the scope of the present text. We conclude this Section with further comments on the Heitler effect (coherence of the Rayleigh peak), that is broken at high pumping, but is eventually restored with narrow-enough filtering. First, regarding the emergence of a thermalization similar to that of incoherent driving, cf. Fig. 2, this is obtained when one enters the Mollow triplet regime. [40] In this case, luminescence has split into a triplet lineshape and, when filtering at resonance (as is the case here), one filters the central peak alone, which is known to correspond to the spontaneous emission of a photon that leaves the state of the dressed two-level system unchanged. [41] As such, the spontaneously emitted photons react to filtering in a similar way than the incoherently pumped two-level system, hence the observed bunching for narrowing filters linewidths. The similarity is only partial, however, as instead of thermalization, with $g^{(2)}_{\Gamma}(0)=2$, the transition is to a super-chaotic state, with $g^{(2)}_{\Gamma}(0)=3$ in the limit of infinite pumping [12] (for the parameters considered here, we find $\max_{\Gamma}g^{(2)}_{\Gamma}(0)\approx 2.2$). More strikingly, when filtering well within the central peak, one then isolates the Rayleigh ($\delta$) peak again and reverts to the low-pumping case, with the statistics becoming uncorrelated, as shown in Fig. 3. Large filtering windows, on the other hand, collect the emission from all three peaks and reproduce the Rabi oscillations, which is the case selected for the Monte Carlo sampling. We explore in more details the opportunities offered by the Mollow triplet in Section 4.

3.3 Effective quantum state

In this Section, we adapt the method of Zubizarreta Casalengua et al. [42] to get the photon distribution $p(k)$ of a frequency-resolved emitter, i.e., the probabilities for it to have $k$ photons, which are, equivalently, the diagonal elements of its effective density matrix. Namely, the frequency-filtered source is regarded as an effective source of its own. Since filtering typically turns the original source into one of another kind, the effective source gets attributed a new annihilation operator $s$ rather than the original $\sigma$. As we consider a photon source and we know that regardless of the original emitter, its filtering can produce an arbitrary number of photons, we assume that $s$ is bosonic. The detector and the effective source are of course related. On the one hand, the normalized correlators are identical, since the detector measures faithfully the correlations of the source:

[TABLE]

and as, on the other hand, there is conservation of energy (assuming an ideal detector):

[TABLE]

since the rate of emission from the effective source is also that of detection. We find from combining Eqs. (11) and (12) the relation between the unnormalized correlators:

[TABLE]

The statistics $p(n)$ of the photon emission from the effective source can now be obtained by inverting the relation $\langle s^{\dagger n}s^{n}\rangle=\sum_{k=n}^{N}({k!}/{(k-n)!})p(k)$ (with $N$ large enough) to provide the probabilities $p(k)$ for the effective source to have $k$ quanta of excitation, for integer $k$. In the cascaded formalism, the correlators of the detectors, $\langle\xi^{\dagger n}\xi^{n}\rangle$ can be computed from the master equation, and are source dependent. The expressions for the population of a detector being fed by an incoherently and coherently driven two-level system are given by López Carreño and Laussy [32] at resonance. Here, we provide a more general version for the detector at an arbitrary frequency. The incoherent case reads111Note that there is a typo in Ref. [32] with an extra factor 4; the correct result is as given here.

[TABLE]

and the coherent excitation reads

[TABLE]

with $\gamma_{ij}\equiv i\Gamma+j\gamma_{\sigma}$. From a Monte Carlo simulation in a time $T$, the detector population is obtained as the ratio between the total number of clicks recorded and $\gamma_{\sigma}T$. This allows to obtain the luminescence spectrum by scanning the detector in frequency (we do not show it but have checked it to be the case). One can now reconstruct the diagonal elements of the effective density matrix that, under an unspecified dynamics, is seen through the detector to yield the recorded photo-detection events. Since this can be achieved from (all) the Glauber correlators and the knowledge of the emitter’s mean population (known from the radiative lifetime), one can recover the effective $p(n)$ in this way. This allows us to access new classes of quantum steady states, tailored by frequency-filtering. We now illustrate how this takes shape in the case already discussed of filtered two-level system emission, starting with the case of incoherent excitation. This is shown in Figs. 4(a–b) for the filtered saturated two-level system, i.e., where the system is held in its excited state by very large pumping, $P_{\sigma}\gg\gamma_{\sigma}$, so that the density matrix reads $p(0)\rightarrow 0$ and $p(1)\rightarrow 1$. With Eq. (14) or from the detected clicks, one can compute the population and reconstruct this quantum state of the emitter, namely, the Fock state $p(n)=\delta_{n,1}$ ($\delta$ being the Kronecker function). The application of a filter turns the system from a two-level emitter to a source able to deliver more than one photon at a time, namely, for narrow enough filtering, $p(n)\approx(1-\theta)\theta^{n}$ for all $n$, with $\theta\approx 0.01$ for the narrowest filter considered here. We have lost two orders of magnitude for the population but the probability to observe two (resp. three) particles is 1% (resp. 0.01%) of that to observe one only, etc., which effectively shows how filtering a single-photon source turns it into a black body with nonzero probability to emit $n$ photons. Its population is smaller than without the filter, as the latter is rejecting some photons, but the statistical distribution of those that go through now corresponds to an altogether different quantum state. A similar situation occurs with coherent excitation (when not filtering so much as to isolate the Rayleigh peak), with, for the case shown in Fig. 4(c,d), $\theta\approx 0.025$. In both cases, one can see in this way at which point filtering prevents a single-photon source to emit non-classical states of light, [43] for instance by comparing $p(1)$ to $3\sqrt{3}/(4e)\approx 47.8\%$, the smallest probability above which a state is non-Gaussian. [44]

4 Tuneable statistics from the Mollow triplet

4.1 Auto-correlations

The Mollow regime that splits the luminescence into three lines provides in this way new natural spectral windows. One of the obvious questions this brings forward is: what is the statistics of the photons emitted by the three peaks? (is it the same as the total emission? One could imagine that all three peaks are emitting antibunched light since they ultimately originate from a two-level system. Our discussions so far in simpler systems should prepare us to find otherwise). The triplet structure, first computed exactly by Mollow [40] but without providing a physical picture as to its origin, can actually be well understood from a simple model, introduced by Cohen-Tannoudji et al.: the “dressed atom”. [45] In this model, the combined atom+laser is considered as a new entity, with a new structure of energy levels, shown in Fig. 5, and in which the transitions between the states account for the photoluminescence. On the basis of this picture, by considering two-photon transitions, one can foresee some correlations between the peaks. The transitions $\left|p\right\rangle\rightarrow\left|\pm\right\rangle\rightarrow\left|q\right\rangle$ for any $p,q=+$ or $-$, that go down the Mollow ladder, could be expected to result in bunching. In contrast, since one cannot chain in this way $\left|p\right\rangle\rightarrow\left|-\right\rangle$ and $\left|+\right\rangle\rightarrow\left|q\right\rangle$ or $\left|p\right\rangle\rightarrow\left|+\right\rangle$ and $\left|-\right\rangle\rightarrow\left|q\right\rangle$, one can expect antibunching. Inspection of all the combinations of two-photon relaxations “suggests” that:

• •

each side peak is antibunched (cases 8 and E in Fig. 5),

• •

the central peak comes with both bunching (2 and 4) and antibunching (9 and F),

• •

the central peak comes with both bunching (0, 1, 6, 7) and antibunching (A, B, C, D) with each side peaks,

• •

the side peaks are bunched together (3, 5).

(One could also go further and consider time-ordering and/or detuning.) Early calculations by Apanasevich and Kilin [46] and Cohen Tannoudji and Reynaud [45] confirmed the side-peaks antibunching in autocorrelation and their mutual bunching in cross-correlations. In a later work with more involved calculations, Schrama et al. [47] have shown that the central and side peaks feature no mutual correlations. This is due to interferences in the order of emission, that could nevertheless be linked to the co-existence of bunched and antibunched emission events. This shows that while intuition is strongly supported by the dressed-atom picture, it does not dispense from exact calculations for cases that could be ambiguous (and it confirms that there is indeed agreement with expectations for cases that cause no ambiguity, such as side-peaks emission).

With the theory of frequency-resolved photon correlations, [11] it is straightforward to compute such correlations exactly, like Mollow did for the luminescence spectrum, without referring to the dressed-states structure. This also allows to consider cases ouside their frequency windows, in fact, the complete landscape of two-photon correlations can be obtained. [12, 29] In the case of the Mollow triplet, [12, 48] it shows that the triplet structure reverberates at the two-photon level, through the apparition of a set of 3 hyperplanes, that obey the “leapfrog” equations

[TABLE]

(The same applies at the $N$-photon level [48]). The triplet structure comes, at any photon-number level, from the three possibilities to join the two dressed initial and final states. Note that while there are four transitions, $N$ is assumed sufficiently large for two transitions to be degenerate. The name of “leapfrog” comes from the fact that, at the $N$-photon level, transitions can occur by jumping over $N-1$ intermediate manifolds. This relaxes energy conservation for the individual photons and restricts it to the combined emission. In this sense, this is a ($N$ photons) version of Göpper Mayer’s two-photon processes. [49] that is however typically difficult to access (a notable exception is the planetary nebulae continuum [50]). The case $N=2$ is shown on the right of Fig. 5. In the case of resonance fluorescence, such correlations have been noted by Apanasevich and Kilin [51] for the $\left|\pm\right\rangle\rightrightarrows\left|\pm\right\rangle$ case (that is, overlooking the $\left|\pm\right\rangle\rightrightarrows\left|\mp\right\rangle$ counterparts, which are equally obvious with the dressed atom picture in mind). Interestingly, part of this school of researchers, who has produced noteworthy works on the problem of photon correlations, [52, 53, 54] has recently expressed some critics on this leapfrog picture, writing that [28] “the concept of the “leapfrog” processes is not justified”, that they “present an alternative explanation” based on “the unnormalized spectral correlation function” which is, they write, “a true measure of spectral correlations” and “which exhibits no signatures of the leapfrog transitions”. From their discussion, one thus understands that the production of strongly-correlated photons away from the peaks, that we predict, is an artifact due to normalization. In principle, one can indeed inflate the vacuum and create what one could regard as an artificial superbunching. This is not what happens with leapfrog emission, however, although according to these authors, nothing of interest takes place away from the peaks, what they illustrate by producing a two-photon spectrum remarkably featureless, in contrast to our two-photon spectrum that is rich from photon correlations flourishing away from the peaks. Their non-normalized spectrum is correct but, we believe, is also not interesting as it merely shows that first order processes smother second-order ones, as is however well-known and expected. We show, in contrast, that the scarce signal from higher-order processes has stronger correlations than those from first-order processes. This will be amply and vividly illustrated through Monte Carlo simulations below. One can also consider placing a cavity in this “featureless” region when not normalized, and observe how the system then keeps emitting strongly correlated photons but now dominating over the other first-order processes, [10, 55] which would not happen would the correlation be an artifact due to normalization. We will show below thanks to the frequency-resolved Monte Carlo simulations how one can anyway see the manifestation of leapfrog processes with the naked eye.

The rest of their discussion is only semantics, in which case we should clarify, as this is apparently needed, that the dressed-atom picture is, precisely, a “picture”, that is, an insightful mental representation that is helpful to visualize the basic mechanism at play, support the intuition and guide further inquiries. This does not preclude exact calculations based on the opaque equations of quantum mechanics. This is possibly why the Mollow triplet is named after Benjamin, not only for his seminal input but also in recognition of the exact expression, although the Cohen Tannoudji and Reynaud approximate picture is the one everybody has in mind when thinking about this problem. We combine both approaches, the sensor technique provides the exact result, while the leapfrog processes provide a physical representation à la Cohen Tannoudji et al.. Thus, in the same way that resonance fluorescence is not spontaneous emission from the dressed atom, the leapfrog emission is not, strictly speaking, spontaneous emission jumping over intermediate manifolds. This is, instead, a complicated process that involves the laser and the two-level system in a sequence of coherent absorption and emission. We have in fact shown [33] how even fluorescence in the low-driving regime does not consist of Rayleigh scattering events but form an intricate interference between emission and absorption, that powers the single-photon emission mechanism by suppressing the fat tails of the Lorentzian and turning the lineshape into a $t_{2}$ distribution instead. In the strongly nonlinear regime, a similar dynamics takes place but it becomes forbidding and certainly not even useful to apprehend the problem in these terms. Note that although inspired by the dressed-atom approach of the problem, ultimately, our computations are exact (and in full agreement with this physical picture). While the dressed atom has proven to be extremely fruitful for their purpose, we find it to be even more so in our $N$-dimensional case, [48] if not mandatory. It is, indeed, very easy based on this concept to understand why some configurations have less strong correlations than others, for instance, one \left|+\right\rangle\mathrel{\vbox{\offinterlineskip\hbox{\rightarrow}\vskip 0.43057pt\hbox{\rightarrow}\vskip 0.43057pt\hbox{\rightarrow}}}\left|+\right\rangle transition with one photon non-degenerate with the two others (of interest for photon-heralding of two-photon emission), is particularly weak. This is because it is resonant with the $\left|-\right\rangle\rightarrow\left|+\right\rangle\rightrightarrows\left|-\right\rangle$ transition, that breaks this channel of relaxation by interposing a real state in the three-photon emission. Configurations \left|+\right\rangle\mathrel{\vbox{\offinterlineskip\hbox{\rightarrow}\vskip 0.43057pt\hbox{\rightarrow}\vskip 0.43057pt\hbox{\rightarrow}}}\left|+\right\rangle whose energies do not intersect with real states, on the other hand, have very strong correlations and are more suitable for heralding purposes. It would seem difficult to make sense of these observable facts, that follow from exact computations, without the leapfrogging concept. In fact we can easily generalise them to arbitrary photon orders and guess which energies are to be avoided in a $N$ dimensional space to harness the best sought configuration of multi-photon emission, without having to undertake any actual computation. Is the concept therefore not justified and should one only be allowed infinite series of Feynman diagrams? We believe that the comments of Shatokhin et al. targetting the leapfrog picture bring very little to the discussion, if not in fact muddying it with confused statements and blurring their actual technical contributions that otherwise concur with our results, and of which we wish to remove nothing, as this approach has its merits.

Back to the general discussion. We show in Fig. 6 the statistics of clicks from photo-detection events of the Mollow triplet in frequency windows spanning from the central peak to the side peaks, including various other windows in between, in particular, the leapfrog window. Note that, here as well, the data is a single-detector observable, that is to say, the different streams shown are not correlated to each others as they have been obtained by the same detector in different runs of the experiment. It would require $5$ detectors to obtain the same result but with correct cross-correlations (this is beyond the scope of the present discussion that will go up to two detectors only, but is of course a topic of interest for applications). As we did in Fig. 2, we show both ticks in a given time window (in black) and with a rescaling of the unit of time so that their densities are equal (in green). Here as well, the relative emission rates mean that longer integration times are required when collecting away from the peaks. The gain in terms of correlation strengths, however, makes it worthwhile to focus on these regions of reduced emission, in a spirit akin to distillation: trading quantity for quality. The frequency windows have been chosen as they correspond to particular cases of interest:

1. $i$

Photons from the central peak. 2. $ii$

Case where $g^{(2)}_{\Gamma}(0)=2$ (usually attributed to thermal light). 3. $iii$

Photons from leapfrog emission. 4. $iv$

Case where $g^{(2)}_{\Gamma}(0)=1$ (usually attributed to coherent or uncorrelated light). 5. $v$

Photons from a side peak.

The central peak is partially thermalized, with a $g^{(2)}(\tau)$ that closely resembles the form of thermal fluctuations, $g^{(2)}(\tau)=1+\exp(-2|\tau|/\tau_{0})$. Upon closer inspection, however, this is an approximation as the exact solution presents small departures, in particular, a differentiable slope at the origin and small ripples that are thinly visible on the theory curve, that we keep separate from the Monte Carlo data for clarity (the quality of their agreement is shown in the rightmost column). Note that the dynamics of coherent driving of a two-level system is considerably more complicated than its incoherent counterpart and we could not find, so far, a general closed-form expression for $g^{(2)}_{\Gamma}(\tau)$ in this case for arbitrary frequencies. Applying the technique of effective-quantum state reconstruction from the correlators, described in Section 3.3, we find that the statistics $p(n)$ fits well with a cothermal distribution with $\approx 80\%$ of thermal emission and $\approx 20\%$ of uncorrelated emission. Overall, the emission of the central peak is thus well described by a mixture of thermal and uncorrelated light. It is, as such, not very interesting per se. Turning now to the second frequency window, $ii$, which features $g^{(2)}_{\Gamma}(0)=2$, characteristic of thermal emission, one can now see more clearly the deviation from the thermal paradigm, with bulging and tails deforming the correlation function. These are well reproduced by the Monte Carlo statistics and we let the reader decide if their statistical acuity lets them, on the small sample of clicks reported here, observe deviations from the thermal paradigm (cf. Fig. 2$viii$)

The most interesting window, $iii$, lies halfway between the central and side peak. This is the frequency at which, according to our interpretation of the theory, [12] two photons can make a leapfrog process from the state $\left|+\right\rangle$ in a given manifold to the state $\left|-\right\rangle$ two manifolds below, jumping over the intermediate manifold, cf. the rightmost transition in Fig. 5. These photons are strongly correlated in several ways. From a photo-detection point of view, they should arise as more occurrences of closely-spaced two-photon clicks than if the emission was uncorrelated. In particular, their rate of coincidences should increase, leading to $g^{(2)}_{\Gamma}(0)\gg 2$, or so-called superbunching. This is both predicted by the exact theory [11, 12] and observed in our Monte Carlo simulations, as seen in Fig. 6. Remarkably, even with as few as 9112 clicks collected in the numerical experiment, we can reconstruct a high-quality signal, revealing the fine details of its structure. Note as well that on the real-time series of clicks, out of the nine photons emitted, four came as two-photon bundles (the fifth and sixth clicks are so closely spaced as almost overlapping; other ticks are single-photon events). The small sample of clicks also shows strong ordering, combining equal spacing and gaps of no emission. While the latter is characteristic of thermal emission, the former is typically characteristic of antibunching. This combination can be seen as the selection through filtering of strongly correlated emission from the emitter, rather than tampering from the filters on the statistics: focusing to this frequency windows allows us to detect the two-photon emission events that occurs, from the dressed-atom picture, at this frequency. It would be rewarding to apply this technique to the filtered emission of a “bundler”, [10] a device that emits the majority, and in some regime, close to 100%, of its light as $N$-photon emission, and for which filtering has been shown to considerably boost the purity of the quantum emission. [56, 55] Also further photon-counting characterization would certainly be enlightening, and preliminary investigations show that the percentage of closely-spaced photons is over an order of magnitude higher in $iii$ than in the others at the exception of $ii$, as compared to which it is only about $3.8$ times larger. We leave further characterizations for future works, but provide a last compelling manifestation of leapfrog emission from the effective quantum state reconstruction approach, cf. Section 3.3. This highlights the frequency window $iii$ as the most dissimilar one as compared to the others, featuring a neat kink at the probability to have two photons, $p(2)$, showing the relative predominance of two-photon emission. Overall, this simulation makes it obvious that the emission in this frequency window suffers from no artifact of post-selection or normalization, but does indeed provide strongly correlated photon streams.

The fourth frequency window, $iv$, chosen for its $g^{(2)}_{\Gamma}(0)=1$ of uncorrelated emission, is also a case that shows strong departures at nonzero $\tau$ due to filtering. This is, here again, well captured by the Monte Carlo clicks and is visibly noticeable on the small sample, that features ordered clumps of uncorrelated clicks. With the last window, $v$, we come back to a case well studied in the literature, of antibunched emission, albeit far from perfect ($g^{(2)}_{\Gamma}(0)\approx 0.42$ and $\min_{\tau}g^{(2)}_{\Gamma}(\tau)\approx 0.37$). The fact that the minimum antibunching is not at zero is another manifestation of frequency filtering, thinly visible on the figure as small oscillations, but not reproduced at this level of signal by the Monte Carlo data. Correspondingly, the $p(n)$ shows increasingly suppressed probabilities to get higher number of photons.

4.2 Cross-correlations

In this final part, we consider cross-correlations, for which the Mollow triplet is also a particularly suitable lineshape. That is to say, we consider two detectors acquiring data simultaneously. The master equation for two detectors upgrade Eq. (5) to:

[TABLE]

with $\xi_{1}$, $\xi_{2}$ the two detectors.The factors $\alpha_{1}=(1-\chi_{0}-\chi_{1})(1-\chi_{2})$ and $\alpha_{2}=\chi_{0}(1-\chi_{3})$, satisfying simultaneously $0\leq\chi_{0},\chi_{1},\chi_{2},\chi_{3}\leq 1$ and $\chi_{0}+\chi_{1}\leq 1$, take into account the several decay channels of the source: a fraction $\chi_{1}$ into free space, a fraction $\chi_{0}$ to the detector $\xi_{1}$ and the remaining fraction $(1-\chi_{0}-\chi_{1})$ to the detector $\xi_{2}$. In analogy with the case of a single detector, the system described by the master equation (17) has five collapse operators:

[TABLE]

and its associated non-hermitian Hamiltonian becomes

[TABLE]

As for the case of autocorrelations, one could similarly demonstrate the equivalence between cross-correlations to any orders as computed through the frequency-resolved photon correlations and those obtained through Eq. (17) above. Also as was done before for single frequency windows, by applying the Monte Carlo techniques to the detectors, one can thus obtain simulated photon emissions, this time in two frequency windows. Computing the correlations from this raw data provides a numerical version of the theoretical correlations. This is shown in Fig. 7 for the joint emission of the two sidebands on the one hand, and then of the two leapfrog windows on the other hand, both when driving the two-level system at resonance or with a detuning.

While we considered a small subset only of the possible autocorrelations in Fig. 6 for the Monte Carlo data, we could still provide a comprehensive theoretical result at least for $g^{(2}_{\Gamma}(0)$ through the color-coded spectrum. For cross-correlations, however, this would require a 2D plot to reproduce the entire two-photon correlation spectrum. [12] Instead, we consider here the case where one detector is fixed and the other one sweeps the rest of the spectrum, providing the cross-correlations. We then place the other detector for Monte Carlo sampling at the location of interest. As before, we show raw photon emission, but with no time-rescaling as the respective frequencies chosen have similar intensities. We also compare two-photon correlations computed from this data (black and blue lines, at resonance and with detuning respectively) with the theoretical result (red lines). The main difference between these cases and the previous ones is that the two streams of photons are now correlated as the detectors are measuring simultaneously. If restricting to one stream only, we recover the previous cases, so in Fig. 7, panels $i$-$i$ and $ii$-$ii$ on the one hand, as well as $iii$-$iii$ and $iv$-$iv$ on the other hand, can be found in Fig. 6, in panels $iii$ and $v$ respectively. The two cases have different parameters, since with two detectors, the simulation is more computer intensive and we chose a case that provides more emission between the peaks. One can check however how the qualitative shapes of the correlations remain the same. At resonance, both streams provide the same type of correlations, but with detuning, they could be different. This is indeed the case for the leapfrog emission, and small departures can be observed between $v$-$v$ and $vi$-$vi$, both in the theoretical line and the Monte-Carlo generated data: the oscillations are more marked for the low-energy window and a depletion is indeed visible in $vi$-$vi$ that disappeared in $v$-$v$. The side peaks, however, feature similar correlations. This shows again the typical richer dynamics away from the peaks.

Both at resonance or with detuning, what is of interest when detecting in different windows simultaneously is their cross-correlations, as shown in the central column of Fig. 7 with panels $i$-$ii$, $vi$-$v$, $iii$-$iv$ and $viii$-$vii$. There are now clear features in these cross-correlations, whereas the same procedure applied to the streams of the previous cases, features no correlations, i.e., one obtains flat lines. At resonance, the cross-correlations are symmetrical in time. The side peaks correlations feature tiny oscillations which are however too small to be observed with the amount of signal we acquired in the numerical experiment, and they are hidden by its fluctuations. With detuning, they can and do become time-asymmetrical, as shown in panels $vi$-$v$ and $viii$-$vii$. In such cases, the order of detection matters, and in both cases, the detection on first detector, $vi$ or $viii$, respectively, makes it more likely to later detect (with around $1/\gamma_{\sigma}$ delay) a photon on the second detector, $v$ and $vii$, respectively. The strength of such correlations, less than 3, is still fairly modest to call this heralding, but this is the basis for such a mechanism to be exploited with proper engineering, such as Purcell-enhancement. Like before, our procedure yields correlated streams of photons of different frequencies, that we have just shown through their agreement with the theory of frequency-resolved photon correlations, simulate the actual photon emission from the system. One could use this raw data to compute numerically, e.g., counting or time-delay distributions, otherwise not easily accessible theoretically. Of course, the scheme could in principle be extended to any number of detectors and allow consequently higher orders of correlation to be computed in this way. In the limit of an infinite number of detectors, each with a given frequency and vanishing spectral width, one would thus simulate the ideal emission of the system. With a finite number of detectors with a finite bandwidth, one would simulate its filtered emission. We believe that a complexity analysis of the correlations would allow to use the emission of the two-level system as a simpler platform than boson sampling [57] to test quantum supremacy by making a laboratory measurement which no classical computer would be able to simulate.

5 Conclusions and perspectives

In conclusion, we have presented computer experiments that simulate numerically the photon emission from a quantum emitter, specifically, a two-level system under both coherent and incoherent driving, at both low and large pumping. Our approach is based on the Quantum Monte Carlo technique [5, 6, 7, 8] applied to the cascaded formalism. [36, 37] We have shown how the correlations computed from the raw data of the simulation match with the theoretical results provided by the theory of frequency-resolved photon correlations. [11] In the simplest case, we have shown how filtering spoils antibunching and turns a two-level system into a thermal source, albeit in a more subtle way than is usually assumed and for which we have provided exact closed-form expressions. We have also shown more generally how frequency filtering provides a resource to tailor and engineer photon statistics, in particular thanks to its selection of strongly correlated processes such as “leapfrog” transitions that consist in the simultaneous emission of a photon bundle between two non-contiguous dressed states in the level structure of the system. [12] This displays rich and potentially useful features that are captured in the Monte Carlo simulation and that would be similarly observed experimentally. An apparent shortcoming is that the signal is scarce in frequency windows that are the most strongly correlated. This is however a direct consequence of dealing with the quantum part of the signal: there is less of it. Frequency filtering acts as a process akin to distillation, with the same consequence of providing quality at the expense of quantity. [29] Nevertheless, quantum engineering can come at the rescue and already the oldest trick of cavity QED—Purcell enhancement—allows, in some regime, to have all the light of the system emitted as strongly correlated photons. [10] Using cavities to Purcell-enhance leapfrog processes, one can devise new generations of heralded $N$-photon sources, or, even more generally, bring the system to emit in any desired distribution of photons. [48] Such configurations remain to be studied in detail and, of course, implemented in the laboratory. This should provide one route for universal multi-photon sources, with heralded $N$ photon sources as the most elementary realization. Since leapfrog processes are energy-conserving $N$-photon relaxations, they also appear particularly suitable for energy-time entanglement emission, that power a class of quantum-cryptographic protocols with technical advantages as compared to those based on entangling in polarization. The latest work from Peiris et al., [58] who is so far leading the laboratory implementation of this emerging branch of quantum optics, focused on the side peaks emission and, as a consequence, failed to break the barrier of a Bell violation. This has been argued no to be a proof of nonlocality anyway [59] due to its 50% post-selection. [60] It is easily computed that leapfrog emission would break the Franson limit, but in the light of the Franson configuration’s loophole, [61] the new challenge is to turn to stricter conditions of Bell violations such as Chained Bell’s inequalities. While this has been recently demonstrated, [62] the tunable statistics from the Mollow triplet and its windows of strong correlations make it a promising platform to further test and advance this line of research. Finally, the combinatoric aspects that quickly make such simple problems numerically forbidding also suggest that a two-level system could be used in the laboratory for tests of quantum supremacy directly from photon detections, without a complex system of beam splitters intervening to bring in the quantum complexity. [63, 64] All these results leave much for room for further works, and we foresee that frequency-resolved photon correlations will become a major theme of photonics. They are relevant even when they are ignored and awareness of the underlying physics should allow to considerably optimize, tune and expand the range of applications of quantum light sources.

Acknowledgments

Funding by the Newton fellowship of the Royal Society, the POLAFLOW ERC project No. 308136, the Ministry of Science and Education of Russian Federation (RFMEFI61617X0085), the Spanish MINECO under contract FIS2015-64951-R (CLAQUE) and by the Universidad Autónoma de Madrid under contract FPI-UAM 2016 is gratefully acknowledged.

Author contributions statement

EdV an JCLC developed the theory, JCLC implemented and performed the numerical simulations, FPL and JCLC composed the figures, FPL wrote the text and supervised the research. All the authors conceived the work, discussed the results and contributed to the manuscript.

Competing Financial Interests statement

The authors declare no competing financial interests.