Frequencyresolved Monte Carlo
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 frequencyresolved photon correlations, extending the range of applications based on correlations of photons of prescribed energy, in particular those of a photoncounting character. We apply the technique to autocorrelations of photon streams from a twolevel 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 twophoton emission events.
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 solid angle surrounding the emitter, this allows to perform a computer experiment of photodetections. From such “clicks” (as we will call a detected photon), one can for instance compute the Glauber correlation functions that measure the deviations of intensity correlations at the 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 after one detection, with no other photon in between ( 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 Carlogenerated photons have also been used to support the introduction of the photon “bundle”, [10] to distinguish the case of photon sources from stronglycorrelated emission at the photon level.
In this text, we apply the quantumjump 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 frequencyresolved photoncorrelations. [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 colorresolved photons. While the quantum Monte Carlo method has been used to compute the power spectrum as well as timeseries of photon emission, [14, 15, 16, 17, 18, 19, 20] its combined use for both time and energyresolved 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 weakcoupling 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 effectivethree 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 selfconsistently 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 twolevel 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 singlephoton emitter, but of photons with wildly fluctuating frequencies (so not indinstinguishable). We illustrate how the loss of antibunching in time deviates notably from the singleexponential 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 frequencyresolved 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 normalordering and timeintegrals in spaces of many dimensions. The technique relies on coupling the system to sensors with strength 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 , without affecting the system’s dynamics. The computation of normalized correlators results in quantities that are independent to first order and exact in the limit (taking, in general, ). 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 socalled “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 twolevel 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 twolevel system while the detector must be truncated high enough to provide a closeenough approximation, which depends on the dynamics of the system, and is therefore a tricky question. Furthermore, when considering crosscorrelations, instead of twolevel systems, one is dealing with harmonic oscillators and the problem becomes numerically forbidding, while the sensors’ Hilbert space scales as which is still tractable for much larger than anything that has been considered so far experimentally. For autocorrelations of the th order with the sensor method, one can also use an harmonic oscillator truncated to excitations instead of twolevel 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 the annihilation operator of a source and that of a sensor probing it, we can describe their joint dynamics by a Liouvillian equation
(1) 
where is the Hamiltonian that describes the internal dynamics of the source, is the decay rate of the source, is the decay rate of the sensor and . The dynamics of an arbitrary operator under the action of this Liouvillian is described with the notations of Ref. [11], by the equation (we set along the paper)
(2) 
where is the unit matrix, are normalordering superoperators for the operators and is a vector of correlators for the th and th powers of the sensor operators , and spanning in normal order all powers of the , operators, i.e., . The notation means that all other terms are of higher order than . The matrix provides the dynamics for the source, , and is independent of the sensor at the lowest order in . At this stage, we do not assume any property of or , which could be bosonic (in which case , , and are unbounded) or fermionic (in which case , , and are 0 or 1). Equation (2) can be integrated, which yields
(3) 
where . This in turn can be solved recursively, down to where the equation selftruncates. This provides terms, all with the same leading order of . 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 with the th power of components of , themselves of order , so that such ratio cancel out to leading order. While the higherorder terms do not have to cancel exactly, as is the case of the leadingorder 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 independent term, in contrast to unnormalised correlators that vanish with . Therefore, in the limit , 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 backaction from the detector to the source. The master equation describing this asymmetric coupling reads
(4) 
The last threeterms of Eq. (4) can be rewritten in the Lindblad form as
(5) 
where is the joined decay operator of the whole system, source and detector, and the interpretation of the factor 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 , for , 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 has the same equation as in the sensor formalism, by computing explicitly the equation for in the cascaded formalism, Eq. (5). This reads, to all orders in the coupling in this case:
(6) 
Remarkably, this equation has the same form as Eq. (2) with . Even though is complex and a vanishing quantity in Eq. (2), with higher order corrections, and 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 , the phase of the coupling , 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 that brings cross terms of and . The sensor formalism, however, has no such joint decay emission and the phase of does not play any role, so that 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 photodetection 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 frequencyresolved 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 can be computed exactly (by recurrence) in terms of lower order ones only, with and . This means that both methods can be applied using twolevel systems as detectors at different frequencies, in order to compute cross correlations, or with a single harmonic oscillator truncated at excitations, to compute the th order monocromatic autocorrelation function . 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 detectors as harmonic oscillators with a high enough truncation to provide converged results. The simulation is conveniently implemented through the quantumjump approach. The dynamics of the system is thus described by a wavefunction 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 , where the evolution of the wavefunction is governed by two elements: a nonHermitian Hamiltonian and random quantum jumps. In a system described by the master equation , the nonHermitian Hamiltonian is constructed as , and the evolution of the wavefunction is given either by
depending on whether the system undertook or not a quantum jump, respectively. The probability that this happened in the interval due to the operator is proportional to the mean value of this operator, namely . For the system described by the master equation (4), the collapse operators are
(7) 
whereas the nonHermitian Hamiltonian is given by
(8a)  
(8b) 
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. (45), in the frequency window determined by the detector, match indeed with the correlations predicted by the theory of frequencyresolved photon correlations. [11] Here, we apply this technique to the driven twolevel 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 crosscorrelations. 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 twolevel system
The simplest system in quantum physics is the twolevel system. We find it inspiring that it is still a Research problem, constantly offering new questions. Manifestly, in a quantum universe, the twolevel 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 , with Hamiltonian
(9) 
for a twolevel system with free energy , and Liouvillian , with the inverse lifetime, one finds a simple enough dynamics of Glauber’s second order correlator:
(10) 
with, in particular, , 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 counterintuitive “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 stronglycorrelated character of the twolevel system emission becomes clear and compelling when computing intensity correlations from the clicks, defined as the density of probability of finding two photons with a time difference . Specifically, from the times of detection , we compute for all with the total number of detected photons (here ) 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 in Fig. 1(e), left. One sees an overall plateau, indicating that photons have the same distribution for longtime separations as if they were emitted by a Poisson process (randomly). But one also observes a clear dip at , 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, , with coincidences, i.e., simultaneous detection of two photons, less likely to occur than other closely spaced detections, with perfect suppression of coincidences when . Since these correlations wash out at long times, one has . The time it takes to reach this plateau is the secondorder 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 , the Monte Carlo data provides a smooth curve in the window of strong correlations. In our simulation, the was and the binning size was taken twice as large, corresponding to the two closelyspaced vertical lines on the right panel of Fig. 1(e), bounding 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 twolevel system
3.1 Incoherent excitation
The effect of a Lorentzian filter on the statistics of emission of an incoherently excited twolevel 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 MonteCarlo simulations of roughly 10 000 clicks each (25 000 for the narrowest filter in case ). 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 () to thermal ones () passing by auxiliary distributions. In the former case, one observes the characteristic antibunching, equallyspaced like distribution of a twolevel 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 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 , as is shown in the eight panels at the bottom of Fig. 2. The case of closetouncorrelated 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 –. There, one should increase the binning and consider larger time windows, as shown in green for case that assumes a binning of instead of for the other cases, and plot the correlations in a time window 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 strongdriving, [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 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 photondetection events, where the spacing appears more regular and between clumps of photons. As far as continuous streams are concerned, this suggests that such stronglyoscillating 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 narrowenough 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 twolevel system unchanged. [41] As such, the spontaneously emitted photons react to filtering in a similar way than the incoherently pumped twolevel system, hence the observed bunching for narrowing filters linewidths. The similarity is only partial, however, as instead of thermalization, with , the transition is to a superchaotic state, with in the limit of infinite pumping [12] (for the parameters considered here, we find ). More strikingly, when filtering well within the central peak, one then isolates the Rayleigh () peak again and reverts to the lowpumping 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 of a frequencyresolved emitter, i.e., the probabilities for it to have photons, which are, equivalently, the diagonal elements of its effective density matrix. Namely, the frequencyfiltered 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 rather than the original . 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 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:
(11) 
and as, on the other hand, there is conservation of energy (assuming an ideal detector):
(12) 
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:
(13) 
The statistics of the photon emission from the effective source can now be obtained by inverting the relation (with large enough) to provide the probabilities for the effective source to have quanta of excitation, for integer . In the cascaded formalism, the correlators of the detectors, 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 twolevel 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 reads^{1}^{1}1Note that there is a typo in Ref. [32] with an extra factor 4; the correct result is as given here.
(14) 
and the coherent excitation reads
(15) 
with . From a Monte Carlo simulation in a time , the detector population is obtained as the ratio between the total number of clicks recorded and . 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 photodetection 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 in this way. This allows us to access new classes of quantum steady states, tailored by frequencyfiltering. We now illustrate how this takes shape in the case already discussed of filtered twolevel system emission, starting with the case of incoherent excitation. This is shown in Figs. 4(a–b) for the filtered saturated twolevel system, i.e., where the system is held in its excited state by very large pumping, , so that the density matrix reads and . 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 ( being the Kronecker function). The application of a filter turns the system from a twolevel emitter to a source able to deliver more than one photon at a time, namely, for narrow enough filtering, for all , with 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 singlephoton source turns it into a black body with nonzero probability to emit 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), . In both cases, one can see in this way at which point filtering prevents a singlephoton source to emit nonclassical states of light, [43] for instance by comparing to , the smallest probability above which a state is nonGaussian. [44]
4 Tuneable statistics from the Mollow triplet
4.1 Autocorrelations
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 twolevel 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 CohenTannoudji 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 twophoton transitions, one can foresee some correlations between the peaks. The transitions for any or , that go down the Mollow ladder, could be expected to result in bunching. In contrast, since one cannot chain in this way and or and , one can expect antibunching. Inspection of all the combinations of twophoton 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 timeordering and/or detuning.) Early calculations by Apanasevich and Kilin [46] and Cohen Tannoudji and Reynaud [45] confirmed the sidepeaks antibunching in autocorrelation and their mutual bunching in crosscorrelations. 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 coexistence of bunched and antibunched emission events. This shows that while intuition is strongly supported by the dressedatom 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 sidepeaks emission).
With the theory of frequencyresolved photon correlations, [11] it is straightforward to compute such correlations exactly, like Mollow did for the luminescence spectrum, without referring to the dressedstates structure. This also allows to consider cases ouside their frequency windows, in fact, the complete landscape of twophoton correlations can be obtained. [12, 29] In the case of the Mollow triplet, [12, 48] it shows that the triplet structure reverberates at the twophoton level, through the apparition of a set of 3 hyperplanes, that obey the “leapfrog” equations
(16) 
(The same applies at the photon level [48]). The triplet structure comes, at any photonnumber level, from the three possibilities to join the two dressed initial and final states. Note that while there are four transitions, is assumed sufficiently large for two transitions to be degenerate. The name of “leapfrog” comes from the fact that, at the photon level, transitions can occur by jumping over intermediate manifolds. This relaxes energy conservation for the individual photons and restricts it to the combined emission. In this sense, this is a ( photons) version of Göpper Mayer’s twophoton processes. [49] that is however typically difficult to access (a notable exception is the planetary nebulae continuum [50]). The case 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 case (that is, overlooking the 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 stronglycorrelated 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 twophoton spectrum remarkably featureless, in contrast to our twophoton spectrum that is rich from photon correlations flourishing away from the peaks. Their nonnormalized spectrum is correct but, we believe, is also not interesting as it merely shows that first order processes smother secondorder ones, as is however wellknown and expected. We show, in contrast, that the scarce signal from higherorder processes has stronger correlations than those from firstorder 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 firstorder processes, [10, 55] which would not happen would the correlation be an artifact due to normalization. We will show below thanks to the frequencyresolved 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 dressedatom 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 twolevel system in a sequence of coherent absorption and emission. We have in fact shown [33] how even fluorescence in the lowdriving regime does not consist of Rayleigh scattering events but form an intricate interference between emission and absorption, that powers the singlephoton emission mechanism by suppressing the fat tails of the Lorentzian and turning the lineshape into a 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 dressedatom 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 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 transition with one photon nondegenerate with the two others (of interest for photonheralding of twophoton emission), is particularly weak. This is because it is resonant with the transition, that breaks this channel of relaxation by interposing a real state in the threephoton emission. Configurations 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 dimensional space to harness the best sought configuration of multiphoton 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 photodetection 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 singledetector 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 detectors to obtain the same result but with correct crosscorrelations (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:

Photons from the central peak.

Case where (usually attributed to thermal light).

Photons from leapfrog emission.

Case where (usually attributed to coherent or uncorrelated light).

Photons from a side peak.
The central peak is partially thermalized, with a that closely resembles the form of thermal fluctuations, . 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 twolevel system is considerably more complicated than its incoherent counterpart and we could not find, so far, a general closedform expression for in this case for arbitrary frequencies. Applying the technique of effectivequantum state reconstruction from the correlators, described in Section 3.3, we find that the statistics fits well with a cothermal distribution with of thermal emission and 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, , which features , 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)
The most interesting window, , 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 in a given manifold to the state 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 photodetection point of view, they should arise as more occurrences of closelyspaced twophoton clicks than if the emission was uncorrelated. In particular, their rate of coincidences should increase, leading to , or socalled 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 highquality signal, revealing the fine details of its structure. Note as well that on the realtime series of clicks, out of the nine photons emitted, four came as twophoton bundles (the fifth and sixth clicks are so closely spaced as almost overlapping; other ticks are singlephoton 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 twophoton emission events that occurs, from the dressedatom 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 photon emission, and for which filtering has been shown to considerably boost the purity of the quantum emission. [56, 55] Also further photoncounting characterization would certainly be enlightening, and preliminary investigations show that the percentage of closelyspaced photons is over an order of magnitude higher in than in the others at the exception of , as compared to which it is only about 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 as the most dissimilar one as compared to the others, featuring a neat kink at the probability to have two photons, , showing the relative predominance of twophoton emission. Overall, this simulation makes it obvious that the emission in this frequency window suffers from no artifact of postselection or normalization, but does indeed provide strongly correlated photon streams.
The fourth frequency window, , chosen for its of uncorrelated emission, is also a case that shows strong departures at nonzero 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, , we come back to a case well studied in the literature, of antibunched emission, albeit far from perfect ( and ). 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 shows increasingly suppressed probabilities to get higher number of photons.
4.2 Crosscorrelations
In this final part, we consider crosscorrelations, 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:
(17) 
with , the two detectors.The factors and , satisfying simultaneously and , take into account the several decay channels of the source: a fraction into free space, a fraction to the detector and the remaining fraction to the detector . In analogy with the case of a single detector, the system described by the master equation (17) has five collapse operators:
and its associated nonhermitian Hamiltonian becomes
(19) 
As for the case of autocorrelations, one could similarly demonstrate the equivalence between crosscorrelations to any orders as computed through the frequencyresolved 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 twolevel 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 through the colorcoded spectrum. For crosscorrelations, however, this would require a 2D plot to reproduce the entire twophoton 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 crosscorrelations. 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 timerescaling as the respective frequencies chosen have similar intensities. We also compare twophoton 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  and  on the one hand, as well as  and  on the other hand, can be found in Fig. 6, in panels and 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  and , both in the theoretical line and the MonteCarlo generated data: the oscillations are more marked for the lowenergy window and a depletion is indeed visible in  that disappeared in . 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 crosscorrelations, as shown in the central column of Fig. 7 with panels , ,  and . There are now clear features in these crosscorrelations, whereas the same procedure applied to the streams of the previous cases, features no correlations, i.e., one obtains flat lines. At resonance, the crosscorrelations 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 timeasymmetrical, as shown in panels  and . In such cases, the order of detection matters, and in both cases, the detection on first detector, or , respectively, makes it more likely to later detect (with around delay) a photon on the second detector, and , 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 Purcellenhancement. Like before, our procedure yields correlated streams of photons of different frequencies, that we have just shown through their agreement with the theory of frequencyresolved photon correlations, simulate the actual photon emission from the system. One could use this raw data to compute numerically, e.g., counting or timedelay 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 twolevel 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 twolevel 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 frequencyresolved photon correlations. [11] In the simplest case, we have shown how filtering spoils antibunching and turns a twolevel system into a thermal source, albeit in a more subtle way than is usually assumed and for which we have provided exact closedform 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 noncontiguous 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 Purcellenhance leapfrog processes, one can devise new generations of heralded 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 multiphoton sources, with heralded photon sources as the most elementary realization. Since leapfrog processes are energyconserving photon relaxations, they also appear particularly suitable for energytime entanglement emission, that power a class of quantumcryptographic 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% postselection. [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 twolevel 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 frequencyresolved 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 FIS201564951R (CLAQUE) and by the Universidad Autónoma de Madrid under contract FPIUAM 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.
References
 1. Ulam, S. M. Adventures of a mathematician (University of California Press, 1991).
 2. Knuth, D. E. 3:16 Bible texts illuminated (AR Editions, 1991).
 3. Ceperley, D. & Alder, B. Quantum Monte Carlo. Science 231, 555 (1986).
 4. Plenio, M. B. & Knight, P. L. The quantumjump approach to dissipative dynamics in quantum optics. Rev. Mod. Phys. 70, 101 (1998).
 5. Zoller, P., Marte, M. & Walls, D. F. Quantum jumps in atomic systems. Phys. Rev. A 35, 198 (1987).
 6. Carmichael, H. J., Singh, S., Vyas, R. & Rice, P. R. Photoelectron waiting times and atomic state reduction in resonance fluorescence. Phys. Rev. A 39, 1200 (1989).
 7. Dalibard, J., Castin, Y. & Mølmer, K. Wavefunction approach to dissipative processes in quantum optics. Phys. Rev. Lett. 68, 580 (1992).
 8. Mølmer, K. & Castin, Y. Monte Carlo wavefunctions in quantum optics. Quantum Semiclass. Opt. 8, 49 (1996).
 9. Hegerfeldt, G. C. The quantum jump approach and some of its applications. Lect. Notes Phys. 789, 127 (2009).
 10. Sánchez Muñoz, C. et al. Emitters of photon bundles. Nat. Photon. 8, 550 (2014).
 11. del Valle, E., GonzálezTudela, A., Laussy, F. P., Tejedor, C. & Hartmann, M. J. Theory of frequencyfiltered and timeresolved photon correlations. Phys. Rev. Lett. 109, 183601 (2012).
 12. GonzálezTudela, A., Laussy, F. P., Tejedor, C., Hartmann, M. J. & del Valle, E. Twophoton spectra of quantum emitters. New J. Phys. 15, 033036 (2013).
 13. Peiris, M. et al. Twocolor photon correlations of the light scattered by a quantum dot. Phys. Rev. B 91, 195125 (2015).
 14. Dum, R., Parkins, A. S., Zoller, P. & Gardiner, C. W. Monte Carlo simulation of master equations in quantum optics for vacuum, thermal, and squeezed reservoirs. Phys. Rev. A 46, 4382 (1992).
 15. Marte, P., Dum, R., Taïeb, R., Lett, P. D. & Zoller, P. Quantum wave function simulation of the resonance fluorescence spectrum from onedimensional optical molasses. Phys. Rev. Lett. 71, 1335 (1993).
 16. Mølmer, K., Castin, Y. & Dalibard, J. Monte Carlo wavefunction method in quantum optics. J. Opt. Soc. Am. B 10, 524 (1993).
 17. Garraway, B., Kim, M. & Knight, P. Quantum jumps, atomic shelving and Monte Carlo fluorescence spectra. Opt. Commun. 117, 560 (1995).
 18. Hegerfeldt, G. C. & Plenio, M. B. Conditional resonancefluorescence spectra of single atoms. Phys. Rev. A 53, 1164 (1996).
 19. Plenio, M. B. Narrow absorption lines induced by electron shelving. J. Mod. Opt. 43, 753 (1996).
 20. Schack, R., Brun, T. A. & Percival, I. C. Quantumstate diffusion with a moving basis: Computing quantumoptical spectra. Phys. Rev. A 53, 2694 (1996).
 21. Tian, L. & Carmichael, H. J. Quantum trajectory simulations of twostate behavior in an optical cavity containing one atom. Phys. Rev. A 46, R6801(R) (1992).
 22. Sánchez Muñoz, C., del Valle, E., Tejedor, C. & Laussy, F. Violation of classical inequalities by photon frequency filtering. Phys. Rev. A 90, 052111 (2014).
 23. del Valle, E., Silva, B., Muñoz, C. S., Carreño, J. C. L. & Laussy, F. Loss of antibunching. arXiv:xxxxx (2017).
 24. López Carreño, J. C., Sánchez Muñoz, C., Sanvitto, D., del Valle, E. & Laussy, F. P. Exciting polaritons with quantum light. Phys. Rev. Lett. 115, 196402 (2015).
 25. Mukamel, S. & Dorfman, K. E. Nonlinear fluctuations and dissipation in matter revealed by quantum light. Phys. Rev. A 91, 053844 (2015).
 26. Kazimierczuk, T. et al. Photonstatistics excitation spectroscopy of a quantumdot micropillar laser. Phys. Rev. Lett. 115, 027401 (2015).
 27. Dory, C. et al. Tuning the photon statistics of a strongly coupled nanophotonic system. Phys. Rev. A 95, 023804 (2017).
 28. Shatokhin, V. N. & Kilin, S. Y. Correlation functions in resonance fluorescence with spectral resolution: Signalprocessing approach. Phys. Rev. A 94, 033835 (2016).
 29. del Valle, E. Distilling one, two and entangled pairs of photons from a quantum dot with cavity QED effects and spectral filtering. New J. Phys. 15, 025019 (2013).
 30. GonzálezTudela, A., del Valle, E. & Laussy, F. P. Optimization of photon correlations by frequency filtering. Phys. Rev. A 91, 043807 (2015).
 31. Sánchez Muñoz, C., Laussy, F. P., Tejedor, C. & del Valle, E. Enhanced twophoton emission from a dressed biexciton. New J. Phys. 17, 123021 (2015).
 32. López Carreño, J. C. & Laussy, F. P. Excitation with quantum light. I. Exciting a harmonic oscillator. Phys. Rev. A 94, 063825 (2016).
 33. López Carreño, J. C., Sánchez Muñoz, C., del Valle, E. & Laussy, F. P. Excitation with quantum light. II. Exciting a twolevel system. Phys. Rev. A 94, 063826 (2016).
 34. Vogel, W. & Welsch, D.G. Quantum Optics (WileyVCH, 3, 2006).
 35. Flayac, H. & Savona, V. Heralded preparation and readout of entangled phonons in a photonic crystal cavity. Phys. Rev. Lett. 113, 143603 (2014).
 36. Gardiner, C. W. Driving a quantum system with the output field from another driven quantum system. Phys. Rev. Lett. 70, 2269 (1993).
 37. Carmichael, H. An open systems approach to Quantum Optics, chap. 6 Photoelectric Detection II, 110 (Springer, 1993).
 38. Carmichael, H. J. Statistical methods in quantum optics 1 (Springer, 2008).
 39. Feller, W. An introduction to probability theory and its applications (John Wiley & Sons, 1968).
 40. Mollow, B. R. Power spectrum of light scattered by twolevel systems. Phys. Rev. 188, 1969 (1969).
 41. Reynaud, S. La fluorescence de résonance: étude par la méthode de l’atome habillé. Annales de Physique 8, 315 (1983).
 42. Zubizarreta Casalengua, E., López Carreño, J. C., del Valle, E. & Laussy, F. P. Structure of the harmonic oscillator in the space of particle glauber correlators. J. Math. Phys. 58, 062109 (2017).
 43. Carreño, J. C. L., Casalengua, E. Z., del Valle, E. & Laussy, F. P. Criterion for single photon sources. arXiv:1610.06126 (2016).
 44. Filip, R. & L. Mišta, J. Detecting quantum states with a positive Wigner function beyond mixtures of Gaussian states. Phys. Rev. Lett. 106, 200401 (2011).
 45. CohenTannoudji, C. N. & Reynaud, S. Dressedatom description of resonance fluorescence and absorption spectra of a multilevel atom in an intense laser beam. J. Phys. B.: At. Mol. Phys. 10, 345 (1977).
 46. Apanasevich, P. A. & Kilin, S. Y. Photon bunching and antibunching in resonance fluorescence. J. Phys. B.: At. Mol. Phys. 12, L83 (1979).
 47. Schrama, C. A., Nienhuis, G., Dijkerman, H. A., Steijsiger, C. & Heideman, H. G. M. Destructive interference between opposite time orders of photon emission. Phys. Rev. Lett. 67 (1991).
 48. López Carreño, J. C., del Valle, E. & Laussy, F. P. Photon correlations from the Mollow Triplet. Laser Photon. Rev. 11, 1700090 (2017).
 49. GöppertMayer, M. Über Elementarakte mit zwei Quantensprüngen. Annalen der Physik 401, 273 (1931).
 50. Spitzer Jr., L. & Greenstein, J. L. Continuous emission from planetary nebulae. Ap. J. 114, 407 (1951).
 51. Apanasevich, P. & Kilin, S. Y. Lightinduced correlations in spontaneous emission. Phys. Lett. A 62, 83 (1977).
 52. Shatokhin, V. & Kilin, S. Secondorder coherence and state reduction in resonance fluorescence with spectral resolution. Opt. Commun. 174, 157 (2000).
 53. Shatokhin, V. N. & Kilin, S. Y. Correlation measurements in resonance fluorescence with spectral resolution and atomic inversion via detection of a spectrally filtered photon. Phys. Rev. A 63, 023803 (2001).
 54. Shatokhin, V. N. & Kilin, S. Y. Atomic state reduction and energy balance relation in spectrally resolved resonance fluorescence. Opt. Commun. 210, 291 (2002).
 55. Sánchez Muñoz, C., Laussy, F. P., del Valle, E., Tejedor, C. & GonzálezTudela, A. Filtering multiphoton emission from stateoftheart cavity qed. arXiv:1707.03690 (2017).
 56. Sánchez Muñoz, C. Generation of Nonclassical States of Light. Ph.D. thesis, Universidad Autónoma de Madrid (2016).
 57. Aaronson, S. & Arkhipov, A. The computational complexity of linear optics. Proceedings of the 43rd Annual ACM Symposium on Theory of Computing 333 (2011).
 58. Peiris, M., Konthasinghe, K. & Muller, A. Franson interference generated by a twolevel system. Phys. Rev. Lett. 118, 030501 (2017).
 59. Jogenfors, J., Cabello, A. & Larsson, J.Å. Comment on “Franson interference generated by a twolevel system”. arXiv:1703.05055 (2017).
 60. Aerts, S., Kwiat, P., Larsson, J.Å. & Żukowski, M. Twophoton Fransontype experiments and local realism. Phys. Rev. Lett. 83, 2872 (1999).
 61. Jogenfors, J., Elhassan, A. M., Ahrens, J., Bourennane, M. & Larsson, J.Å. Hacking the Bell test using classical light in energytime entanglement–based quantum key distribution. Science Advances 1, e1500793 (2015).
 62. Tomasin, M. et al. Highvisibility timebin entanglement for testing chained Bell inequalities. Phys. Rev. A 95, 032107 (2017).
 63. Broome, M. et al. Photonic boson sampling in a tunable circuit. Science 15, 794 (2013).
 64. Crespi, A. et al. Integrated multimode interferometers with arbitrary designs for photonic boson sampling. Nat. Photon. 7, 545 (2013).