Engineering Gaussian states of light from a planar microcavity

Quantum fluids of light in a nonlinear planar microcavity can exhibit antibunched photon statistics at short distances due to repulsive polariton interactions. We show that, despite the weakness of the nonlinearity, the antibunching signal can be amplified orders of magnitude with an appropriate free-space optics scheme to select and interfere output modes. Our results are understood from the unconventional photon blockade perspective by analyzing the approximate Gaussian output state of the microcavity. In a second part, we illustrate how the temporal and spatial profile of the density-density correlation function of a fluid of light can be reconstructed with free-space optics. Also here the nontrivial (anti)bunching signal can be amplified significantly by shaping the light emitted by the microcavity.


Introduction
Generation and manipulation of nonclassical states of light has been at the heart of quantum optics ever since its development in the early days [1]. A widely applied criterion to quantify the nonclassicality of a state is provided by the intensity correlations of the electromagnetic field g (2) (0) = :Î 2 : / Î 2 , with ':' denoting normal operator ordering. Whenever g (2) (0) < 1, there is no classical analog of the quantum statistics and one is therefore led to conclude that the state has intrinsic quantum properties.
The photon blockade is one of the most general mechanisms to realize these antibunched photon states. A strong nonlinear spectrum of the cavity makes it energetically forbidden for a second photon to enter once there is a first photon inside [2]. To date, a plethora of platforms exists in which this nonclassical feature of a light field has been demonstrated: by using a trapped atom in a cavity [3], with quantum dots [4,5] and, more recently, in superconducting circuits [6,7]. See Refs. [8,9,10,11] for recent reviews on these topics.
Exciton-polaritons in planar microcavities, arising from the strong coupling between a cavity photon and a quantum-well exciton, provide a versatile platform to both generate and detect nontrivial photon features (see Refs. [12,13] for recent reviews). The statistics of the polaritonic field is directly accessible through the photons that escape from the cavity, making these systems particularly interesting for measuring photon correlations. Moreover, the properties of the fluid, such as the density, velocity and phase, can be directly manipulated by the external laser field.
In spite of these tremendous advantages, the major bottleneck of exciton-polariton systems has turned out to be the relatively weak interparticle interactions compared to currently achievable linewidths, ensuring that genuine quantum effects are often hidden by the strong incident laser field. Translated into the intensity fluctuations of the cavity output photon field, this means that g (2) (0) ≈ 1, the value for a coherent photon field. Nevertheless, albeit very small, there is a correction to this trivial value which stems from polariton-polariton interactions [14]. On the first level of approximation, this effect can be understood in terms of the creation of pairs of quantum fluctuations that propagate through the fluid in opposite directions. The statistics of these pairs of fluctuations is described in terms of a two-mode squeezed Gaussian state.
It was shown in a number of recent studies [15,16,17,18] that strongly antibunched photon statistics is not necessarily a consequence of a strong cavity nonlinearity. Certain setups consisting of two coupled cavities can exhibit a strong suppression of g (2) (0), despite having interactions that are substantially weaker than the typical dissipation rate. Dubbed the unconventional photon blockade, this remarkable effect was attributed to a destructive interference between two independent pathways for injecting a second photon in the cavity. In a later stage all these different ideas were somehow unified and clarified in the framework of optimally amplitude-squeezed Gaussian states [19]. Therefore, squeezed coherent states can display sizeable antibunching, provided the average value of the field and the squeezing are tuned in an optimal way [20]. The theoretically predicted antibunched signal was recently observed in an experiment with microwave photons on a superconducting circuit [21].
The crucial point here is the ability of the system to show a significant squeezing without the coherent part overwhelming the quantum fluctuations. In weakly interacting spatially extended systems, the main reason why quantum effects are weak is the dominant contribution of the coherent field over the fluctuations. For a spatially homogeneous system, however, the coherent field and quantum fluctuations are easily separated in Fourier space. Significant antibunching can then be achieved after a suitable attenuation of the k = 0 component of the light field. Previous proposals in the context of the unconventional photon blockade relied on a setup that consists of two coupled (0D) microcavities. The setup that we introduce in this work, illustrated in Fig. 1, utilizes the photons generated from quantum fluctuations in a 2D planar microcavity to engineer an interference between two squeezed intracavity modes, with opposite momentum (k, −k), and the coherent pumping field to produce the desired antibunched statistics.
The proposed setup consists of engineering an appropriate filtering scheme that exactly tailors single-mode states of the form discussed in Ref. [19]. With the introduction of a few linear optical elements, like beam splitters, phase shifters and attenuators, the optimal conditions for maximal suppression of g (2) (0) can be closely approached. We also compute the temporal profile g (2) (τ ) of delayed correlations to illustrate the expected sustain times of the antibunched signal, which is essentially limited by the photon lifetime.
While the cavity output field in the unconventional photon blockade is Gaussian within good approximation, it can still display highly nonclassical features due to the severely reduced intensity fluctuations. Consequently, the UPB provides an attractive mechanism for generating sizeable single photon sources, which can then be used as input for a computation scheme based on linear optics. However, the downside of the original proposal of two coupled cavities is the required fine tuning of system parameters like cavity coupling, photon nonlinearity and laser detuning [20]. In our scheme, this is in part circumvented by placing the interference and selection scheme after the microcavity, thus separating the squeezing and interference stage for increased control.
Since the antibunching in this setup originates from genuine particle interactions inside a planar microcavity, we devote the second part of this manuscript to investigating the spatial profile of correlations in the quantum fluid of light itself. By collecting and interfering all modes that escape from the microcavity, rather than isolating a single one, the real-space image of correlations can be reconstructed. Moreover, by carefully shaping the bundle of light before interference, it is possible to manipulate the (anti)bunched features in the spatialtemporal correlation pattern. Recently, the profile of density fluctuations has proven its importance in the context of analog gravity, as it encodes information about Hawking pair emission at a sonic hole horizon [22,23], as was recently investigated in a cold-atom experiment [24].
The structure of our paper is as follows. We start by giving in Sec.2 an overview of photon correlations in quantum fluids of lights in ideal planar geometries. We then continue with explaining how to implement the unconventional photon blockade by filtering out specific modes in Sec. 3. Finally, in Sec. 4, we illustrate how the spatial-temporal profile of densitydensity correlations can be reconstructed with free-space optics. Technical details on the Bogoliubov approximation and on the time-dependence of operators are given in Appendices A and B, while Appendix C is devoted to a brief discussion of the principal imperfections and noise sources that may disturb the unconventional blockade effect in realistic planar microcavity devices.

Quantum fluctuations in fluids of light
In this section, we give an overview of the main features of the photon dynamics in nonlinear planar cavities under a coherent and monochromatic illumination. While our model is fully general, our discussion will be focused on the experimentally most relevant case of semiconductor microcavities with strong light-matter coupling, where the elementary luminous particles, the so-called exciton-polaritons, have a mixed light-matter character. For these systems, a wide theoretical and experimental literature has appeared that investigates the photon dynamics from the many-body perspective of the quantum fluids of light [12,13]. The interested reader should not find any difficulty in transferring our results to different material platforms and to different frequency domains.
We first set the stage by introducing the model that will be used in the following sections, we then review the dynamics of small collective excitations around the steady state as introduced in [25] and finally we summarize the main features of quantum correlations within a linearized approximation, as they were first presented in Ref. [26].

The model
We consider quasi-resonantly driven photons in a planar microcavity irradiated by a monochromatic and spatially plane-wave laser with frequency ω L at normal incidence, as depicted in Fig. 1(a). After performing a unitary transformation to remove the time-dependence of the drive, the Hamiltonian in terms of the polariton field operator Ψ(r), with r ≡ (x, y), is found as (we set = 1 throughout) The polaritonic modes have a dispersion LP (k) with a cut-off frequency at LP (0) and a lowenergy behavior of the form LP (k) ≈ LP (0) + k 2 /2m with an effective mass m. Photons are quasiresonantly injected into the cavity at a constant rate by the coherent laser, represented by the drive term with amplitude F . The third-order optical nonlinearity is assumed to be defocusing and is described in the model as a contact potential with interaction constant g > 0. In the polariton case, this represents the repulsive interactions between quantum-well excitons. With typical values of the interaction constant of g ∼ 1µeV · µm 2 and polariton masses of m ∼ 10 −4 − 10 −5 m e (with m e the electron rest mass), we find that the dimensionless interaction constant mg ≈ 10 −4 1, so that the polariton system is clearly in the weakly interacting regime. This legitimates a linearized treatment of quantum fluctuations within a Bogoliubov approximation, as performed in the following.
Photons inside the cavity have a finite lifetime before they escape by tunneling through the cavity mirrors. The photonic losses can be described in the Born-Markov approximation by introducing them into the master equation for the density matrix in the form of a dissipator of the Lindblad form Here, τ = 1/γ is the polariton lifetime (typically of the order of 10 − 100 ps) andn =Ψ †Ψ . Equivalently, the system dynamics can be formulated in terms of a quantum Langevin equation [27,28] for the polariton field operatorΨ(r, t) . In this framework, the equation of motion reads [12] i∂ tΨ (r, where we have defined the detuning of the laser frequency from the bottom of the bare LP dispersion δ = ω L − LP (0). The non-unitary time evolution due to Markovian photon losses is represented by the imaginary loss term inside the brackets and by the noise operatorsξ(r, t), which assume Gaussian statistics. In the low-temperature regime (k B T ω L ) under consideration here, their variance is given by The quantum Langevin equation (4) will be the starting point of our analysis.

The linearized equations of motion
In the spatially uniform configuration considered here, we can conveniently parametrize the photon field operator aŝ Ψ(r, t) = ψ 0 (t) +φ(r, t) = ψ 0 (t) QWs DB DB Figure 1: A schematic image of the selection and interference scheme that we propose. (a) A sketch of a planar semiconductor microcavity in the strong light-matter coupling regime. A quantum well (QW) is placed between two distributed Bragg mirrors (DB) and is pumped with a coherent laser beam. Photon interactions, mediated by the optical nonlinearity of the QW, lead to a small depletion of the condensate in terms of quantum fluctuations with nonzero momentum. The excitations are formed in pairs with opposite momentum and form a two-mode squeezed Gaussian state. A fluctuation with momentum k can leave the cavity as propagating radiation at an angle ϑ k , with sin ϑ k = ck/ω L , from the perpendicular axis. (b) A schematic image of the setup, consisting of linear optical building blocks, that we propose to engineer squeezed coherent states as output. The coherent field of the pumping laser is attenuated and interfered with the two-mode squeezed Gaussian state of the quantum fluctuations to construct the photon fieldσ k . (c) The output stateσ k can be approximately parametrized as a squeezed coherent state. The squeezing parameter r k is momentum dependent and can be found through (29). The phase between squeezing (θ k ) and displacement (ζ) can be varied with the phase shifters from (b), as given in (30), allowing the output field to go from an amplitude to a phase squeezed state. (d) A HBT setup to measure photon correlations of the output state. Detecting simultaneous clicks of detectors 1 and 2 allows for the measurement of the delayed intensity fluctuations of the photon fieldσ k and gives access to quantity (31). 'BS' stands for a (50:50) beam splitter and 'PS' for phase shifter.
where the classical field ψ 0 represents the coherent component and the quantum fieldφ(r, t) describes the fluctuations around it. From a many-body perspective, the classical field ψ 0 can be seen as a condensate coherently created in the k = 0 mode of the cavity by the incident laser beam. Owing to polaritonpolariton interactions, the condensate is slightly depleted and a fraction of the condensate particles is converted into photons with nonzero momentum k, as reflected by the fluctuation operatorsφ k . In the many-body language, this corresponds to the so-called quantum depletion of the condensate [29].
As a first step to solve the quantun Langevin equation (4), we can perform a mean-field approximation where the creation of these quantum fluctuations is neglected and only the condensate mode is assumed to be populated. In the steady state, the condensate density n 0 = |ψ 0 | 2 can then be found as the solution of the polynomial equation with the interaction-renormalized detuning For given pumping conditions (defined by the detuning δ and the amplitude F ), Equation (8) can have either one or two stable solutions for the mean-field density n 0 . When δ < √ 3γ/2 the resonator is in the optical-limiter regime. Alternatively, when δ > √ 3γ/2 the system will exhibit a bistable behavior with a corresponding hysteresis curve [12]. For the analysis that follows, we will always implicitly determine F by fixing a value of n 0 .
Making use of the decomposition (7), an equation of motion for the fluctuation operatorŝ φ k can then be derived out of the quantum Langevin equation (4). We will restrict to the lowest-order approximation of the interaction term in (4), which yields a set of linear equations for the zero-mean fluctuationsφ k . Effects of higher-order interaction terms were recently studied in the context of Beliaev-Landau-type scatterings in the quantum fluid [30].
At our level of approximation, the motion of the linearized quantum fluctuations is determined by where the Bogoliubov matrix B k involves the usual interaction-induced off-diagonal coupling µ = gn 0 but a shifted single-particle dispersion The eigenvalues ±ω k of B k governing the dispersion relation of the linearized fluctuations have the same analytical form as the collective Bogoliubov excitations of a dilute Bose gas at equilibrium but an important difference arises from the modified single-particle dispersion (11), which can either be gapped and strictly positive for ∆ < 0 or display regions of negative values for ∆ > 0. (c) If ∆ > 2µ, the diffusive modes shift to higher momentum and form a ring in momentum space. The diffusive regions are marked with blue shades. For clarity, we plot ω k − γ/2, which represents the total decay rate of a Bogoliubov mode with momentum k.
Restricting to the most straightforward cases, a negative effective detuning ∆ < 0 is found in the δ < 0 optical-limiter regime or on the high-density branch of the hysteresis loop in the δ > √ 3γ/2 bistable regime. A positive effective detuning ∆ > 0 is instead found on the low-density branch of the hysteresis loop and is characterized by the Bogoliubov dispersion being purely imaginary in some regions, so that the dynamics of the elementary excitations is diffusive-like (see Appendix A). Note that this regime is only parametrically stable when 2µ < γ, as we explain in Appendix B. In contrast to the equilibrium case where general arguments impose that the Bogoliubov dispersion is gapless [29], here the ε k=0 = 0 condition is only recovered in the special case ∆ = 0, i.e. at the end-point of the high-density branch. For clarification, we show the quasiparticle dispersion for different mean-field regimes in Fig.  2.
As a consequence of the photonic losses, the usual time-dependent Bogoliubov transformation, employed to diagonalize the equations of motion (10), acquires a damped exponential time dependenceφ The expressions for the time-dependent Bogoliubov coefficients η k (t) and ζ k (t) in the various mean-field regimes are presented in Appendix B. It should be noted, though, that the usual interpretation of the Bogoliubov operators as bosonic quasiparticles breaks down for the case of a diffusive dispersion. This is discussed in detail in Appendix A. The freedom to tune ∆ in the setup (by merely changing the pump frequency ω L ) allows us to explore parameter regimes, inaccessible to the conservative dynamics of dilute Bose gases at equilibrium, that can lead to novel, exotic physics. Previous work in this context has addressed superfluidity features [25] and the related drag force of a driven-dissipative fluid flowing past a defect [31]. Remarkably, it was pointed out in [32] that the diffusive modes for ∆ > 0 can even give rise to a negative effective drag force.

Correlations in the steady state
From the stochastic equations of motion for the quantum fluctuations, presented in (10), we can derive equations for the evolution of the quadratic correlation functions. For this we define the momentum distribution n k = φ † kφ k and the pair correlation c k = φ kφ−k to find The steady state of these equations is readily evaluated by setting the left-hand side to zero and leads to the equal-time correlation functions in the stationary regime [26] n k = 1 2 The nonequal-time correlations can also be found by making use of the Bogoliubov transformation (13). Thanks to the stationarity of the state, only the relative time difference Based on the quantum regression theorem [27], the delayed correlation functions can be obtained by evolving the later operator over a time τ ≥ 0 with (13), which gives in terms of the equal-time correlations n k and c k from (16).

Antibunched emission from squeezed quantum fluctuations
After introducing the general theoretical framework, we can start discussing how strongly antibunched light can be obtained by shaping the output of a weakly nonlinear coherently driven planar microcavity. We start this section by reviewing the basic concepts of squeezed coherent states and derive theoretical bounds on the maximal amount of antibunching that can be obtained in the family of Gaussian states that we can engineer as output from the cavity. Following the lines of Ref. [19], we then propose a first example of a selection and interference scheme that is able to manipulate the output of a many-mode planar cavity and obtain antibunched light by letting three k modes symmetrically located at ±k and 0 to mutually interfere. We finally characterize the intensity fluctuations and the level of antibunching that this scheme is able to produce.

Basics of Gaussian states
In the most general case, a single-mode squeezed state of the modeâ is represented by the density matrixρ is the squeezing operator with ξ = re iθ . Here, ρ n th is assumed to be a thermal density matrix with mean population n th = tr[ρ n thâ †â ].
Conversely, a Gaussian state is entirely characterized by its mean value and second moments and there exists a one-to-one map [1] that allows to extract squeezing parameters ξ and n th from them, Additionally, the modeâ can be displaced with a coherent field α =ᾱe iζ , which leads to the new density matrixρ is the displacement operator. The displacement field is then found back fromρ α,ξ,n th by relating A genuine signature of the non-classical nature of a state of light is provided by the intensity fluctuations of the photon field. The correlation function withn =â †â can be shown to obey g (2) (0) ≥ 1 for any classical state of light. Consequently, a violation of this inequality is a manifest indication of quantum correlations in the photon state. It was shown in Ref. [19] how optimally amplitude-squeezed Gaussian states of type (21) can strongly violate the inequality. Specific relations were derived between displacement α, squeezing ξ and thermal density n th to attain the theoretical lower bound on g (2) (0). In general, the equal-time second-order correlation function of a displaced, squeezed Gaussian state of form (21) reads Here c =ce iθ and we have set θ = 2ζ, so that the squeezing takes place exactly in the amplitude quadrature, thus obtaining optimal antibunching conditions. When the second moments n and c of the fluctuations are fixed, we can vary the displacement field α to find optimal antibunching conditions from (24). After straightforward algebra we find that settingᾱ establishes optimal antibunching. The intensity fluctuations of this particular displacement field are then given by g

Engineering squeezed coherent output states
In the squeezing language of this section, the non-trivial quadratic correlations found in (16) reflect the fact that the two modes at ±k are in a two-mode squeezed state. Similar to the single-mode case, the two-mode squeezing operator is defined asŜ 2 (ξ) = exp [−ξâb + ξ * â †b † ], with in our caseâ ≡φ k andb ≡φ −k . The main idea of this section is to illustrate that the intrinsic squeezing of the quantum fluctuations in the oppositely propagating modes can be utilized to construct a strongly antibunched output beam. Our proposal is based on the fact that the emission angle ϑ k of photons escaping a planar cavity is directly related to their cavity in-plane momentum k via sin ϑ k = ck/ω L , where c is the speed of light in vacuum and ω L is the laser angular frequency. Quantum fluctuations with in-plane momentum k will therefore lead to emission of pairs of photons at angles ±ϑ k , which can conveniently be isolated and later interfered to obtain a single-mode squeezed state as show in Fig. 1a). We illustrate in a schematic way (see Fig. 1b), the setup we propose to achieve this goal. We create a coherent population of polaritons by shining laser light onto the sample. The laser is tuned to the k = 0 lower polariton energy so as to only excite resonantly k 0 polaritons. Quantum fluctuations with momentum k > 0 will lead to emission of pairs of photons at angles ±ϑ k . The proposed experiment consists of combining the photons originating from the ±k modes onto a 50:50 beam splitter. One simple way of realizing this interference uses lenses to image the Fourier plane on two pinholes (or equivalently on the cores of two single-mode fibers) that only transmit the desired Fourier components, and to recombine them on a free-space (or fiber) beam splitter. We also suggest that imaging the Fourier space directly on an optical fiber bundle or on a spatial light modulator that selectively transmits the chosen modes would provide a more compact and elegant experimental implementation of this interference. After recombining these two modes, the light is interfered with a coherent field obtained by suitably attenuating and phase-shifting the pumping laser. Photon correlations are finally measured in a standard Hanbury Brown and Twiss (HBT) setup.
The final output state after selection and interference is then found to be of the form with ζ and ϕ ± the accumulated relative phases in the arms. By evaluating its expectation valueᾱe iζ and its quadratic correlation functions, where n k and c k are given in (16), it is straightforward to see that the modeσ k can be regarded as a squeezed coherent state of form (21). Thanks to the k → −k symmetry of the setup, the two-mode nature of the squeezing results here in exactly the same statistics as expected for a single-mode squeezed state. We can next determine the effective squeezing parameter r k e iθ k and thermal density n th,k for the output modeσ k . After straightforward algebra, we find from expressions (19,20) the thermal density and squeezing of output stateσ k (27) n th,k = n k + 1 with n k and c k found within the Bogoliubov framework and given in (16). Additionally, we see that the squeezing phase is found as Therefore, the difference between θ k and the displacement phase ζ can be tuned by varying ϕ + , ϕ − and ζ with the phase shifters in the setup from Fig. 1b). In the schematic image of the squeezed state shown in Fig. 1c), this corresponds to rotating the ellipse, which allows to switch between an amplitude and a phase squeezed state. In Fig. 3(a-b) we show how the parameters r k and n th,k from expression (29) depend on the selected momentum k in the two cases of respectively negative (Fig. 3b) and positive (Fig. 3a) values of the interaction-renormalized detuning ∆ defined in (9). Since excitations become generally more particle-like and have a larger frequency at large k, we expect the squeezing parameter r k to drop to zero in this limit. However, also their number n th,k goes to zero in the same limit as it is generally less likely to excite fluctuations with higher momenta. We can anticipate at this point that the decay of n th,k leads, in principle, to an asymptotic perfect antibunching of the photon statistics in the output state.
While both n th and r monotonously drop to zero in the ∆ < 0 case (as can be observed in Fig. 3b), the presence of a set of diffusive modes in the ∆ > 0 case from Fig. 3a (see Appendix B) leads to a more versatile behavior. As is reflected by the peak in n th at nonzero momentum, these diffusive modes are parametrically amplified. In addition, these modes are also strongly squeezed, which we conclude from the enhanced squeezing parameter r for the same momentum values as the peak in n th .

Optimizing the antibunching
The freedom to vary at will the attenuation level and the phase shift experienced by the coherent laser field in the setup from Fig. 1b) permits to approach the condition (25) for optimal antibunching. A measurement of the temporal correlation of the intensity fluctuations of the output stateσ k (27) gives access to the quantity where n k (τ ) and c k (τ ) are given in (17) and the total phase η = ϕ + + ϕ − − 2ζ. In Fig. 1d we provide a schematic image of a Hanbury-Brown-Twist (HBT) setup to illustrate how this quantity is typically measured. Let us first start by analyzing g (2) (0) (i.e at zero time delay) for the case where the phases are tuned such that squeezing is exactly in the amplitude quadrature, as realized by setting This condition amounts to rotating the major axis of the ellipse in Fig. 1c) into the direction perpendicular to the displacement vector. In Fig. 3(c,d) we show how the optimal displacement field α opt (see (25)) and the minimal value of g (2) (0) (see (26)) depend upon the selected momentum k in the setup from Fig. 1b). At high momenta g (2) (0)| min always drops to zero, meaning that we can, in principle, approach a perfect antibunching by selecting higher momenta. This can be understood by noticing that 2 2 Figure 3: (a,b) The squeezing parameter r and the thermal occupation n th as a function of momentum k for a steady state with a positive (a) and a negative (b) interactionrenormalized detuning (8). (c,d) The optimal displacement amplitude α opt (21) and corresponding g opt = g (2) (0)| min (26) as a function of momentum for the same parameters as above.
(e,f) The temporal profile of g (2) (τ ) after selecting various momenta as indicated below. The displacement amplitudeᾱ and phase η have been chosen to fulfil the optimal antibunching conditions (as derived from (c,d)) at τ = 0. The insets show g (2) (0) upon varying the total phase between squeezing and displacement (corresponding to rotating the ellipse in Fig. 1c)) for the showed momenta.
g (2) (0)| min ≈ 8 √ n th for small n th in the case of optimal squeezing and displacement, as was explained in Ref. [19]. While we have experimental control to tune α to its optimal value, the squeezing parameter r is set by the nature of the nonlinear processes inside the cavity. We can verify that r = r opt in general, but, interestingly, we still find that g (2) (0)| min → 0 in the limit n th → 0. Note also that, unfortunately, α opt approaches zero in this limit as well, such that the total photon flux is expected to become vanishingly small.
For intermediate values of k, the behavior is different according to the sign of ∆. While in the ∆ < 0 case, a monotonously decaying behavior is observed in Fig 3d) for both α opt and g (2) (0)| min as a function of selected momentum k, the presence of parametrically amplified modes results in a strong increase of n th,k and r k in the ∆ > 0 case, shown in Fig 3c) . The overall effect on g (2) (0)| min is, however, detrimental, as this quantity is pushed back towards its value g 2) (0)| min ≈ 1 for a classical coherent field.
Let us now move to the full temporal dependence of g (2) (τ ), as expressed in (31). In Fig.  3(e-f) we analyze the delayed second-order correlations g (2) (τ ) for the displacement amplitudẽ α =α opt that optimizes g (2) (0) (see expression (25)) for different values of k/ √ mγ. In general we see from Fig. 3(e,f) that g (2) (τ ) shows a damped oscillatory behavior with an oscillation frequency set by ω k (i.e. the Bogoliubov frequency of the selected Fourier component), and a damping equal to γ/2 − ω k (i.e. the lifetime of the associated Bogoliubov mode). Moreover, we can shift the offset of the oscillation by varying the phase η, which we illustrate in the insets of Fig. 3(e-f). In the main panels (a,b), the offset of g (2) (τ ) has been chosen such that maximal antibunching is achieved at τ = 0, corresponding to setting η = η opt from (32) to realize the result presented in Fig. 3(c,d). The opposite choice of the relative phase would instead give the typical enhanced intensity fluctuations of a phase-squeezed state. As already mentioned, a stronger τ = 0 antibunching can in principle be obtained by selecting higher momentum modes with the setup from Fig. 1b), but in Fig.3(e,f) we see that the oscillation frequency increases as well: the rapid fluctuations of g (2) (τ ) between low and high values mean that the initial antibunching signal can be easily washed away by the finite response time of a realistic detector, typically on the order or even longer than the photon lifetime in the cavity.
Another experimental difficulty may arise from the requirement of a spatially homogeneous fluid inside the microcavity. Disorder along the cavity plane may lead to unwanted scattering that could spoil the signatures of the coherent pair-creation processes, as we investigate in Appendix C.1. Additionally, the interaction with a thermal phonon bath may lead to dephasing and a reduced squeezing of the output photons; this is discussed in Appendix C.2. Finally, various sorts of uncontrolled relaxation mechanisms could lead to the building up of a thermal polariton population, see Appendix C.3.

Manipulating and probing the photon statistics with a wideaperture lens
In the previous section we have introduced a first example of an optical interference scheme that is able to convert the (weak) squeezing of in-cavity modes into a (sizeably) antibunched single-mode output beam. The proposed setup required the isolation of three k components and the subsequent manipulation by a series of linear optical elements. The present section exploits in-cavity interference between all k modes in real space to enhance the antibunching statistics of the output beam while keeping its qualitative spatial profile. The proposed configuration is depicted in Fig. 4a): A spatial image of the cavity field can be reconstructed by placing a system of two wide-aperture lenses in a confocal configuration in front of the microcavity. In the focal plane between the two lenses, a spacedependent attenuation element, e.g. based on a Spatial Light Modulator (SLM) sketched in Fig. 4b), provides the required k selection mechanism that is necessary to reduce the amplitude of the k = 0 mode with a factor F, while keeping other k components intact. A most remarkable feature of this alternative scheme is that the antibunching statistics results here from the sum of contributions of all k modes. In the following, we will illustrate how the present setup, providing the freedom to vary the coherent amplitude, offers a wide flexibility in the design of the spatio-temporal pattern of photon statistics.

Intensity correlations in position space
From ansatz (7), it is immediate to see that the real-space two-point correlation functions of fluctuations in a spatially uniform setup can be obtained from (17) as the Fourier transforms of the quantities n k (τ ) and c k (τ ) with x = |r − r |, reflecting the homogeneity of the setup. From equations (16), one sees that at large k the quadratic correlation functions behave as n k ∼ k −4 and c k ∼ k −2 , which are universal scaling laws for gases interacting with a contact potential [29]. While this poses no issues for the first-order correlation function n(x, τ ), the pair correlation c(x, τ ) suffers from an ultraviolet divergence in two or more spatial dimensions in the limit of vanishing separation x. In other words, the convenient introduction of a zeroranged contact interaction has the inconvenient consequence that pair correlations are only reproduced correctly for separations x much larger than the true range of the interaction potential. Interpolation between the two-body problem with the correct scattering potential at short distances and the many-body result at large distances has proven to be a way out of this issue [33]. In our specific case, the unavoidable finite aperture of all optical elements naturally imposes an ultraviolet cut-off in k.

The density-density correlation function
In general, the density-density correlation function at nonzero spatial separation and time delay is equal to Under the assumption of Gaussian fluctuations for the photon field, the density correlation function can be expressed by Wick theorem in terms of the two-point correlation functions and a coherent displacement field, where δn = n(0, 0) is the density of noncondensed particles and ψ f = Fψ 0 (see Expr. 7) is the attenuated condensate mode. From expression (36) one sees that the density correlation function suffers from the same ultraviolet divergence as the pair-correlations when a zero-ranged contact interaction is used. In our discussion, we therefore focus on nonzero separations x, significantly larger than the true potential range, such that our results do not suffer from this issue. In practice, we see that for each non-zero point (x, τ ), the spatial-temporal profile of g (2) (x, τ ) converges to a well-defined value for a sufficiently large cutoffand solely the point (0, 0) is suffers from the ultraviolet divergence. For this analysis, we choose the cutoff high enough such that the profiles are converged. In any practical experiment, a cutoff in momentum space is always introduced by the finite aperture of the lenses used in the imaging system.
The analysis of our numerical calculation for g (2) (x, τ ) as a function of x and τ is discussed in the next sections for the physically most remarkable cases. We will show the results for a polariton system with dimensionless interaction constant mg = 10 −4 , but we emphasize that our results stand regardless of the value of the interaction constant g. A lower g merely requires a stronger attenuation F of the k = 0 mode, provided the mean-field interaction energy µ = gn 0 , with n 0 the density of particles in the condensate, remains invariant.

Lightcone-like correlations in the high-density regime
We first analyze the ∆ < 0 case, which is obtained in the optical-limiter regime under the additional condition µ > γ, |∆|. Except for the energy gap in the very low-k range, the quasiparticles that arise in this case share similar features with the familiar phononic modes in an equilibrium condensate because interactions (quantified by µ) dominate over losses (quantified by γ). We therefore expect that the elementary excitations, once created, can travel through the fluid at the speed of sound during a time roughly corresponding to the average photon lifetime τ = γ −1 . The spatial-temporal profile of the density-density fluctuations is then expected to exhibit features that we can relate to a lightcone-like propagation of the quasiparticles. In Fig. 4d) we show the profiles of a fluid pumped with µ = 5γ and ∆ = −γ for different attenuation levels F of the k ≈ 0 modes, which give different amplitudes of the displacement fields ψ f .
It is illuminating to first discuss the case with F = 0, i.e the situation with a completely attenuated coherent field. In that case we are only looking at the photon statistics that originate from the quantum fluctuations, without interference of the coherent field. As the quasiparticles are created in pairs with opposite momenta, we expect to observe bunched photon statistics in the spatial-temporal profile of g (2) (x, τ ). We see in Fig. 4d that for F = 0 (upper panel) all statistics in the profile exhibits bunching, with higher values at (x, τ )-points that can be related to creation processes of quasiparticle pairs. Most notably, an oscillation with a frequency |∆| in time is seen, which is the frequency of low-momentum quasiparticles.
By varying the attenuation F of the k = 0 mode with the SLM (see Fig. 4b)) to add a coherent field ψ f to the signal, we can drastically change the appearance of the spatialtemporal profile of g (2) (x, τ ). Upon increasing F (i.e. transmitting a fraction of the condensate mode with the SLM, rather than blocking everything) we observe how the bunching of the Figure 4: (a) The confocal two-lens setup to measure density-density correlations of a fluid of light. A spatial light modulator (SLM) is placed in the common Fourier plane of the two lenses to perform the desired k-space selection. Delayed correlations between photon detections separated by a time interval τ allow to measure the g (2) intensity correlation function defined in (35). (b) The spatial profile of the SLM that is used for the shaping of the beam. All modes in the 2D plane are transmitted, except a small disk centered around k = 0, where the coherent field is situated. White corresponds to transmission, black to full blocking and gray to attenuating with a factor F, in order to transmit a coherent field ψ f = Fψ 0 . (cd) Spatial-temporal profiles of the density-density correlation function g (2) (x, τ ) for varying (top to bottom) filtering fraction F (see (36)) and different (left/right) pumping parameters. Red shades correspond to bunching and blue to antibunching. For the parameters of (c), the parametrically amplified modes give rise to a spatial pattern of alternating bunching and antibunching, which turns into complete bunching for F → 0. The quasiparticle dispersion of this mean-field configuration is plotted in Fig. 2(c). For the parameters of (d) We notice the appearance of an approximate sound cone x = τ /2c (white dashed lines) of antibunched correlations. Also here, when F → 0, the antibunching turns into bunching. The quasiparticle dispersion for this case is plotted in Fig. 2(a). quasiparticle pairs turns into antibunching in a sound-like band due to interference with the condensate mode (see the middle panel with F = 0.009). For a larger fraction of ψ f (see panel in (d) with F = 0.018) the profile of g (2) (x, τ ) remains largely the same in shape, but the variation from g (2) = 1, the statistics of a fully coherent state, reduces.
The expected antibunched statistics at short times and distances, stemming from interparticle repulsion, is apparent for sufficiently high displacement fields ψ f (including the standard case of no attenuation). At later times we see how the same-point antibunching quickly disappears on a time-scale of the order of a fraction of 1/µ and transforms into a propagating antibunching feature that travels through the fluid (see blue band) at a velocity close to 2c, with c = µ/m being the speed of sound (white dashed lines).
Taking inspiration from well-known results on quantum quenches and correlation functions in conservative systems [34,35,36,37], a (somewhat hand-waving) physical interpretation of this result is the following. Detection of a photon at the initial time and position creates a Bogoliubov quasiparticle in the photon fluid, which then travels away with the group velocity of the mode it is emitted into. In turn, the presence of this Bogoliubov quasiparticle modifies the probability of detecting a second photon at space-time points located along its world-line. Correlations are peaked on the lightcone, but the fact that they display a sizeable intensity both inside and outside the lightcones is due to the strong k-dependence of the Bogoliubov group velocity, that increases in a faster-than-linear way at large k because of the quadratic cavity dispersion, and tends to zero for small k because of the energy gap characteristic of the driven-dissipative condition.

Spatial pattern formation with diffusive modes
When the fluid is pumped on the lower-intensity branch of a bistable regime, a range of diffusive-like Bogoliubov modes is present. Provided losses still dominate over interactions (2µ < γ), the homogeneous state remains dynamically stable and quasiparticles generated by quantum fluctuations eventually decay. In the most interesting case 2µ < ∆, shown in Fig.  4c), the set of diffusive modes (indicated by a blue band) has a disk shape in k space around a nonzero momentum. Since the lifetime τ k = 1/(γ − 2Γ k ) (with Γ k = ω k ) of these modes can be substantially longer than non-amplified modes for which τ k = 1/γ, we expect that they may leave an important imprint on the photon statistics.
In Fig. 4c), we show the spatial-temporal profiles of g (2) (x, t) for a fluid pumped with mean-field parameters µ = 0.4γ and ∆ = 3γ and for varying displacement field ψ f , which we engineer again from the condensate mode with the SLM (see Fig. (a-b)). For these parameters we have a disk of diffusive-like modes centred around nonzero momentum k c ≈ 2.3 √ mγ (see the corresponding spectrum in Fig. 2c)). Due the parametric amplification of these modes, we see that a standing-wave-like pattern appears in g (2) (x, t), with a wavelength corresponding to the parametrically amplified wave vectors. Remarkably, the vanishing real part of the frequency of these modes implies a zero group velocity v g k = ∂ k ω k , so that the spatial pattern persists in time, practically without moving.
Upon varying F, we can again switch from a profile with alternating bunching and antibunching regions in space (lower panels of Fig. 4c, with F = 0.008, 0.016) to a profile with only bunching when the condensate is attenuated completely (panel with F = 0). There is an optimum at about F ≈ 0.008, which stabilizes a temporal band with minimum density correlations g (2) (x, t) at a separation of x ≈ 2(mγ) −1/2 . In all cases, we see that the spatial structure, as imprinted by the parametrically amplified modes, is well preserved in time. The temporal duration of the interesting correlations is now substantially longer, thereby facilitating measurement with realistic photo-detectors with nonzero photon collection time. On the other hand, the strong suppression of g (2) (x, t) at nonzero x cannot be used to generate strongly antibunched photon statistics, since it relates to correlations between two spatially separated points.
Another interesting feature of Fig. 4c) is the presence, at short times, of small ripples on top of the otherwise very stable space-time structure, which then quickly propagate away and disappear. As for the additional features that were visible on top of the lightcones in Fig. 4d), we attribute these features to the presence of modes with a nonzero group velocity. Just outside the parametrically amplified disk, we can even verify that the modes exhibit a diverging group velocity, as can be seen on Fig. 2(b-c).

Conclusion
In this work we have theoretically investigated the peculiar nature of quantum fluctuations displayed by the light field in nonlinear planar microcavities. The fluctuations are generated as a result of pair-creation processes of quasiparticles with nonzero momentum and we have proposed free-space optics configuration to manipulate their shape and intensity.
In the first part of this work, we have shown how an appropriately designed selection and interference scheme allows to translate the intrinsic two-mode squeezing due to optical nonlinearities into a single-mode squeezed output state with strongly antibunched photon statistics. Even though our results can be placed within the framework of the Unconventional Photon Blockade [15,16,17,18], the planar microcavity geometry differs from the usual two-cavity geometry typically considered in this literature.
In the second part we have illustrated how the spatial-temporal structure of the densitydensity correlation function can be reconstructed and shaped with free-space optics. By filtering the beam in the far-field, as to attenuate the k = 0 coherent component, the nontrivial spatial-temporal profile of the correlation function can be manipulated. Upon removing the condensate mode completely, the statistics generated by the fluctuations only leads to a strongly bunched signal. When we interfere an attenuated coherent mode with the fluctuations, we reconstruct a reinforced copy of the spatial profile of density-density correlations inside the fluid of light. In this framework, we have first discussed the high-density regime, where the correlation function exhibits an antibunching features at the coincidence point and then a lightcone-like behavior away from it. We then showed that the presence of a set of parametrically amplified modes with nonzero momentum leads to an oscillating spatial pattern that is well preserved in time.
From a wide perspective, the results of this manuscript suggest a new avenue to generate interesting quantum states of light using planar microcavity devices displaying only a weak nonlinearity.
Trento. MW acknowledge financial support from the FWO-Odysseus program. IC was funded by the EU-FET Proactive grant AQuS, Project No. 640800, and by Provincia Autonoma di Trento, partially through the project "On silicon chip quantum optics for quantum computing and secure communications (SiQuro)". AI was funded by an ERC Advanced investigator grant (POLTDES). SR acknowledges support from the ETH Zürich Postdoctoral Fellowship Program.

A The Bogoliubov transformation revisited
The Bogoliubov approximation provides an approximate description of a quantum field in terms of a coherent condensate with Gaussian quantum fluctuations on top. Although the formalism for a driven-dissipative fluid is very similar to the equilibrium case, we would like to draw the attention to a few notable differences.
For equilibrium atomic condensates the quantity ω k from (12) is generally known as the quasiparticle spectrum and it indicates the frequency at which a particular Bogoliubov modê χ k oscillates [29]. Caution must be taken when this view is generalized to out of equilibrium systems. The reason is that, in contrast with an equilibrium condensate, the bare-particle dispersion ε k of non-equilibrium ones is not necessarily a positive function of k. In a drivendissipative quantum fluid the condensate phase is set by the detuning δ, a tunable parameter in experiment, while at equilibrium it is fixed by the chemical potential µ, such that in that case the gapless phonon condition ∆ = 0 holds exactly (see (9)).

A.1 Evolution of the quasiparticle operators
Following expression (10), the evolution of the particle operatorsφ k can be formulated as We can always write such that V T are the left eigenvectors of B k with eigenvalues ±ω k . Notice that B k is in general not Hermitian and that therefore left and right eigenvectors do not necessarily coincide, nor must they form an orthogonal basis. Without any loss of generality, we can now define new quasiparticle operators Ξ

A.2 Regular modes
If we want theχ (±) k to be operators that satisfy bosonic commutation relations, they need to fulfil two conditions. First of all they must be each others Hermitian conjugate, which leads such that we can choose to write Secondly, the bosonic commutation relation [χ −k ] = 1 tells us that After evaluation and making use of (38), the parameters u k and v k are found as, At this point we have derived the standard text book definition of the Bogoliubov transformation without making any assumptions other than the Bogoliubov operators being bosonic [29].

A.3 Diffusive-like modes
Diffusive-like modes are characterized by having ε k < 0, such that ω k purely imaginary and we can write ω k = iΓ k with Γ k real. By making use of relation (38) we can now evaluate This clearly contradicts (42), meaning that either condition (40) or condition (42) cannot be satisfied. Henceforth, the Bogoliubov operatorsχ k of diffusive-like modes cannot be bosonic. We can choose another parametrization (not unique) for the transformation matrix U k from (38) where we choose the normalization s k r * k − r k s * k = i, such that With this choice of parametrization we derive after straightforward algebra from (38) that B The time evolution of the particle operators We focus on the linear part of time evolution of an excitation. In general, we find that equation (10) can be solved by introducing time-dependent operators of the form Here, η k (t) and ζ k (t) are time-dependent Bogoliubov coefficients. When ∆ is tunable, we can distinguish three different regimes that are separated in momentum space [25,32].
• ∆ < 0: ε k is positive for any k and the quasiparticle spectrum ω k has a gap given by |∆|(|∆| + 2µ). In the limiting case ∆ → 0 the gap closes and we retrieve the familiar linear spectrum of an equilibrium condensate. Particle-like (hole-like) excitations oscillate with a frequency ω k (−ω k ), but are damped by γ, reflecting the overall finite lifetime of particles. The time-dependent Bogoliubov transformation is then found as • 0 < ∆ < 2µ: A disk of modes k < √ 2m∆ appears where ε k < 0. Modes in this region have a purely imaginary frequency ω k = i|ω k | so that they are damped or amplified at a rate Γ k = |ω k |. Therefore one branch of excitations will be strongly damped in time with a lifetime 1/(γ + 2Γ k ), while excitations on the other branch may have a much longer lifetime 1/(γ − 2Γ k ) and are parametrically amplified. Modes in this region are traditionally called diffusive-like [25]. We derive that their time evolution is governed by with s k and r k given in (47). Note that, in order to have only exponentially damped diffusive modes, we must ensure that γ > 2µ.
• 2µ < ∆: Same as above, but in this case the diffusive modes are found on a ring 2m(∆ − 2µ) < k < √ 2m∆, while modes in the inner disk k < 2m(∆ − 2µ) oscillate with a real frequency ω k , like usual. For these modes the time evolution is found as with u k and v k given in (43).

C Main noise sources
In this Appendix, we discuss the influence of the dominant noise sources in the setup. First, we assume a homogeneous distribution of polaritons in the plane of the microcavity; this can be distorted in the presence of cavity disorder. Second, interaction with a phonon bath may lead to pure dephasing of polaritons, thus altering the squeezing properties of the light. Eventually, this together with other relaxation mechanisms may result in a population of thermal polaritons in the microcavity.

C.1 Disorder
We may give an estimate for the effects of disorder by considering the Fourier transform V k of a random potential, V (r) = 1 √ V k V k e ik·r , which is applied to the planar microcavity. We then find that the random potential enters into the evolution of the mean field in the rotating frame as When we restrict to evaluating the linear response of the mean field to disorder, we find that the non-uniform polariton field can be formulated as ψ(r) = ψ 0 + 1 √ V k δψ k e ik·r . After substitution in (52) and collecting terms up to linear order in δψ k and V k , we find a set of linear equations for each mode with the response matrix and k = k 2 /2m − δ + U n 0 . By solving (53) we derive the response of the density distribution to the disorder potential in the linear regime δn k = |δψ k | 2 = V 2 k n 0 2 k + γ 2 /4 ω 2 k + γ 2 /4 2 (55) with ω k given in (12). For this qualitative analysis we consider ω k ≈ k , such that δn k ≈ V 2 k n 0 /( ω 2 k + γ 2 /4 . When we compare this with the momentum distribution from coherent pair creation (16), we conclude that V 2 ∆V c g · µ/2, where V 2 is the variation of the disorder potential and ∆V c is the correlation volume of disorder. Plugging in numbers, we find that V (r) should not vary more than ∼ 40µeV over a scale of 1 µm for an interaction constant g = 10µeV · µm 2 , in line with [14], and µ = 300µeV . Probably, this is the primary challenge to implement our proposal in experiment, based on values reported in Ref. [38].

C.2 Pure dephasing
Pure dephasing arises when the polariton fluids interacts with a thermal bath of phonons, present in the material. In the Markov approximation, we find that this amounts to including jump operators of the form ∆V c Ψ † (r)Ψ(r) ≈ ∆V c n 0 + n 0 V k,k where ∆V c is the correlation volume of the phonons, which we estimate for simplicity as ∆V c = λ 2 dB , the de Broglie wavelength of phonons. Notice that this may be somehow more complicated when the full functional form of the phonon distribution is considered. Upon explicitly evaluating the Lindblad equation with dissipators of form (56) and integrating over space, we find that ∂ t n k deph = γ deph n 0 λ 2 dB (57) Therefore, the dominant effect of dephasing will be a scattering with phonons of polaritons from the condensate, thereby ending up in nonzero momentum modes. This is quantified by a (Markovian) rate γ deph . Notice that long-wavelength phonons may not allow for a simplified Markovian treatment. Still, it is expected that (56), describing a build-up of incoherent polaritons from scattering out of the condensate, will be the main influence of pure dephasing.
The finite-momentum density resulting from these spurious scattering processes builds up linearly in time and is damped with photon decay rate γ of losses from the cavity. Therefore, we obtain in the long time limit that δn k deph = γ deph γ n 0 λ 2 dB (we take that the coherent density 16 is a number of order 1). Here, n 0 λ 2 dB is the number of condensate particles within a volume λ 2 dB and can be a relatively large number for a typical residual temperature of the order of µ = gn 0 , the chemical potential. As such, if γ deph is too large, but still γ deph < γ, we propose to employ a pulsed excitation scheme to circumvent this issue. Given the complexity of the dephasing process, it is difficult to obtain an accurate estimate of γ deph or to extract it from the literature. Still, it is clear that the condition γ deph < γ/(n 0 λ 2 dB ) must be satisfied for a CW pumping scheme to be employed.

C.3 Other noise sources
As a consequence of other spurious relaxation processes, an incoherent population of polaritons can build up at the bottom of the excitation branch. This may result in an extra density of polaritons n inc at nonzero k, not generated by coherent pair-creation, which is to be added to n k in (16) and therefore reduces the squeezing of the output light. One way to circumvent this problem is also to employ a pulsed excitation scheme. Then, the polariton population is expected to build up on the time scale of the polariton lifetime, while thermalization into the bottom of the lower polariton branch requires some relaxation process, which typically occurs on a much longer time scale.