Conventional and unconventional photon statistics

We show how the photon statistics emitted by a large variety of light-matter systems under weak coherent driving can be understood, to lowest order in the driving, in the framework of an admixture of (or interference between) a squeezed state and a coherent state, with the resulting state accounting for all bunching and antibunching features. One can further identify two mechanisms that produce resonances for the photon correlations: i) conventional statistics describes cases that involve a particular quantum level or set of levels in the excitation/emission processes with interferences occurring to all orders in the photon numbers, while unconventional statistics describes cases where the driving laser is far from resonance with any level and the interference occurs for a particular number of photons only, yielding stronger correlations but only for a definite number of photons. Such an understanding and classification allows for a comprehensive and transparent description of the photon statistics from a wide range of disparate systems, where optimum conditions for various types of photon correlations can be found and realized.

Quantum optics was born with the study of photon statistics [1]. Following Hanbury Brown's discovery of photon bunching [2,3] and Kimble et al.'s observation of antibunching [4], there has been a burgeoning activity of tracking how pairs of detected photons are related to each other. From the various mechanisms that generate correlated photons, one that turns an uncorrelated stream into antibunched photons has attracted much attention across different platforms for its practicality of operation (with a laser) and appealing underlying mechanism [5][6][7][8][9][10][11][12][13]. This so-called "blockade" effect describes how the occupation of an energy level by a particle forbids another particle to occupy the same level. At such, it is reminiscent of Pauli's exclusion principle [14,15] and arXiv:1901.09030v2 [quant-ph] 24 May 2019 indeed the first type of blockade involved electron in the so-called "Coulomb blockade" [16][17][18]. However, while Pauli's principle relies on the antisymmetry of fermionic wavefunctions, one can also implement a "bosonic blockade" with photons exciting an anharmonic system [19], based on nonlinearities in the energy levels due to selfinteractions. The idea is simple: when photons are resonant with the frequency of the oscillator, a first photon can excite the system, but due to photon-photon interactions, a subsequent photon is now detuned from the oscillator's frequency. If its energy is not sufficient to climb the ladder of states, it cannot excite the system, that remains with one photon. In this way, one can turn a coherent-that is, uncorrelated-stream of photons, with its characteristic Poissonian fluctuations, into a more ordered stream of separated photons, effectively acting as a "photon turnstile". The quality of such a suppression of the Poisson bursts can be measured at the two-photon level with Glauber's second-order correlation function g (2) (τ ), that compares the coincidences in time to those expected from a random process of same intensity. Correlations decrease from 1, with no blockade, towards 0 as the nonlinearity U increases. In the limit where U becomes infinite, putting the second excited state arbitrarily far and realizing a two-level system (2LS), a second photon is strictly forbidden and g (2) becomes perfectly antibunched. For open bosonic systems, the ratio of the interaction to the decay rate is an important variable for the blocking to be effective. The "blockade" regime is reached when interactions start to dominate the dissipation. This can be marked as the onset of antibunching: one photon starts to suppress the next one [20].
The driven damped anharmonic system is an important model, not least because it is one of the few cases to enjoy an exact analytical solution [21]. While much of the mechanism is contained in this particular case, compound systems-where the anharmonic system is coupled to a single-mode cavity-have also attracted considerable attention. This describes for instance interacting quantumwell excitons coupled to a microcavity mode [22]. The effect is then known as "polariton-blockade", after the eponymous light-matter particles that constitute the elementary excitations of such systems. This configuration was first addressed theoretically by Verger et al. [23] who studied the response of the cavity around the lowerpolariton resonance, predicting antibunching indeed, although of too small magnitude with the parameters of typical systems to be observed easily. This spurred interest in polariton boxes and other ways of confining polaritons to enhance their interactions [24]. A few years later, Liew and Savona computed a much stronger antibunching from a seemingly distinct compound system with the same order of nonlinearity, namely, coupled cavities [25]. This so-called "unconventional polariton blockade" was quickly understood as originating not from the particular configuration of coupled cavities with weak Kerr nonlinearities but from a subtler type of blocking, due to destructive interferences between probability amplitudes whenever there are two paths that can reach the excited state with two photons [26]. This result has generated considerable attention, although it was later remarked [27] that it was a known effect [28,29], observed decades earlier [30] where it received a much smaller followup. Besides, Lemonde et al. [27] further clarified how unconventional photon blockade is connected to squeezing rather than single-photon states, which had been presented as one of the main interest of the effect. Recently, both conventional [31,32] and unconventional [33,34] blockades have been reported in solid-state systems, where the 2010 revival of the idea had triggered intense activity.
In this text, we provide a unifying picture of the two types of polariton blockades. We show how they typically sit next to each other in interacting coupled light-matter systems along with other phenomenologies, that produce superbunching instead of the blockade antibunching. In particular, we show that they are both rooted in the single-component system, either a 2LS or an anharmonic oscillator, with strong photon correlations produced by interfering the emitter's incoherent signal with a coherent fraction. We will nevertheless highlight how the two blockades are intrinsically different mechanisms with different characteristics. Most importantly, the conventional blockade, based on dressed-state blocking, yields photon antibunching at all orders in the number of photons, i.e., g (N ) → 0 for all N ≥ 2, while the unconventional blockade can only target one N in isolation, producing bunching for the others. Another apparent similarity is that both types of blockades produce the same state in what concerns the population and the twophoton correlation g (2) at the lowest order in the driving, but differences occur at higher orders, namely, at the second-order for g (2) and at the third-order for the population. Differences exist already at the lowest order in the driving for g (3) and higher-order photon correlations, making it clear that the two mechanisms differ substantially and produce different states, despite strong resemblances in the quantities of easiest experimental reach. The state produced by both types of blockade at the twophoton level results from a simple interference between a squeezed state and a coherent state. While the squeezing is typically produced by the emitter, the coherent fraction can be either brought from outside, idoneously, as a fraction of the driving laser itself-a technique known as "homodyning"-or be produced internally by the driven system itself, a concept introduced in the literature under the apt qualification of "self-homodyning" [28]. We will thus highlight that, to this order, essentially the same physics-of tailoring two-photon statistics by admixing squeezed and coherent light, discussed in Section IIItakes place in a variety of platforms, overviewed in Section II. We will further synthesize this picture by unifying the cases where the nonlinearity is i) strong, namely, provided by a 2LS (Section IV) or on the contrary ii) weak, namely, provided by an anharmonic oscillator (Sec-tion V), and how these are further generalized in presence of a cavity where self-homodyning becomes a compelling picture since a cavity is an ideal receptacle for coherent states. This brings the 2LS into the Jaynes-Cummings model (Section VI) and the anharmonic oscillator into microcavity polaritons (Section VII), respectively. There are many variations in between all these configurations, that the literature has touched upon in many forms, as we briefly overview in the next Section. For our own discussion, while we have tried to retain as much generality as possible for the variables that play a significant role, we do not include for the sake of brevity all the possible combinations, which could of course be done would the need arise for a given platform. Section VIII summarizes and concludes.

II. A SHORT REVIEW OF BLOCKADES AND RELATED EFFECTS
We will keep this overview very short since a thorough review would require a full work on its own. A good review is found in Ref. [35]. This Section will also allow us to introduce the formalism and notations. For details of the microscopic derivation, we refer to Ref. [23].
A fairly general type of photon blockade is described by the Hamiltonian + Ω a e i a t a + Ω b e i b t b (1c) + h.c.
where ω c is the free energy of the modes c = a, b, both bosonic, g describes their Rabi coupling, giving rise to polaritons as eigenstates of line (1a), U c are the nonlinearities of the respective modes, here again for c = a, b, and Ω c describes resonant excitation at the energy c . This is brought to the dissipative regime through the standard techniques of open quantum systems, namely, with a master equation in the Lindblad form ∂ t ρ = i [ρ, H] + c=a,b (γ c /2)L c ρ, where the superoperator L c ρ ≡ 2cρc † − c † cρ − ρc † c describes a decay rate of mode c at rate γ c (we do not include incoherent excitation nor dephasing, which are other parameters one could easily account for in the following analysis from a mere technical point of view).
Particular cases or variations of Eq. (1) have been studied in a myriad of works, even when restricting to those with a focus on the emitted photon statistics. This ranges from cases retaining one mode only [36] to the most general form of Eq. (1) [37][38][39][40][41][42] with further variations (such as using pulsed excitation [43]). A first consideration of the effect of field admixture on the photon statistics was made by Flayac and Savona [39], who found that that the conditions for strong correlations are shifted rather than hampered. This touches upon, in the framework of input/output theory, the mechanisms of mixing fields that we will highlight in the following, where we will show that beyond being altered, correlations can be drastically optimized (becoming exactly zero to first order for antibunching and infinite for bunching). In a later work [40], they further progressed towards fully exploiting homodyning by including a "dissipative, one-directional coupling" term, which allowed them to achieve a considerable improvement of the photon correlations, especially in time, with suppression of oscillations and the emergence of a plateau at small time delays. This is due to the same mechanism than the one we used with identical consequences but in a different context [44] (a twolevel system admixed to an external laser). Emphasis should also be given to the bulk of work devoted to related ideas of mixing fields by the Vucković group, starting with their use of self-homodyning to study the Mollow triplet in a dynamical setting [45]. Initially used as a suppression technique to access the quantum emitter's dynamics by cancelling out the scattered coherent component from their driving laser [46,47], they later appreciated the widespread application of their effect and its natural occurence in other systems [48], where it had passed unnoticed, as well as the benefits of a tunability of the interfering component [49], which they proposed in the form of a partially transmitting element in an onchip integrated architecture that combines a waveguide with a quantum-dot/photonic-crystal cavity QED platform. As such, they have made a series of pioneering contributions in the effect of homodyning for quantum engineering and optimization, which is promised to a great followup [50]. In the following, we will provide a unified theory of the mixing of coherent and quantum light that has been developed and implemented throughout the recent years by Fischer et al. [45,48]. The possibility and benefits of an external laser to optimize photon correlations also appeared in a work by Van Regemortel et al. [51], with a foothold in the same ideas. The effect of tuning two types of driving was emphasized by Xu and Li [52], who reported among other notable results how changing their ratio can bring the system from strong antibunching to superbunching, an idea which we will revisit from the point of view of interfering fields through (controlled) homodyning or (self-consistent) selfhomodyning. Similar ideas have then been explored and extended several times in many variations of the problem [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68] which all fit nicely in the wider picture that we will present. The microcavity-polariton configuration with interactions in one mode only (describing quantum well excitons, the other being a cavity mode) has been studied mainly from the (conventional) polariton blockade point of view [23,69], in which case, the (much stronger) unconventional antibunching has been typically overlooked. We will focus in the following on this case rather than on the possibly more popular two weakly-interacting sites. First, because this allows a direct comparison with the Jaynes-Cumming limit, second, because this configuration became timely following the recent experimental breakthrough with polariton blockade [31,32]. Whatever the platform, it needs be also emphasized that while many works have focused on singlephoton emission as the spotlight for the effect (which is dubious when antibunching is produced from the unconventional route), others have also stressed different applications or suggested different contextualization, such as phase-transitions [37] or entanglement [70], and there is certainly much to exploit from one perspective or another.
Among other configurations that cannot be accommodated by Eq. (1) as they add even more components, one could mention examples from works that involve additional modes (three in Refs. [71][72][73]), different types of nonlinearity, e.g., a 2 b † in Refs. [74][75][76], a four-level system in Ref. [77]), two two-level system in a cavity in Ref. [78] and up to the general Tavvis-Cummings model [79], pulsed coherent control of a two-level system in Ref. [80], two coupled cavities each containing a two-level system [81,82] up to a complete array [83]. It seems however clear to us that the phenomenology reported in each of these particular cases would fall within the classification that we will establish in the remaining of the text, i.e., they can be understood as an homodyning effect of some sort.

III. HOMODYNE AND SELF-HOMODYNE INTERFERENCES
We will return in the rest of this text to such systems as those discussed in the previous Section-all a particular case or a variation of Eq. (1)-to show that the twophoton statistics of their emission can be described to lowest order in the driving by a simple process: the mixing of a squeezed and a coherent state. In this Section, we therefore study this configuration in details.
We first consider the mixture of any two fields as obtained in one of the output arms of a balanced beam splitter (cf. Fig. 1a), which is fed by a coherent state on one of its arms, with a complex amplitude α = |α|e iφ , and another field of a general nature, described with annihilation operator d, on the other arm. The field that leaves the beam splitter is a mixture of the two input fields, whose annihilation operator can be written as s = α + d, where we are leaving out the normalization (1/ √ 2) and π 2 -phase shift in the reflected light as reasoned in Appendix A. Within this description, any normally-ordered correlator of the resulting field can be expressed in terms of the inputs as From this expression, we can compute any relevant observable of the mixture. For instance, the total population is with n d ≡ d † d. Apart from the sum of both input intensities, there is a contribution (last term) from the firstorder interference between the coherent components of each of the fields or mean fields. Similarly, the secondorder coherence function, which is defined as can be readily obtained from the correlators in Eq. (2). In this text, we omit the time t in all expressions, as it is considered to be large enough for the system to have reached the steady state under the continuous drive, and will also set the delay τ = 0, thus focusing on coincidences. This simplifies the notation g s (t → ∞, τ = 0). We will also consider N -th order coherence functions, also at zero delay: g These correlators can always be written as a polynomial series of powers of the amplitude of the coherent field α: where c k (φ) are coefficients that depend on the phase of the coherent field φ, and mean values of the type d †µ d ν with µ+ν ≤ N 2 . In particular, the 2nd-order correlation function, Eq. (4b), can be rearranged as with I m ∼ |α| m [28,[84][85][86], where 1 represents the coherent contribution of the total signal, and the incoherent contributions read Here, the notation "::" indicates normal ordering,  for the interference between a coherent and a quadraturesqueezed field, as in the setup shown in panel (a). The resulting g (2) s is shown in panel (b) as a function of the coherent field |α| and squeezing parameeter r (colour map in log scale). The relative phase between φ between the squeezed and coherence state is the one that optimises antibunching, i.e., φ = θ/2. (c) The n-norm (up to g (6) , i.e., n = 5) as defined in Eq. (14). Dashed black lines in both panels mark the minimum of g (2) s , showing that the best antibunching is no guarantee of good two-photon emission. (d) Cut of g (2) s along the horizontal dashed line in panel (b) (|α| = 0.3) and its decomposition given by Eqs. (6) and (7). The black dotted line shows g namely, to root nonclassical effects observed in optical bistability with many atoms in a cavity, to the physics of a single atom coherently driven, i.e., resonance fluorescence. This is at this occasion of unifying squeezing and antibunching from two seemingly unrelated platforms under the same umbrella of self-homodyning that he introduced this terminology. These concepts have thus been invoked and explored to some extent from the earliest days of the field, but it is only recently that they seem to start being fully understood and exploited [40,44,45,48,51]. Using such interferences to analyse the squeezing properties of a signal of interest (here, the field with annihilation operator d) through a controlled variation of a local oscillator (here, the coher-ent field α), was first suggested by W. Vogel [85,86], subsequently implemented to show squeezing in resonance fluorescence [87], and recently impulsed in a series of works from the Vucković group, as previously discussed. The physical interpretation of the contributions to the decomposition are as follows: • The numerator of I 0 is the normally-ordered variance of the signal intensity, that is, :(∆n d ) 2 : = :n 2 d : − n d 2 with n d = d † d and ∆n d = n d − n d . Therefore, I 0 < 0 indicates that the field d has sub-Poissonian statistics, which in turn contributes to the sub-Poissonian statistics of the total field s.
• The numerator of I 1 represents the normallyordered correlation between the fluctuation-field strength and intensity, d † d 2 − d † d d = :∆d ∆n d : , which have been referred to as anomalous moments [85,86] and been recently measured [88]. A squeezed-coherent state has such correlations.
• The numerator of the last component, I 2 , can also be written in terms of the quadratures of the field d.
Having I 2 < 0 necessarily implies that the state of light has a squeezing component. This can be proved by noting that :X 2 d : = X 2 d − 1/4 (the same for :Y 2 d :). Then, rearranging the numerator of Eq. (7c) leads to: where X d,φ = (e iφ d † + e −iφ d)/2 is the quadrature with the same phase of the coherent field (given by the angle φ). If I 2 < 0, the dispersion of X d,φ must be less than 1/2, but since X d,φ and its orthogonal quadrature, X d,φ+π/2 , must fulfil the Heisenberg uncertainty relation, ∆X d,φ ∆X d,φ+π/2 ≥ 1/4, then ∆X d,φ+π/2 > 1/2. This necessarily implies that there is a certain degree of squeezing in d. Nevertheless, the opposite statement is not true. A state with a non-zero degree of squeezing can have I 2 ≥ 0, for instance, if the relative direction between the coherent and squeezing contributions fulfils θ − 2φ = π/2 (a straightforward example is provided by the displaced squeezed state). Furthermore, if d = 0, the numerator of I 2 simplifies to 4|α| 2 ( :X 2 d : − |α| 2 ).
An analogous procedure allows us to decompose the third-order coherence function, g  As an illustration which will be relevant in what follows, let us consider the interference between a coherent and a squeezed state, as shown schematically in Fig. 1(a). The coherent state can be written as |α = D a (α) |0 , where α = |α|e iφ is the amplitude of the coherent state as before, and D a (α) = exp αa † − α * a is the displacement operator of the field with annihilation operator a. Similarly, the squeezed state may be written as |ξ = S d (ξ) |0 , where ξ = re iθ is the squeezing parameter and S d (ξ) = exp ξd † 2 − ξ * d 2 is the squeezing operator of the field with annihilation operator d. Thus, the state that feeds the beam splitter is |ψ in = |α , ξ . The interference at the beam splitter mixes these two states, and the state of the light that leaves the beam splitter is a two-mode squeezed state that is further squeezed and displaced (the detailed transformation is given in Appendix A). Since we are only interested in the output of one of the arms of the beam splitter, we take the partial trace over the other arm and we end up with a displaced squeezed thermal state [27], where now the displacement and squeezing operators correspond to the operator s = a + d, and n th , the thermal population, can be obtained from the population of the squeezed state, n d , which for a balanced beam splitter follows the relation n th = ( 1 + n d − 1)/2. The second-order correlation for the total signal is computed as in Eq. (4), taking the averages using the state in Eq. (9), namely g (2) s = Tr ρ s s † 2 s 2 / Tr ρ s s † s 2 , which can be decomposed as in Eq. (6) into where n s = |α| 2 + sinh 2 (r). Here, d = 0 but also Eq. (10b) is exactly zero because, for a squeezed state, the correlators d †µ d ν vanish when µ + ν is an odd number. Useful expression for the decomposition of g (2) s and g (3) s in terms of the incoherent component and the two-photon coherence are given in Appendix A.
Inspection of Eqs. (10) shows that the only way in which g (2) s < 1, that is, the statistics of the total signal can be sub-Poissonian, regardless of the value of the squeezing parameter, is for I 2 to be negative, which implies that the phases of the displacement and the squeezing must be related by |θ − 2φ| < π/2. We take for simplicity the minimizing alignment, θ = 2φ, which means that the coherent and squeezed excitations are driven with the same phase, since the phase of the squeezed state is θ/2. Using such a relation, the interference yields the correlation map shown in Fig. 1(b) as a function of the amplitude of the coherent |α| and squeezing r intensities. The black dashed line shows the optimum amplitude of the coherent state that minimizes g (2) s for a given squeezing, which is given by |α| min = e r cosh(r) sinh(r) .
This goes to zero although at the same time as the population goes to zero. Figure 1(d) shows a transverse cut of the correlation map in (b) along the purple long-dashed line, which corresponds to |α| = 0.3. The decomposition and total g (2) s are shown as a function of the squeezing parameter, with minimum g (2) s = 0.26 at r ≈ 0.078. Without the interference with the coherent state, the squeezed state can never have sub-Poissonian statistics. In fact, in such a case the correlations become independent of the phase of the squeezing parameter: which diverges at vanishing squeezing r → 0 (with also vanishing signal n d = sinh 2 (r)), and is minimum when squeezing is infinite r → ∞.
There is a great tunability from such a simple admixture since g (2) s of the light at the output of the beam splitter can be varied between 0 and ∞ simply by adjusting the magnitudes of the coherent field and the squeezing parameter. In particular, the most sub-Poissonian statistics occurs when coherent light interferes with a small amount of squeezing r < |α min |, in the right intensity proportion, given by Eq. (11). Counter-intuitively, g (2) s 1 occurs when the squeezed light itself is, on the opposite, super-Poissonian (even super-chaotic g (2) d > 2) [89]. This is a fundamental result that we will find throughout the text in order to find the conditions for and manipulate sub-Poissonian statistics and antibunching in various systems under weak coherent driving.
An important fact for our classification of photon statistics is that, since the sub-Poissonian behaviour is here due to an interference effect, the set of parameters that suppresses the fluctuation at the two-photon level does not suppress them at all N -photon levels, which means that the multi-photon emission cannot be precluded simultaneously at all orders. In other words, the condition in Eq. (11) that minimizes g (2) s , also minimizes the two-photon probability in the interference density matrix 2| ρ s |2 , at low intensities. But this is not the same condition that minimizes any other photon probability n| ρ s |n . This incompatibility is revealed by the n-norm, as defined in Ref. [90], which is the distance in the correlation space between signal s and a perfect single-photon source: In Figure 1 s , which lies in a highfluctuation region when the higher order correlation functions are taken into account. Further increasing n renders the correlation map completely red which means that multiphoton emission is not suppressed even if we have g (2) s close to zero. This is a feature typical of antibunching that arises from a two-photon interference only, that suggests that their use as single-photon sources may be an issue in the context of applications for quantum technology where higher-photon correlations may jeopardize two-photon suppression. This is related to the fact that this antibunching stems from a Gaussian state, which is the most classical of the quantum states.
The discussion presented above and, in particular, the decomposition of the second-order correlation as in Eq. (6), are not limited to the particular case of interfering pure states set as initial conditions. This can also be applied to the dynamical case of a single system which provides itself and directly a coherent compo-nent α along with another, and therefore quantum, type of component. Calling s the annihilation operator for a particular emitter which has such a coherent-but not exclusively-component in its radiation, one can thus express its emission as the interference (or superposition) of a mean coherent field s and its quantum fluctuations, with operator d = s − s . That is, one can always write Following the terminology previously introduced in the literature for a similar purpose [28], we call this interpretation of the emission a self-homodyne effect. Since g (2) s is also given by Eq. (7) with the simplification brought by the fact that d = 0, by replacing α → s and d → s − s , we obtain the general expressions in terms of s †n s m for the emission of a single-emitter s, interfering its own components: For completeness, in Appendix C we present the possible models (Hamiltonians and Liouvillians) that produce a coherent state and a squeezed state in a cavity, and give the correspondence between the dynamical parameters (such as the coherent driving and the squeezing intensity) and the abstract quantities α and r.
In the following Sections we will use the self-homodyne decomposition Eqs. (16) and this understanding in terms of interferences between coherent and quantum components, which in our cases may or may not be of the squeezing type, to analyse some statistical properties (anti-and super-bunching) of systems in their lowdriving regime, and contrast their statistics with conventional blockade effects. We will focus on cases that are both fundamental and tightly related to each other, namely, the 2LS (resonance fluorescence) in Section IV, the anharmonic oscillator in Section V, the Jaynes-Cummings Hamiltonian in Section VI and microcavity polaritons in Section VII.

IV. RESONANCE FLUORESCENCE IN THE HEITLER REGIME
We first consider the excitation of a two-level system (2LS) driven by a coherent source in the regime of low excitation-commonly referred to as the Heitler regime. Such a system is modeled by the Hamiltonian ( = 1) This is the particular case of the general Hamiltonian (1) when only one mode is considered and U → ∞. Here, the 2LS has a frequency ω σ and is described with an annihilation operator σ that follows the pseudospin algebra, whereas the laser is treated classically, i.e., as a complex number, with intensity Ω σ (taken real without loss of generality) and frequency ω L . The dynamics only depends on the frequency difference, ∆ σ ≡ ω σ − ω L . The dissipative character of the system is included in the dynamics with a master equation describes the decay of the 2LS at a rate γ σ . The steadystate solution (computed as indicated in Appendix D) can be fully written in terms of two parameters: the 2LS population (or probability to be in the excited state) n σ ≡ σ † σ , and the coherence or mean field α ≡ eff (Eq. (E16)) when the driving is low enough. The red solid line indicates when the effective squeezing parameter fits properly the actual statistics and becomes dashed when the approximation fails. In both cases, the range of validity is Ωσ 0.1γσ. where As a consequence of the fermionic character of the 2LS, it can only sustain one excitation at a time. Therefore, all the correlators different from those in Eq. (19) vanish, and in particular the N -photon correlations of the twolevel system are exactly zero, namely g (N ) σ = 0 for N ≥ 2. We call this perfect cancellation of correlations to all orders conventional blockade or conventional antibunching (CA), as it arises from the natural Pauli blocking scenario. To investigate the components of the correlations that ultimately provide the perfect sub-Poissonian behaviour of the 2LS, we separate the mean field from the fluctuations of the signal (σ = α + ), in analogy with Eq. (15). Following Eqs. (16), g (2) σ can be decomposed as in Eq. (6) with, These are presented in Fig. 2(a) as a function of the intensity of the driving laser. The decomposition shows that, although the photon correlations of the 2LS are always perfectly sub-Poissonian, or antibunched [92], the nature of their cancellation varies depending on the driving regime [44]. In the high-driving regime, the coherent component is compensated by the sub-Poissonian statistics of the quantum fluctuations (I 0 < 0) since lim Ωσ→∞ α = 0 and fluctuations become the total field, → σ. In contrast, in the Heitler regime the coherent component is compensated by the super-Poissonian but also squeezed fluctuations (I 2 < 0). The Heitler regime is, therefore, an example of the type of self-homodyne interference that we discussed in Sec. III.
Let us then analyse more closely the fluctuations by looking into their correlation functions: where n = n σ −|α| 2 is the contribution from the fluctuations to the total population of the 2LS. More details are given in Appendix E. In particular, the N -photon correlations from the fluctuations alone is given by, that in terms of the physical parameters reads In Fig. 2(b) we plot g (2) confirming that fluctuations are sub-Poissonian or super-Poissonian depending on whether the effective driving defined as Ω eff ≡ Ω σ / 1 + (2∆ σ /γ σ ) 2 is much larger or smaller than the system decay γ σ , respectively (the figure is for the resonant case).
In the Heitler regime, we need to consider only the magnitudes up to leading order in the effective normalised driving p ≡ 2Ω eff /γ σ . The main contribution to the intensity n σ = |α| 2 + n , in the absense of pure dephasing, comes from the coherent part |α| 2 of the signal. Fluctuations only appear to the next order, having, up to fourth order in p: The coherent contribution corresponds to the elastic (also known as "Rayleigh") scattering of the laser-photons by the two-level system, while the fluctuations originate from the two-photon excitation and re-emission [93]. In the spectrum of emission, this manifests as a superposition of a delta and a Lorentzian peaks with exactly these weights, |α| 2 and n , both centered at the laser frequency, with no width (for an ideal laser) and γ σ -width, respectively [44,94,95]. Fluctuations have no coherent intensity by construction, = 0. At the same time, their second momentum is not zero but exactly the opposite of the coherent field one: 2 = −α 2 , thanks to the fact that σ 2 = α 2 + 2 = 0. This means that both contributions, coherent and incoherent, are of the same order in the driving p when it comes to two-photon processes and can, therefore, interfere and even cancel each other. This is precisely what happens and is made explicit in the g (2) σ -decomposition above. The strong twophoton interference (I 2 ) can compensate the Poissonian and super-Poissonian statistics of the coherent and incoherent parts of the signal (1 + I 0 ). Since quadrature squeezing is created by a displacement operator, or a Hamiltonian, based on the operator 2 , this situation corresponds to a high degree of quadrature squeezing for the fluctuations. Further analysis of the fluctuation quadratures (details of the calculation are in Appendix E) shows that their variances behave similarly to a squeezed thermal state in the Heitler regime, which allows us to derive an effective squeezing parameter r eff = p 2 to describe the state to lowest order in the driving. This is plotted in Fig. 2(b) as a function of the driving, with the line becoming dashed when the interference can no longer be described in terms of a squeezed thermal state. Note that the total signal has no squeezing at low driving, only fluctuations do, because the coherent contribution is much larger.
Resonance fluorescence by itself always provides antibunching, due to the perfect cancellation of the various components. However, one can disrupt this by manipulating the coherent fraction, simply by interfering the signal σ in a beam splitter with an external coherent state |β . This allows to change the photon statistics of the total signal s = Tσ + iRβ, where T 2 and R 2 are the transmittance and reflectance of the beam splitter. Actually, since the decomposition affects correlators to all orders, Eq. (22), one can target the N -photon level instead of the 2-photon one. Namely, one can decide to set the N -photon coherence function to zero. As a particular case, the 1-photon case cancels the signal altogether, which is obtained by solving the condition (2) s as a function of both F and φ (in units of π) with the colour code in inset (logscale). Tuning the external laser allows to choose between various resonant conditions. n s = T 2 |α + iRβ 1 /T| 2 = 0 (because n = 0 to second order in Ω σ ). We will show that the possibility to target one N in isolation of the others introduces a separate regime from conventional blockade. Given their relationship and in line with the terminology found in the literature, we refer to this as unconventional blockade and unconventional antibunching (UA).
With this objective of tuning N -photon statistics and in order to avoid referencing the specificities of the beam splitter which do not change the normalized observables, let us define β ≡ Rβ/T ≡ |β |e iφ and parametrise its amplitude as a fraction F (always a positive number) of the laser field exciting the 2LS: With this, the coherence function g (N ) s of the interfered field in the Heitler regime is given by: Since g (1) s = 1, this expression also provides the population n s by considering the case N = 1. One can appreciate the considerable enrichment brought by the interfering laser by comparing n s (with the interfering laser, F = 0) to Eq. (19(a)) (without, F = 0) and even more so by comparing the N -photon correlation function, which is identically zero without the interfering laser, and that becomes Eq. (26) with the interfering laser. Interestingly, there is now another condition that suppresses the correlations and yields perfect antibunching at a given N -photon order, in addition to the one obtained in the original system without the interfering laser (CA). The new conditions exist for any detuning and are given by: Focusing on the resonance case for simplicity, we have F N = 2N and always the same phase, φ N = π, which corresponds to the field iβ N = −N |α|. The total coherent fraction changes phase for all N : s ) vanishes due to a first order (or one-photon) interference at the external laser parameter F 1 = 2, which translates into iβ 1 = −i|α|. The external laser completely compensates the coherent fraction of resonance fluorescence, in this case α = i|α| (with |α| = 2Ω σ /γ σ ). This situation corresponds to a classical destructive interference, which equally occurs between two fully classical laser beams with the same intensity and opposite phase.
In the case of highest interest, that of two-photon correlations g (2) s , we find a destructive two-photon interference at the intensity F 2 = 4, which corresponds to an external laser iβ 2 = −2|α|, that fully inverts the sign of the coherent fraction in the total signal: α + iβ 2 = −α. This coherent contribution leads to perfect cancelation of the two-photon probability in a wavefunction approach [96] (see the details in Appendix F 2). Note that this does not, however, satisfy all other N -photon interference conditions and g (N ) s with N > 2 do not vanish. This is a very different situation as compared to the original resonance fluorescence where g (N ) σ = 0 for all N > 1. One, the conventional scenario, arises from a an interference that takes place at all orders. The other, the unconventional scenario, results from an interference that is specific to a given number of photons.
All these interferences can be seen in Fig. 3(a) where we plot them up to N = 4. When there is no interference with the external laser, F = 0, antibunching is perfect to all orders recovering resonance fluorescence. At the onephoton interference, the denominator of g (N ) s becomes zero and the functions, therefore, diverge. This produces a superbunching effect of a classical origin, as previously discussed: a destructive interference effect that brings the total intensity to zero. In this case, the external laser removes completely by destructive interference the coherent fraction of the total signal. Therefore, the statistics is that of the fluctuations alone, what we previously called g (N ) , given by Eq. (22). We have already discussed how, in the Heitler regime, fluctuations become super chaotic and squeezed. We can see, on the left hand side of Fig. 1, that in the limit of Ω σ → 0, they actually diverge. Such a superbunching is thus linked to noise. The resulting state is missing the one-photon component and, consequently, the next (dominating) component is the two-photon one. Nevertheless, there is not a suppression mechanism for components with higher number of photons so that the relevance of such a configuration for multiphoton (bundle) emission remains to be investigated, which is however better left for a future work. We call this feature unconventional bunching (UB) in contrast with bunching that results from a N -photon de-excitation process that excludes explicitly the emission of other photon-numbers. This superbunching, as well as the antibunching by destructive interferences, will be reappearing in the next systems of study. The Heitler regime is, therefore, a simple but rich system where all the squeezing-originated interferences already occur although we need an external laser to have them manifest.
We now turn to the subtle point of which quantum state is realized by the various scenarios. To lowestorder in the driving, the dynamical state of the system can be described by a superposition of a coherent and a squeezed state, insofar as only the lower-order correlation functions (namely, population and g (2) ) are considered. This is shown in Table I, where we compare g (N ) for 1 ≤ N ≤ 3 (with N = 1 corresponding to the population) for the fluctuations in the Heitler regime vs the corresponding observables for a squeezed thermal state, on the one hand, and the laser-corrected configuration vs the displaced squeezed thermal state on the other hand. As explained in Appendix E, such a comparison can be made by identifying the squeezing parameter and thermal populations to various orders in a series expansion of the quantum states with the corresponding observables from the dynamical systems. One finds: where we have defined Γ 2 σ = γ 2 σ + 4∆ 2 σ . By definition, fluctuations have a vanishing mean, i.e., = 0 so we must choose α = 0. On the other hand, for the corrected emission, since one is blocking the two-photon contribution (at first order, this gives g (2) = 0), the comparison with a displaced thermal state is obtained by imposing the condition for g (2) to vanish at first order (r = |α| 2 and θ = 2 φ). The results are compiled in the table up to the order at which the results differ. Through the typical observables that are the population and g (2) , one Squeezed Thermal Heitler fluctuations Displaced Squeezed Thermal Laser-corrected configuration I. Two-level system. Comparison of first-(population), second-and third-order photon correlations, i) between a Squeezed Thermal state and the fluctuations in the Heitler regime and ii) between a Displaced Squeezed Thermal state and the fluctuations in the laser-corrected Heitler regime, to various orders in the driving Ωσ.
g (N ) Displaced Squeezed Thermal AO antibunching Displaced Squeezed Thermal Laser-corrected configuration Anharmonic oscillator. Comparison of first-(population), second-and third-order photon correlations, i) between a Displaced Squeezed Thermal state and the anharmonic oscillator with ∆ b = ∆− (optimal antibunching) and ii) between a Displaced Squeezed Thermal state and the laser-corrected for the optimal g (2) configuration (F2,2 and φ2,2). Selected parameters: can see how the system is indeed well described to lowest order in the driving by a coherent squeezed thermal state (displaced if there is a laser-correction). However to next order, there is a departure, showing that the Gaussian state representation is an approximation valid up to second-order only. In fact, for three-photon correlations, the disagreement occurs already at the lowestorder in the driving, and is of a qualitative character, as is also shown in the table. Therefore, such a description is handy but breaks down if a high-enough number of photons or a too high-pumping is considered.

V. ANHARMONIC BLOCKADE
To show that the effects of conventional (selfhomodyne interference at all N ) and unconventional (self-homodyne interference at a given N only) blockades take place in a general setting and are not specific of strong quantum nonlinearities (such as a 2LS), we now address the case of a single anharmonic oscillator, that describes an interacting bosonic mode with a Kerr-type nonlinearity, which can be very weak. With driving by a coherent source (a laser) at frequency ω L , its Hamiltonian reads where the cavity operators are represented by b † and b, ∆ b = ω b − ω L is the detuning between the cavity and the laser, U denotes the particle interaction strength (that provices the nonlinearity) and the driving amplitude is given by Ω b . This is the particular case of the general of an anharmonic oscillator, as a function of the detuning ∆ b From light to dark blue, U increases (the exact values are U/γ b = 0.01, 0.1, 1, 2, for fixed γ b = 1). A zoom-in of the smallest case is shown in the inset.
Hamiltonian (1) when only one mode is considered and U remains finite and, generally, small. The level structure of this system (at vanishing driving) is given by the simple expression E (N ) = N ω b + N (N − 1)U . The condition for the laser frequency to hit resonantly the N -photon We restrict our analysis of the dynamicsρ = −i [H ao , ρ] + (γ b /2)L b ρ, with γ b the decay rate of the mode, to the case of vanishing pumping, i.e. Ω b γ b . Solving the correlator equations in this limit gives the population the 2nd-order Glauber correlator as well as the higher-order correlators This shows that, when scanning in frequency, g has two extrema, one minimum and one maximum, as can be seen in Fig. 4, whose positions are given by with respective optimum antibunching (−) and bunching (+) Both of these features are linked to the level structure: the antibunching condition is that of resonantly driving the first rung, E (1) (note that ∆ − ∼ 0, especially when U γ b ), and the bunching condition, that of driving the second rung, E (2) (∆ + ∼ −U/2). In both cases, all other rungs are off-resonance and will remain much less occupied. Therefore, these effects are of a conventional nature, as we have defined it in the previous section: CA and conventional bunching (CB), respectively. The difference with resonance fluorescence is that here, CA is not a perfect interference at all orders (g (N ) = 0 for N > 1) but an approximated one. For instance, g (2) ≈ (γ b /U ) 2 (to leading order in U/γ b ), is only zero in the limit U → ∞, when the system converges to a 2LS. On the other hand, there was not CB in resonance fluorescence due to the lack of levels N > 1. Here, we see it appearing for the first time.
The decomposition g  together with its decomposition, as a function of the non-linear interactionnonlinearity strength U following the conventional antibunching frequency ∆ b = ∆−(U ). and negative (for ∆ b < −U/2) values, resulting in super-Poissonian statistics or on the contrary favouring antibunching. The various terms and the total g (2) b are shown in Fig. 5 as a function of U , for the case ∆ b = ∆ − (U ) that maximizes antibunching, showing the evolution from Poissonian fluctuations in the linear regime of a driven harmonic mode to antibunching as the two-level limit is recovered with I 0 → 1 and I 2 → −2 (cf. Fig. 2(a)).
As previously, the statistics can be modified by adjusting the coherent component of the original signal b with an external laser β = |β|e iφ . The resulting signal is then described by the operator s = Tb + i Rβ, with coherent contribution now given by s = T b + i Rβ. To simplify further the calculations, we choose β = T R β , where β is also written in terms of the driving amplitude and an adimensional amplitude F: Additionally, we shift the phase φ → φ + π so in the limit of high U , all the results are consistent with the previous case. Then, the total population becomes: where we have usedΓ 2 Here as well, we can compare the enrichment brought by the interfering laser by comparing Eqs. (30) and (37) for populations and Eqs. (31) and (38) for second-order correlations, with and without the interference, respectively. In this case, higher-order correlators could also be given in closedform but is too awkward to be written here. The cases g (k) s for 2 ≤ k ≤ 4 are shown in Fig. 6 as a function of the parameters of the interfering laser. By comparing this to Fig. 3 for the 2LS, one can see that the anharmonic system is significantly more complex, with resonances for the correlations that occur for specific conditions of the phase for each F that leads to unconventional forms of antibunching or superbunching, rather than to be simply out-of-phase previously. This makes salient the punctual character of the unconventional mechanism: each strong correlation at any given order must be realized in a very particular way: the one that matches the corresponding interference.
The maximum bunching (UB) accessible with the interfering laser is reached when the coherent-fraction population goes to zero (1-photon suppression) for which the conditions on the phase and amplitude read Those conditions are exactly the same as Eq. (27) for N = 1. Analogous conditions for the multi-photon cases can be found solving g (N ) s = 0. For the case N = 2, we find four different roots: Since these should be, by definition, real, this imposes another constrain on φ. Although the expression for real-valued φ to make F real cannot be given in closed form, they are readily found numerically. It is possible to get four real solutions, that are however degenerate. There are only two different conditions for φ since the real part is the same for each pair of roots, i.e., Re(F 2,1 ) = Re(F 2,4 ) and Re(F 2,2 ) = Re (F 2,3 )). This yields two physical solutions. For instance, for U = γ b and ∆ b = ∆ − (the case shown in Fig. 6), g (2) s vanishes at F 2,1 ≈ 0.615 and φ 2,1 ≈ 0.659 π for one solution and at F 2,2 ≈ 2.907 and φ 2,2 ≈ 0.860 π for the other one. Similar resonances in higher-order correlations could be found following the same procedure.
Regarding the quantum state realized in the system, similar conclusions can be drawn for the anharmonic oscillator than for the two-level system (previous Section, cf. Tables I and II). Specifically for this case, the system can be described by a displaced squeezed thermal state, properly parameterized, but to lowest-order in the driving and for the population and the two-photon correlation only. Departures arise to next-order in the pumping or to any-order for three-photon correlations and higher. The main differences is that the anharmonic oscillator case has to be worked out numerically, so the prefactors are given by the solutions that optimize the antibunching, for the system parameters indicated in the caption. The same result otherwise holds that the Gaussian-state description is a low-driving approximation valid for the population and two-photon statistics. The same holds for the systems studied in the following Sections, although this point will not be stressed anymore.

VI. JAYNES-CUMMINGS BLOCKADE
Now that we have considered the two-level system on the one hand (Section IV) and the bosonic mode on the other hand (Section V), we turn to the richer and intricate Physics that arises from their coupling. We will show how the themes of the previous Sections allow us to unify in a fairly concise picture the great variety of phenomena observed and/or reported in isolation.
We thus consider the case where a cavity mode, with bosonic annihilation operator a and frequency ω a , is coupled to a 2LS, with operator σ and frequency ω σ , with a strength given by g. Such a system is described by the Jaynes-Cummings Hamiltonian [97,98], where we also include both a cavity and a 2LS driving term by a laser of frequency ω L , with respective intensities Ω a and Ω σ and relative phase φ. We assume g and Ω σ to be real numbers, without loss of generality since the magnitudes of interest (G (N ) a ) are independent of their phases. The relative phase φ between dot and cavity drivings is, on the other hand, important. We also limit ourselves in this text to the case where the frequencies of the dot and cavity drivings ω L are identical and the analysis could be pushed further to the case where this limitation is lifted. The dissipation is taken into account through the master equation ∂ t ρ = i [ρ, H jc ] + (γ σ /2)L σ ρ + (γ a /2)L a ρ where γ a is the decay rate of the cavity. We solve the steady state in the low-driving regime, i.e., when Ω a γ a , γ σ , as indicated in Appendix D and find the populations: with matching upper/lower indices (including ±) and withΓ 2 = a, σ). Similarly, we find the two-photon coherence function from the cavity: where∆ ij ≡ i∆ a +j∆ σ ,γ ij = iγ a +jγ σ ,Γ 2 ij ≡γ 2 ij +4∆ 2 ij and χ = Ω σ /Ω a is the ratio of excitation. The range of χ extends from 0 to ∞ so that it is convenient to use the derived quantityχ = 2 π atan(χ) which varies between 0 and 1. Equation (43) is admittedly not enlightening per se but it contains all the physics of conventional and unconventional photon statistics that arises from selfhomodyning, including bunching and antibunching, for all the regimes of operations. It is quite remarkable that so much physics of dressed-state blockades and interferences can be packed-up so concisely. We plot a particular case of this formula as a function of ω a and ω L in Fig. 7(a), namely, only driving the cavity (Ω σ = χ = 0). Its reduced expression, Eq. (I1), is given in Appendix I. The general case is available through an applet [99] and we will shortly discuss other cases as well. The structure that is thus revealed can be decomposed in two classes, as shown in panel (b): the conventional statistics that originates from the nonlinear properties of the quantum levels, in solid lines, and the unconventional statistics that originates from interferences, in dashed lines. Both can give rise to bunching (in red) and antibunching (in blue). We now discuss them in details.

A. Conventional statistics
Conventional features arise from the laser entering in resonance with a dressed state of the dissipative JC ladder [100,101], which energy is the real part of The first rung E ± yields the CA lines in Fig. 7(b). This corresponds to an increase in the cavity population, as shown in Fig. 7(c) as two white lines, corresponding to the familiar lower and upper branches of strong coupling. The system effectively gets excited, but through its first rung only. The second rung blocks further excitation according to the conventional antibunching (CA), or photon-blockade, scenario, so that with the increase of population goes a decrease of two-photon excitation, leading to antibunching. This is in complete analogy with the CA that appears in the Heitler regime of resonance fluorescence. This is not an exact zero in g (2) a in the low driving regime (the imaginary part of the root does not vanish) because the conditions for perfect interference are no longer met having a strongly coupled cavity with a decay rate. We recently showed in Ref. [44] that even in the vanishing coupling regime, g → 0, when the cavity acts as a mere detector of the 2LS emission, perfect antibunching is spoiled due to the finite decay rate (γ a representing the precision in frequency detection). This is due to the fact that the cavity is effectively filtering out some of the incoherent fraction of the emission while the coherent fraction is still fully collected. The interference condition in the g (2) a decomposition, 1 + I 1 = −I 2 = 2, is no longer satisfied (see Fig. 2 of Ref. [44]).   On the other hand, driving resonantly the second rung, E ± , leads to conventional bunching (CB), shown as red solid lines in Fig. 7(b). These quantum features are well known and also found with incoherent driving of the system in the spectrum of emission [101], they are not conditional to the coherence of the driving. This also corresponds to an increase in the cavity population although this is not showing in Fig. 7(c), where only first order effects appear.

B. Unconventional statistics
We now turn to the other features in g (2) a that do not correspond to a resonant condition with a dressed state: these are, first, a superbunching line at ω L = 0 (dashed red in Fig. 7(b)) and second, two symmetric antibunched lines (dashed blue). All correspond to a self-homodyne interference that the coherent field driving the cavity can produce on its own, without the need of a second external laser. In this case, this also involves several modes (degrees of freedom) and more parameters than in resonance fluorescence, so the phenomenology is richer, but can be tracked down to the same physics. We call them again unconventional antibunching (UA) and unconventional bunching (UB) in full analogy with the Heitler regime of resonance fluorescence and in agreement with the literature that refers to particular cases of this phenomenology as "unconventional blockade" [26] (the term of "tunnelling" has also been employed but the underlying physical picture might be misleading [10]).
We first address antibunching (UA). This is found by minimizing g (2) a in regions where there is no CA, which yields (for the particular case χ = 0): which is the analytical expression for the the dashed blue lines in Fig. 7(b) (we remind that ∆ i ≡ ω i − ω L for i = a, σ). The most general case when both the emitter and cavity are excited is given in Appendix H. In the minimization process, we also find the condition for CA, due to the first-rung resonance, but this can be disconnected from UA beyond the fact that CA is already identified, because of UA also admits an exact zero, which is found by either solving g (2) a = 0 or setting to zero the two-photon probability in the wavefunction approximation [96] (see Appendix F 2). This gives the conditions on the detunings as function of the system parameters [102]: These conditions are met in Fig. 7(a) at the lowest point where the blue UA line intersects the (e) cut (and on the symmetric point ω a < 0). When the laser is at resonance with the 2LS (∆ σ = 0) and cavity losses are large (γ a γ σ ), this occurs when the cooperativity C ≡ 4g 2 γaγσ = 1. This type of UA interference is secondorder, so it is not apparent in the cavity population at low driving, Fig. 7(b). One has to turn to two-photon correlations instead. Note also that UA requires a cavityemitter detuning that is of the order of g.
Since this is an interference effect, we perform the same decomposition of g (2) a in terms of coherent and incoherent fractions, as in previous Sections, given by Eq. (6), and show the terms that are not zero in Fig. 7(d-e). The full expressions are in Appendix G. The term I 1 is exactly zero to lowest order in the driving and only the fluctuation-statistics I 0 and the two-photon interference I 2 play a role, like in the Heitler regime of resonance fluorescence. We observe that, in this decomposition, there is no difference between the CA and UA, since both occur approximately when the statistics of the laser and fluctuations, 1 + I 0 = 2, are compensated by their two-photon interference, I 2 = −2, again as in the Heitler regime. The fundamental differences between these two types of antibunching will be discussed later on. Before that, we discuss the last feature: the unconventional bunching at ω L = 0.
The reason for the super-bunching peak labelled as UB in Fig. 7(b) is also the same as in resonance fluorescence: the cancellation of the coherent part, in this case, of the cavity emission, and the consequent dominance of the fluctuations only, which are super-Poissonian in this region. Therefore, contrary to the CB, this superbunched statistics is not directly linked to an enhanced N -photon (for any N ) emission and it does not appear one could harvest or Purcell-enhance it, for instance, by coupling the system to an auxiliary resonant cavity. Since it is pretty much wildly fluctuating noise, the actual prospects of multi-photon physics in this context remains to be investigated. In any case, the conditions that yield the super-Poissonian correlations can thus be obtained by minimising the cavity population n a or, from the wavefunction approximation detailed in Appendix (F 2), by minimising the probability to have one photon, given by Eq. (F12), which coincide with the coherent fraction to lowest order in Ω a . One cannot achieve an exact zero in this case but the cavity population is clearly undergoing a destructive interference, as shown by the black horizontal line in Fig. 7(c). The resulting condition links the laser frequency with the 2LS one: which reduces to simply ∆ σ = 0 (laser in resonance with the 2LS) if i) the dot and cavity are driven with a π/2phase difference or ii) the laser drives the cavity only (χ = 0). So far, we have focused on the particular case of Eq. (43) where Ω σ = 0 (i.e., Eq. (I1)). This is the case dominantly studied in the literature and the one assumed to best reflect the experimental situation. It is also for our purpose a good choice to clarify the phenomenology that is taking place and how various types of statistics cohabit. It must be emphasised, however, that while the physics is essentially the same in the more general configuration, the results are, even qualitatively, significantly different in configurations where the two types of pumping are present. This is shown in Fig. 8. While conventional features are stable, being pinned to the level structure, the unconventional ones that are due to interferences are very sensible to the excitation conditions and get displaced or, in the case of QD excitation only, even completely suppressed. If one is to regard conventional features as more desirable for applications, this figure is therefore again an exhortation at focusing on the QD excitation configuration.
While we have focused on the two-photon statistics, both the conventional and unconventional effects occur at the N -photon level, in which case they manifest through higher-order coherence functions g (N ) , and their N thorder behaviour is one of the key differences between conventional and unconventional statistics. Regarding conventional features, resonances happen at the N -photon level when N photons of the laser have the same energy than one of the dressed states (and only one, thanks to the JC nonlinearities): ω L = Re{E (N ) ± }/N . The blockade that is realised is a real blockade in the sense that all the correlation functions are then depleted simultaneously. In Fig. 9, the counterpart of Fig. 7(a) is shown for g (3) a and g (4) a and shows how more conventional features appear with increasing N but otherwise stay pinned to the same conditions, while the number of unconventional features stays the same, but their positions drift with N , so that one cannot simultaneously realise g (N ) a < 1 for all N . This is an important difference between a convex mixture of Gaussian states, which is a semi-classical state, and a state beyond this class, which is genuinely quantum, as previously mentioned. The latter requires the ability to imprint strong correlations at several and possibly all photon-orders. This suggests that CA could be more suited than UA for quantum applications. Note how with the 2LS direct excitation, shown in the second row of Fig. 9, one only finds conventional statistics, with magnified features such as broader antibunching in the photon-like branch and narrower one in the exciton-like branch. The N -photon resonances are neatly separated for large-enough detuning, which is the underlying principle to harness rich N -photon resources [103].
We now turn to another noteworthy regime, out of the many configurations of interest that are covered by Eq. (43), namely, the transition from weak to strong coupling. The so-called strong-coupling, when g > |γ a − γ σ |/4, is one of the coveted attributes of light-matter interactions, leading to the emergence of dressed states and to a new realm of physics. It is also, however, an ill-defined concept in presence of detuning [101] and one would still find the dressed-state structure of Fig. 9 in the largely detuned regime when driving the 2LS, even up to large photon-order [22]. The restructuration of the statistics when crossing over to the weak-coupling regime is explored in Fig. 10(a), where we track the impact on g (2) a of changing the coupling g, on the cut in Fig. 7(e) that intersects from top to bottom CB, UA (twice), UB and CA. One can see how the features converge as the coupling is reduced, with the conventional ones disappearing first, which is expected from the disappearance of the underlying dressed states, that are responsible for the conventional effects. The unconventional antibunching, on the other hand, is more robust and can be tracked well into weak coupling where all effects ultimately vanish at the same time as they merge. Conventional antibunching is the most robust feature, as can be seen by tracking for instance the UB peak at the point where it is the most isolated from the other features, namely, at resonance where ω L = ω a = ω σ = 0. Spanning over the two main parameters that control strong coupling, the coupling strength g (in units of γ σ ) and the rates of dissipation rates γ a /γ σ , one sees that the strong bunching is not always sustained but can be instead overtaken by conventional antibunching, which is the well-defined blue line in the figure (given by Eq. (45)). The region where the UB peak is well-defined can be identified by inspecting the second derivative of g at ω L = 0 and is shown in Fig. 10(b) as a dashed black line. The white line that separates the antibunching region from the bunching one corresponds to the critical coupling strength g P between the cavity and the 2LS that leads to g (2) a = 1 (its expression is given in the Appendix, Eq. (I3)). The strong-weak coupling frontier g/γ σ < |γ a /γ σ − 1|/4 is indicated with a dotted green line as a reference, illustrating again the lack of close connection between strong-coupling and the  photon-statistics features.
We conclude our discussion of the Jaynes-Cummings system with the second main difference between conventional and unconventional statistics, namely their resilience to higher driving. All our results are exact in the limit of vanishing driving, that is to say, in the approximation of neglecting Ω terms of higher-orders than the smallest contributing one. For non-vanishing driving, numerically exact results can be obtained instead (and can be made to agree with arbitrary precision to the analytical expressions, as long as the driving is taken low enough, what we have consistently checked). A characteristic of the unconventional features is that, being due to an interference effect for a given photon number only, it is fragile to driving, unlike the conventional features which display more robustness. This is illustrated in Fig. 10(c) for the case of cavity driving Ω a , where we compare the analytical result from Eq. (43) or, in this case, Eq. (I1), in black, with the numerical solution for Ω a = 0.25γ a , so still fairly small. One can see how the conventional features are qualitatively preserved and quantitatively similar to the analytical result, while the unconventional antibunching has been completely washed out.
One could consider still other aspects of the Physics embedded in Eq. (43). We invite the inquisitive and/or interested readers to explore them through the applet [99] which is helpful to get a sense of the complexity of the problem. Instead of discussing these further, we now turn to another platform of interest that bears many similarities with the Jaynes-Cummings results.

VII. MICROCAVITY-POLARITON BLOCKADE
Microcavity polaritons [22] arise from the strong coupling between a planar cavity photon and a quantum well exciton, both of which are bosonic fields with annihilation operators a and b, respectively. These fields are coupled with strength g and have frequencies ω a and ω b . Moreover, the excitons, being electron-hole pairs, have Coulomb interactions that we parametrise as U/2. Thus, the Hamiltonian describing the polariton system is given by where ∆ a,b = ω a,b − ω L are the frequencies of cavity/exciton referred to ω L , which is the frequency of the laser that drives the photonic/excitonic field with amplitudes Ω a,b . The phase difference between them is represented by φ = φ a −φ b and since the absolute phase can be chosen freely, φ a,b are fixed to be φ and 0. The dissipative dynamics of the polaritons is given by a master equa- where γ a and γ b are the decay rates of the photon and the exciton, respectively. As compared to the Jaynes-Cummings Hamiltonian (41), the polariton Hamiltonian substitutes the 2LS by a weakly interacting Bosonic mode, b → σ with nonlinearities b † b † bb thus slightly displacing the state with two excitations while the 2LS forbids it entirely. In the case where U → ∞, the Jaynes-Cummings limit is recovered, but in most experimental cases, U/γ a is very small. In all cases, in the low driving regime (Ω a → 0) the steady-state populations of the photon and the exciton are given by the same expressions as in the Jaynes-Cummings model, Eq. (42) with σ → b, since the 2LS converges to a bosonic field in the linear regime. The differences arise in the two-particle magnitudes (cf. Eq. (43)): where we have used the short-hand notation γ + = γ a + γ b , ∆ ± = ∆ a ± ∆ b and Γ 2 c = γ 2 c + 4∆ 2 c for c = a, b, + as well as ∆ ij ≡ i∆ a + j∆ b ,γ ij = iγ a + jγ σ ,Γ 2 ij ≡γ 2 ij + 4∆ 2 ij ,Ũ ij = iU + j∆ b andj denotes negative integer values (j = −j).
Note that a major conceptual difference with the Jaynes-Cummings model is that it now becomes relevant to consider the emitter (in this case, excitonic) two-photon coherence function, g b , while in the Jaynes-Cummings case one has the trivial result g (2) σ ≡ 0. The exciton statistics enjoys noteworthy characteristics, as we shall shortly see.
We repeat in Fig. 11 the same plots for the polariton system as for the Jaynes-Cummings case (Fig. 7). The applet [99] also covers this more general case. The cavity population is exactly the same, as already mentioned, and all other panels bear clear analogies. The two-photon coherence function converges to the Jaynes-Cummings one in the infinite interation limit (lim U →∞ g (2) a ) but is distinctly distorted for high-energy laser driving in the positive photon-exciton detuning region, and features an additional UA and CB couple of lines in the negative detuning region. The decomposition of g (2) a as in Eq. (6) can be made (the expressions are however bulky and relegated to Appendix G) and are shown in Fig. 11(d) and (e).

A. Conventional statistics
Like in the Jaynes-Cummings model, one can identify the conventional antibunching (CA) and bunching (CB) by mapping the observed features to an underlying blockade mechanism, namely, the positions at which N -photon excitation occurs, which is when the laser is resonant with one of the states in the N -photon rung. The first rung that provides CA is given by the same Eq. (44), with N = 1, since this corresponds to the linear regime where both systems converge. One finds, therefore, that the two CA blue lines in Fig. 11(a), marked in solid blue in (b), are the same as in the Jaynes-Cummings model. They coincide as well with the white regions in Fig. 11(c) where the cavity emission is enhanced. This is the standard one-photon resonance, with a blockade of photons into higher rungs due to the non-linearity now introduced by the interactions (instead of the 2LS).
Higher rungs are different from the Jaynes-Cummings model, but their effects otherwise follow from the same principle and they are similarly obtained by diagonalizing the effective Hamiltonian in the corresponding N -excitation Hilbert subspace, that is, in the basis {|N, 0 , |N − 1, 1 , . . . |0, N } (where each state is characterised by the photon and exciton number). At the two-photon level, one is interested in the second rung, which contains three eigenstates. The expressions for the general eigenenergies are rather large but we can provide here the first order in the interactions U in the strong coupling limit (g γ a , γ b ): with R = g 2 + (ω a − ω b ) 2 /4 the normal mode splitting typical of strong coupling. In this limit, E (1) The CB lines are positioned, therefore, according to the conditions for two-photon excitation by the laser: in increasing order, as they appear in Fig. 11(a), marked with solid red lines in Fig. 11(b). The upper CB line, corresponding to E (2) + , is the faintest one in the cavity emission due to the fact that it has the most excitonic component. It is monotonically blueshifted with increasing U and becomes linear in the density plot as E

B. Unconventional statistics
We now shift to the unconventional features in polariton blockade. Superbunching, or UB, is found by minimization of n a and, therefore, also occurs for the same condition as the Jaynes-Cummings model ∆ b = χg cos φ. The maximum superbunching is obtained at one of the crossings of UB and CB. Now turning to the more interesting Unconventional Antibunching (UA) features, they are found, in the polariton case as well, from the minimization of g (2) a . Since the equations are quite bulky, only the case of cavity excitation (Ω b = 0) is included here (the full expression can be consulted at App. H). The UA curve is given by the solution of a and (b) its structure in terms of conventional and unconventional features. CA (CB) is given by the resonant condition with the first (second) polariton manifold, cf. Eqs. (44) (for N = 1) and (50). UA is given by the interference condition Eq.   8.63 g). The decomposition in Ij components is also shown, with the same conventions as in Fig. 7. Parameters: γa = 0.1 g , γ b = 0.01 g and U = 10γa. and the conditions for perfect antibunching come from solving the equation and subsequently imposing that every parameter must be real (or the more restrictive case: real and positive) that lead to additional restrictions. b . In the latter case, all the unconventional features have disappeared. The bottom row identifies the features through the structure of conventional and unconventional lines. Parameters: γa = 0.1 g, γ b = 0.01 g, χ = 0 and U = 10 γa.
As already noted, the polariton case adds a third CA line as compared to its Jaynes-Cummings counterpart. The correspondence between both cases is still clear, but this is largely thanks to the large interaction strength chosen in Fig. 11, namely, U/γ a = 10. This choice will allow us to survey quickly the polaritonic phenomenology based on the more thoroughly discussed Jaynes-Cummings one. Figure 12, for instance, shows the polaritonic counterpart of Fig. 8 on its left panel but for one case of mixed-pumping only, highlighting the considerable reshaping of the structure and the importance of controlling, or at least knowing, the ratio of exciton and photon driving. The right panel of Fig. 12 provides g (2) b , which, if compared to Fig. 11, shows that the main result is to remove all the unconventional features and retain only the conventional ones. The peaks are also less in the excitonic emission, producing a smoother back- ground. Another dramatic feature of the excitonic correlations, which is apparent from Eq. (49), is that it is independent from the ratio χ of driving, i.e., the same result is obtained if driving the cavity alone, the exciton alone, or a mixture of both, in stark contrast of the cavity correlations (cf. Fig. 11(a) and 12(a) where the only difference is that half the excitation drives the 2LS in the second case rather than going fully to the cavity in the first case). This could be of tremendous value for spectroscopic characterization of such systems since it is typically difficult to know the exact type of pumping, while experimental evidence shows that both fields are indeed being driven under coherent excitation [104]. If measuring the excitonic correlations, there is no dependence from the particular type of coupling of the laser to the system. On the other hand, excitonic emission is much less straighftorward of access. Note, finally, that one could similarly consider the lower and upper polariton statistics, but they are even less featureless, with correlations of the signal that merely follow the polariton branches (their expression is given in the appendix for completeness). Finally, in Fig. 13, we focus on the effect of the interaction strength and how to optimise the observation of antibunching. We have already emphasized that for clarity of the connection between the Jaynes-Cummings and the polariton case, we have considered a value of U/γ a substantially in excess even of the few cases themselves largely in excess of the bulk of the literature [105][106][107]. While it is not excluded that such a regime will be available in the near future, it is naturally more relevant to turn to the most common experimental configuration where U/γ a 1. We show such a case in Fig. 13(a), where U/γ a = 0.1. We see how, as a result, the CB and CA lines of the positively-detuned case collapse one onto the other. The UA line previously in between has, in the process, disappeared. The CA and CB however do not cancel each other but merge into a characteristic dispersive-like shape, shown in panel (b), the observation of which, predicted over a decade ago [23], has been a long-awaited result for polaritons and has indeed been just recently reported from two independent groups [31,32]. While this shape has been regarded as an intrinsic and fundamental profile of polariton blockade, our wider picture shows how it arises instead from different features brought to close proximity by the weak interactions. The difficulty in reporting polariton blockade lies in the weak value of antibunching, which is largely due to the fact that no optimisation over the full structure of photon correlations, that was unknown till this work, has been made for the driving configuration that yields the best antibunching for given system parameters.

VIII. SUMMARY, DISCUSSION AND CONCLUSION
We have connected an hitherto disparate and voluminous phenomenology of photon statistics in the light emitted by a variety of optical systems, into a unified picture that identifies two classes of conventional and unconventional features, covering both the cases of antibunching and bunching, which leads us to a classification of CA, UA, CB and UB. One class (conventional), linked to real states repulsion, occurs at all orders and for all photon numbers while the other (unconventional) occurs for a given photon-number with no a priori underlying level structure.
To lowest order in the driving, the dynamical response can be described by an interferences between a squeezed component and a coherent component, and thus, in this picture, one can understand the photon statistics emitted by many optical systems as simply arising from the particular way each implementation finds to produce some squeezing on the one hand and some coherent field on the other hand, and interfere them during its emission. In agreement with the previous literature, we call this phenomenon "self-homodyning". With this understanding, one can bring considerable tailoring of photon cor-relations by modifying the relative importance of coherence vs squeezing, which is conveniently achieved by superimposing a fraction of the driving laser to the output of the system ("homodyning"). This Gaussian-state approximation holds to lower orders in the driving for the first correlation functions only. For instance, for the case of the laser-corrected two-level system (Table I), the squeezed-coherent Gaussian description holds up to the second-order in the driving for the population and to the first-order in the driving for g (2) . Deviations occur for these observables to higher-orders in the driving while higher-order correlation functions already differ to the lowest order in the driving. Such deviations seem to arise from the non-Gaussian nature of the quantum fluctuations in these highly non-linear systems. This remains to be studied in detail.
Such a general picture can explain under a unified mechanism a wealth of observations that could otherwise appear to be peculiarities that are specific to a particular configuration. To take one recent example from a group that has been leading in the development and applications of the type of homodyning and self-homodyning that we have studied above, in Ref. [79], Trivedi et al. study the generalization of the Jaynes-Cummings system to N emitters: the so-called Tavis-Cummings Hamiltonian. Here, it is found that driving resonantly the eigenstates [108] produces conventional antibunching, flanked by unconventional antibunching for laser frequencies detuned from the one-and two-photon resonances. This is the counterpart of the situation of Fig. 7(d) (resonance) and (e) (detuning), both also shown in Panel (a), where increasing N has the effect of bosonizing the interacting (matter-like) part of the system or decreasing the effective nonlinearity, similarly to decreasing g for N = 1. Interestingly, it is reported that while for the case of resonance, antibunching is spoiled with an increasing number of emitters N , in presence of a detuning, one of the antibunching peaks is, on the opposite, enhanced with increasing N . This apparently puzzling behaviour is easily understood once the conventional and unconventional nature of the respective antibunching lines are recognized. In the resonant case, antibunching is always conventional, and as such it is spoiled by the bosonization of the system due to it increasing number of emitters [109], or by reducing the coupling. Since both weaken the nonlinearity in the level structure, this destroys the conventional blockade that is based on it. With detuning, on the other hand, one finds not only conventional but also unconventional antibunching, cf. Fig. 7(b). Their CA is also spoiled with increasing N , as reported, but their UA, however, increases, which can be expected since it is due to an self-homodyning interference between the coherent and incoherent parts of the emission at the two-photon level, as explained above, and this does not suffer from a reduced nonlinearity (or increasing N ). It can in fact be also optimized (i.e., reduced) like all types of UA and as a result, should even reach g (2) = 0 to lowest order for a proper choice of the detuning, that will depend on N in a way that remains to be computed. Since we have shown, however, that the interference nature of UA makes it sensitive to dephasing, and that detuning [110] results in fast oscillations in autocorrelation times, with a narrowing plateau of antibunching, one can also expect this antibunching to be particularly fragile and difficult to resolve when including a realistic model for its detection. This is consistent with the finding of the Authors that inhomogeneous broadening quickly spoils UA. Finally, they also find in both detuned and resonant cases the unconventional bunching, as the large bunching central peak that is a typical feature of the general mechanism (cf. Fig. 7). This is therefore the super-chaotic noise due to self-homodyning stripping down the emission to its mere fluctuations. As such, the interpretation in terms of two-photon bound states that is offered in Ref. [79] and in other works [10,47] should be further analyzed and quantified. We suspect the emission in UB to be less efficient for multiphoton Physics as compared to leapfrog emission [103], due to the lack of a suppression mechanism for higher photon-number processes, and despite the large values of the correlation functions that they produce.
As a conclusion, our picture brings considerable simplification in the interpretation and identification of the various phenomena observed in a plethora of systems, in particular with respect to connecting them between each other throughout platforms. It clarifies the value but also the limitation of a description in terms of Gaussianstates. It also opens a new route to control and fine-tune such photon correlations and make a more informed and better use of them for quantum applications.
(2) s, min < 1/2, it is required r < log ( √ 6 − 1)/2 ≈ 0.186, which implies g σ (τ ) = 1 for a continuously driven system, so that g (2) σ = 0 also implies g (2) σ < g Thus, for either low squeezing signal (r 1) or almost symmetrical BS (T − R ≈ 0), the output signal could be described as follows. The commutator [S o , S s ] vanishes for any possible values, so the exponential simplifies into S d (ξ) = S o (ξ o ) S s (ξ s ) S os (ξ os ). With the previous results, the output state can be written as: where |ξ os represents a two-mode squeezed state. In the Fock basis, this state can be written as (from Ref. [111], Chapter 7): where r os = |ξ os | = RT r. The corresponding density matrix for this pure state reads ρ out = |ψ out ψ out |. Tracing out output o, we obtain the density matrix for output s only (our signal of interest): ρ s = Tr o {ρ out }. For the next step, we will use the properties of the trace to move clockwise the operators acting over the output o subspace and where1 o is the identity operator. Furthermore, any operator that only acts on the s-subspace can be taken out of the trace. With this, Computing the partial trace: from the second to the third equality we have used m|p o = δ m,p and p|n o = δ p,n . The resulting density matrix has the form of a thermal state ρ th with mean population p th ≡ s † s th = sinh 2 r os . To sum up, the output field detected at a single arm of the system corresponds to a displaced squeezed thermal state. So, to sum up, with parameters α s = T|α|e iφ , ξ s = r s e iθs = R 2 e i(θ+π) and p th = sinh 2 (RT r). Even though T and R appear as parameters, last equation only works for R ≈ T. We restrict ourselves to the case of 50:50 beam splitter (T 2 = R 2 = 1/2). Thermal population can be expressed in terms of squeezed population of the input signal n d = sinh 2 r: From ρ s we can compute the observables for the mixed signal: We provide the third-order coherence function of signal s, in terms of two interfering fields, s = α + d. As in the case for g (2) s , the highest order contributions in powers of α can be gathered into the coherent term (given by 1), yielding: and zero everywhere else. In the main text, we discuss first the case of resonance fluorescence, which corresponds to having only the 2LS operator σ and no cavity mode a (taking g, Ω a = 0 here). Second, we solve the Jaynes-Cummings model with both cavity and dot driving with a phase difference between the sources, which corresponds to setting Ω σ /Ω a = χe −iφ here. Similarly, for the polariton model where the matter field is an exciton (boson), we have: and, again, the remaining matrix elements are zero. These equations can be solved numerically, by choosing a high enough truncation in the number of excitations, in order to obtain the steady state (∂ t C {m,n,µ,ν} = 0) for any given pump. However, we are interested in an analytical solution when applying vanishing driving limit (Ω σ → 0 or Ω a → 0). In this case, it is enough to solve recursively sets of truncated equations. That is, we start with the lowest order correlators, with only one operator, that we write in a vectorial form for convenience (using the JC model as an example): v 1 = ( a a † σ σ † ) T . Its equation, ∂ t v 1 = M 1 v 1 + A 1 + h. o. t., provides the steady state value v 1 = −M −1 1 A 1 + h. o. t., to lowest order in Ω a (with h. o. t. meaning higher order terms). We proceed in the same way with the two-operator correlators v 2 = ( a 2 a †2 a † a σ † σ σ † a . . .) T , only, in this case, we also need to include the steady state value for the one-operator correlators as part of the independent term in the equation: We, in particular, aim at obtaining photon correlators of the type a †N a N that follow a †N a N ∼ (Ω a ) 2N , to lowest order in the driving Ω a . The normalized correlation functions g (N ) a = a †N a N / a † a N are independent of Ω a to lowest order, their computation requiring to solve the 2N sets of recurrent equations and taking the limit lim Ωa→0 g (N ) a .
Appendix E: Homodyne interference with resonance fluorescence: correlations and squeezing from the fluctuations and the total signal The correlations from the fluctuations of resonance fluorescence, with operator = σ − α, can be accessed using the technique of homodyne detection explained in Sec. (III). In this case, we feed one of the beam splitter arms with resonance fluorescence (d → σ) and the other with a coherent field (a → β). The correlators of the output of the arms as defined in Appendix A, s = iRβ + Tσ, are given by: Since σ †p σ q = 0 for p, q > 1, this simplifies to For instance, the coherent fraction and total population of the output field are: Therefore, we can chose the coherent field to compensate exactly the coherent component of the 2LS α = σ , i.e., setting β = i T R σ so that we hace only the transmitted fluctuations s = T . In such a case the correlators simplify even further, †n m = s †n s m /T n+m = (−1) m+n α m−1 α * (n−1) nm n σ − (n + m − 1) α 2 .
With this general expression, we obtain the population and coherence functions, Eq. (22) of the main text. We want to recover and analyse light properties from the original 2LS (before the beam splitter). In order to do this, the factor T 2 on the population should be eliminated (making T 2 n σ → n σ , T σ → σ ). Note that this change is not inconsistent, given that the beam splitter divides by two the incoming signal and thus merely attenuating it by a factor T 2 . Moreover, since the photon correlations are normalized objects, a global attenuation in the unnormalized correlators do not change them. However, suppressing the coherent contribution of the emission is not the only possibility. We can also tune the coherent contribution by choosing β = e iφ |β |, where the amplitude is parametrize as |β | = R T |β|. Thus, we are broadening the range of possible output configurations [44]. So, for the most general case, the N -particle correlators have the following form: After this general expression, the amplitude |β | is usually expressed in a more suitable way, referencing it to the driving intensity of the laser: |β | = Ωσ γσ F. Other two important quantities are the mean and variance of the quadratures, X s,χ = 1 2 (e iχ s † + c.c.) and their dispersion ∆X 2 s,χ = X 2 s,χ − X s,χ 2 , respectively: The mean value only depends on the total coherent contribution s = T(iβ + α).
The maximum and minimum of the (normal-ordered) quadrature variance for a single-mode can be inferred independently of specific nature of the field: : ∆X 2 s : max/min = ∆X 2 s max/min − where the sign correspond to maximum and minimum, respectively. While the variance is strictly a positive-valued quantity, its normal-ordered counterpart is not. This latter indicates the deviation of the variance from the vacuum value (which is 1 4 ) so values below 0 will reveal some degree of quadrature squeezing. Likewise the angle of squeezing is generically given by: After substituting the correlators (E4) on (E7) and, then, using the steady-state solution given in Eq. (18): and the angle of squeezing will be: It is not surprising that factor F does not appear since all the squeezing properties exclusively come from the fluctuations. The strength of this effect is reduced by the factor T 2 as the input signal σ (α + ) is divided by the beam splitter, which can naturally absorbed into Ω 2 σ . Now we are interested in the low driving regime (Ω σ → 0), so the previous expressions at the lowest order in Ω σ simplify In this limit, both dispersion are symmetrical (but with opposite signs). We can regard these expression as a limit of low squeezing (and coherent intensity |α| 2 ) from a displaced squeezed thermal state (although the coherent part does not contribute to the variance). Such states has got the following dispersion when r → 0: : where the superscript on the average reminds that the observable corresponds to an displaced squeezed thermal state. We have approximated 1 + 2 n th as 1 since the thermal population grows like Ω 4 σ (which comes from the first order of the incoherent population). Comparing (E11) with (E12), we found that the incoherent population in the Heitler regime behaves like a squeezed thermal state with squeezing parameter and the effective thermal population p th : From these two parameters an effective g (2) , namely g eff , can be obtained for the fluctuations. Supposing that, in the low excitation regime, fluctuations would behave similar to an squeezed thermal state, then g (2) should have the same form. Fixing |α| = 0 in Eq. (A10b) and taking the limit r 2 → 0 and p th → 0 (both go to 0 with the same power dependence), we get which, after substituting Eqs. (E13)-(E14), reads Appendix F: Wavefunction approximation method at vanishing pumping regime In the context of this paper, the wavefunction approximations [96] consist of assuming that the state of the system composed by two fields, with annihilation operators ξ and c following either fermionic or bosonic algebra, can be approximated by a pure state, which in the Fock state basis reads, where C nm are the probability amplitude of having n photons in the field described with operator ξ and m photons in the field described with operator c; and the summation is done until the allowed number of photons: one for a fermionic field and N for a bosonic one. Given that the dynamics of the system is given by the master equation where H is the Hamiltonian of the system and we have assumed that the dissipation is given by "jump operators" j k at ratesΓ k , the dynamics of the wavefunction is given by Schrödinger equation where H eff is a non-hermitian Hamitonian constructed as H eff = H − i kΓ k j † k j k , and the coefficients evolve as, In the following sections we make explicit both the effective Hamiltonians and the differential equations for the coefficients for all the systems considered in the main text.
1. Two-level system in the Heitler regime The Hamiltonian describing the excitation of a sensor (a cavity) by the emission of a 2LS, which in turn is driven in the Heitler regime by a laser, is given by the Hamiltonian (41). To complete the analogy of beam splitter setting and be consistent with the main text, both driving and coupling for the sensor has to be defined in terms of coherent source amplitude |β| and BS parameters T and R: Ω a → iR|β|, g → Tg. The system and driving source are not necessarily at resonance so we define the detuning as ∆ σ = ω σ − ω L . The effective Hamiltonian that describes the dynamics in the wavefunction approximation is where H is the Hamiltonian in Eq. (41). Replacing the effective Hamiltonian in Eq. (F5) on the expression in Eq. (F3), we obtain the differential equations for the coefficients of interest: where we have assumed that the driving to the 2LS is low enough so that the states with three or more excitation can be safely neglected. Assuming that the coherent field that drives the sensor can be written as a fraction of the field that drives the two-level system (|β| = g TΩσ Rγσ F), very similar as Eq. (25), and to leading order in the coupling and the driving intensity of the two-level system, the solution to Eq. (F6) is [4γ 2 σ + F 2 e 2iφ (γ σ + 2i∆ σ ) (γ σ + Γ + 2i∆ σ ) − 4Fe iφ γ σ (γ σ + Γ + 2i∆ σ )] .
The population of both the 2LS and the cavity, and the g (2) a can be obtained from the coefficients in Eq. (F7). However, to recover some information from the unfiltered signal (Γ → ∞), this limit has to be performed carefully and the previous expressions need some manipulation first. A new set of coefficients is defined as: C ij = Γ 2Tg i C ij so any explicit dependence of sensor parameters disappears, resulting in a non-vanishing (finite) solution after the proper limit are taken. After the substitution and taking the limit, the new coefficients are Now, these solutions provide useful information about the equivalent filtered signal: n a ≈ |C 10 | 2 , P 20 = |C 20 | 2 (probability of two-photon state) and g (2) a ≈ 2|C 20 | 2 /|C 10 | 4 . The cancellation of the coefficient C 20 ,and therefore of g (2) a , yields the condition on the attenuation factor which can only be satisfied-F is an attenuation factor, and thus a real number-when the relative phase between the driving field and the 2LS coherent contribution is either 0 or π (opposite phase), in agreement with Fig. 3(a-d) in the main text. Note, as well, that the cancellation of the coefficient C 10 , and therefore of the population of the cavity, is obtained when F = 2e −iφ γσ γσ+2i∆σ which is a real number for the same phases for which Eq. (F9) is a real number.

Exciton-polaritons blockade
The Hamiltonian describing exciton-polaritons is given in Eq. (48), and the dynamics is complemented with a master equation that takes into account the decay rate of the cavity with rate γ a and of the excitons with rate γ b . Thus, the effective Hamiltonian that described the dynamics in the wavefunction approximation is Replacing the Hamiltonian in Eq. (F13) in Eq. (F4), we find that the differential equations for the relevant coefficients are i∂ t C 10 = e iφ Ω a + ∆ a − i γ a 2 C 10 + g C 01 i∂ t C 11 = e iφ Ω a C 01 + Ω b C 01 + ∆ a + ∆ σ − i γ a + γ σ 2 where we have assumed that the driving is low enough for the states with three or more photons to be neglected. The steady-state solution of the coefficients in Eq. (F14) is obtained when the derivatives in the left-hand side of the equation vanish. Thus, assuming that the coefficient of the vacuum dominates over all the others, i.e., C 00 ≈ 1, and to leading order in the driving of the cavity, the coefficients are C 01 = 2Ω a χ(2∆ a − iγ a ) − 2e iφ g 4g 2 + (γ a + 2i∆ a )(γ b + 2i∆ b ) , where we have used N = [4g 2 + (γ a + i∆ a )(γ b + ∆ b )][(γ a + i∆ a )(γ 11 + i∆ 11 )(U + 2∆ b − iγ b ) + 4g 2 (U + 2∆ 11 − iγ 11 )], and χ,γ ij and∆ ij share the same definition as described above changing σ by b.

Microcavity polariton
The perfect antibunching conditions can be derived straightaway from the equation g (2) a = 0, where g (2) a is given by the expression (I2). Then, we clear ∆ a from the previous equation, which leads to: where N is defined as: Nevertheless, by definition ∆ a must a real quantity. Taking the real part of this last expression leads to the equation for the curve which follows the UA effect. Moreover, the cancellation of its imaginary part (same as the JC model) turns out to be a second condition to exactly reach g a = 0, which can be cleared for any chosen parameter. Further analysis shows that more restrictions emerge from the fact of selecting only real-valued (physical) parameters. As an illustration, we present here the case of cavity excitation (χ = 0). Splitting both real and imaginary parts from Eq. (H4): First expression provides an implicit equation for the 3 distinct curves of UA shown in Fig. 11 whereas the latter gives the exact location where g (2) a vanishes.

. (I3)
A smaller coupling g < g P produces antibunched light while a larger coupling g > g P produces bunched light.