Wigner negativity in the steady-state output of a Kerr parametric oscillator

The output field from a continuously driven linear parametric oscillator may exhibit considerably more squeezing than the intracavity field. Inspired by this fact, we explore the nonclassical features of the steady-state output field of a driven nonlinear Kerr parametric oscillator using a temporal wave packet mode description. Utilizing a new numerical method, we have access to the density matrix of arbitrary wave packet modes. Remarkably, we find that even though the steady-state cavity field is always characterized by a positive Wigner function, the output may exhibit Wigner negativity, depending on the properties of the selected mode.


I. INTRODUCTION
As opposed to a quantum field confined in a cavity, a field propagating in free space is characterized by a continuum of modes. In order to give a Schrödinger picture description of the propagating output state, wave packet modes are introduced. In quantum optics experiments, it is possible to infer information about a state stabilized inside a cavity by looking at the output field in a wave packet mode with a bandwidth matching the inverse cavity decay time. The only requirement imposed on this wave packet function is that it is square-integrable and normalized to unity, which guarantees that the filtered signal corresponds to a single bosonic mode. Whereas for a linear nondriven system, there exists a unique wave packet mode which maps the cavity state into the output field [1], this is not the case in the presence of a continuous drive. In this case, nothing restricts us from going beyond the cavity natural bandwidth with arbitrary function profiles.
Depending on the selected wave packet, properties of the cavity field and its corresponding output field can be quite different. An old and well-known example of this difference is given by the degenerate parametric oscillator (PO). This is commonly employed to generate squeezed states, i.e., states exhibiting quadrature fluctuations below the vacuum level. In the steady state, the maximum attainable squeezing inside the PO cavity is equal to 50% of the vacuum noise. Nonetheless, the output field squeezing can largely surpass that, and ideally achieve 100% squeezing in a narrow frequency bandwidth around the cavity frequency [2].
Squeezed states are a prominent example of nonclassical states of light [3][4][5]. A signature of the nonclassical behavior of a squeezed state is oscillations in the photon-number distribution [6][7][8]. Such oscillations are not by themselves a signature of a nonclassical state, since classical states can also display population oscillations. A useful nonclassicality criterion which studies the departure of a photon-number distribution from a classical probability distribution was introduced by Klyshko [9]. More commonly, the nonclassical nature of a quantum state is defined in terms of so-called quasiprobability distributions [10]. Examples of these are the Husimi Q function, the Wigner function, and the Glauber-Sudarshan P function. A generally accepted definition of nonclassicality is given in terms of the P function: if it is negative or more singular than a Dirac δ function the state is considered nonclassical [11,12]. However, the singular nature of a nonclassical P function makes it difficult to access experimentally. Another widely used nonclassicality criterion is negativity of the Wigner function [13,14], which in contrast can be directly measured in experiments [15]. In addition, negativity of the Wigner function plays an essential role in quantum computation with continuous variables, as it is a distinguishing feature of states that are resourceful for a computational advantage [16,17]. In this work we will focus on the Wigner function as an indicator of nonclassical behavior.
The process that generates squeezing in a PO is (degenerate) parametric down-conversion, in which a pump photon of frequency 2ω splits into two photons each with frequency ω. This process occurs in a medium with a second-order nonlinear susceptibility [18]. Real-life nonlinear materials are however not restricted to second order, and higher-order nonlinearities may also need to be taken into account. A higher-order nonlinearity is also a requirement for the generation of Wigner-negative states. Including a third-order (Kerr) nonlinearity results in the Kerr parametric oscillator (KPO) which has recently found applications in microwave quantum optics for the dissipative stabilization [19][20][21] and the adiabatic preparation [22][23][24] of quantum states of light (cat states), which are useful for quantum computing [25,26]. Additionally, by utilizing superconducting circuits it is possible to explore the so-called Kerr single-photon regime in which the strength of the nonlinearity surpasses the inverse cavity decay time by an order of magnitude, enabling the observation of previously undetected quantum effects [23,[27][28][29].
In this work we study the steady-state output from the KPO cavity with respect to temporal wave packet modes. We find that the KPO cavity and output fields can have markedly different nonclassicality properties, mainly in terms of the Wigner function. For comparison, we also review the PO output field in terms of temporal wave packet modes. While the squeezing spectrum of the output is previously known, we analytically solve for the full output field density matrix, shedding new light on this well-studied system by investigating how its features depend on the properties of the temporal wave packet. Since an analytical solution is not straightforward for the KPO output field, we use numerical simulations using the "input-output with quantum pulses" formalism introduced by Kiilerich and Mølmer [30,31]. We establish the nonclassical character of the KPO output field through its photon-number distribution using the Klyshko nonclassicality criterion. We find that the presence of the Kerr nonlinearity leads to larger population oscillations than for the PO. We also find that under circumstances that render the steady-state KPO cavity field Wigner-positive, the leaked output field may nonetheless exhibit Wigner negativity, and the magnitude of a Klyshko coefficient directly correlates with the amount of Wigner negativity. Moreover, by numerical optimization we find the temporal wave packet that maximizes the Wigner negativity.
The paper is structured as follows: In Sec. II we describe our system model. In order to make a comparison with the KPO, we review the PO output field in Sec. III before proceeding with the KPO output in Sec. IV. There, we compare the population statistics between the cavity steady state and the output state in order to understand how Wigner negativity can appear in the output despite the cavity being Wigner-positive. In Sec. V we take a closer look at the nonclassical properties of the KPO output field for different nonlinearity strengths. Finally, Sec. VI summarizes our results.

II. THE MODEL
Our system is a parametrically driven nonlinear cavity, and we are interested in its steady-state emission. In a frame rotating at the cavity resonance frequency ω, the Hamiltonian of the driven cavity is (h = 1) where β is the complex two-photon drive (parametric pump) amplitude, K is the Kerr nonlinearity, andĉ (ĉ † ) is the annihilation (creation) operator for the photons inside the cavity. This is a ubiquitous quantum optics model which describes squeezing and parametric amplification [18,32]. The twophoton drive results from the nonlinear interaction of two modes, typically denoted as the signal and the pump. Using the the so-called parametric approximation in which the pump field is assumed to be classical, the pump mode operator can be substituted by a c-number. As an alternative to having a second-order nonlinear crystal in an optical cavity, this Hamiltonian can also be obtained via a time-dependent boundary condition such as a movable mirror [33] or a tunable Josephson inductance in a microwave circuit [34][35][36].
The nonunitary dynamics of the cavity state is well described by the Lindblad quantum master equation. This equation follows from the interaction of a system (in our case, the cavity) and its environment under the weak-coupling approximation and assuming that the environment is memoryless. These approximations are often collectively referred to as the Born-Markov approximation [37]. Considering the environment to be in a thermal state at zero temperature, the master equation is where {·, ·} denotes the anticommutator and γ is the singlephoton loss rate of the cavity. In typical quantum optics applications, single-photon loss is the main decay channel of cavity photons into the the continuum of electromagnetic modes which comprise the cavity environment. Here, we consider photon emission into a one-dimensional transmission line. The state of the fieldâ out that leaks out of the cavity can be inferred from the input-output boundary condition [38], which relatesâ out with the drive field and the cavity field operatorĉ. Considering the two-photon drive as a weakly coupled input channel, the cavity output field depends solely on the cavity state asâ

Wave packet modes
The cavity output field provides information of the cavity state through the input-output relation (3). However, the cavity field corresponds to a single bosonic mode, while the output is comprised by a continuum of frequencies. Therefore, in order to do a faithful comparison of the cavity and output states, it is necessary to define a bosonic mode out of this continuum. This is done in terms of wave packet modes. For convenience, we are going to restrict to a temporal description in the remainder of this work, but the corresponding frequency representation is simply related by a Fourier transform.
We are going to characterize the emission from the cavity in the Fock space defined by the wave packet creation operator with f (t ) satisfying the normalization condition ∞ 0 dt | f (t )| 2 = 1 in order forÂ f to fulfill the bosonic commutation relation [Â f ,Â † f ] = 1. For simplicity, we restrict f (t ) to be a real-valued function. This creation operator defines a symmetric wave packet in frequency space around the cavity resonance frequency.
Experimentally, detection in a given temporal profile is implemented by means of a pulsed local oscillator in homodyne detection or by processing the measurement signal by means of a digital filter corresponding to the function f . For this reason, we sometimes refer to the wave packet profile as a filter function.

III. PARAMETRIC OSCILLATOR
Before moving to the Kerr parametric oscillator, we start by studying the linear, or Gaussian, degenerate parametric oscillator (PO) with K = 0 in Eq. (1). It is known that the steady state of the field inside the resonator exhibits quadrature squeezing. This means that the minimum value of the variance of the generalized quadrature operatorX θ = (be −iθ + b † e iθ )/ √ 2 is smaller than the standard quantum limit set by vacuum noise. Hereb is a generic bosonic field operator [b,b † ] = 1; i.e., it can equally refer to the cavity operatorĉ or the filtered outputÂ f . For our chosen normalization of the generalized quadratures, the vacuum noise variance is 1/2. We label the minimum value of the variance as and consequently, squeezing is indicated by s < 1/2. The variance is defined as and it attains its minimum value for a particular value of θ . For our setup, it is dependent on the phase of the drive β in Eq. (1). But for simplicity, from now on we are going to restrict ourselves to a real-valued two-photon drive amplitude β ∈ R without loss of generality. In our particular system we have for the PO steady state ĉ ss = 0. Consequently, first moments of the quadrature field also vanish. In general, it is always possible to displace the field in order to meet this condition, so from now on we will ignore first-order moments. Then, the variance can be rewritten as ( X θ ) 2 = b †b + b 2 e −2iθ /2 + b †2 e 2iθ /2 + 1/2. Taking into account that b †2 = b 2 * and noting thatb †b is a positive-semidefinite operator, we see that there is squeezing if and only if Re[ b 2 e −2iθ ] < − b †b . We see that squeezing can only appear if the expectation value b 2 has a large enough magnitude. This expectation value is associated with off-diagonal terms in the density matrix: coherences between number states |n and |n − 2 .
The amount of squeezing of the cavity field increases as we approach the so-called pump threshold β th = γ /2, where the system becomes unstable, i.e., the mean number of photons diverge. However, the maximum attainable squeezing never surpasses 50% of the vacuum noise, i.e., s = 1/4. On the other hand, the field leaking out of the cavity can achieve perfect squeezing s → 0 near threshold in a very small bandwidth around the cavity frequency [2]. Below we will throw light on this phenomenon in a new way by looking at the filtered output in the Schrödinger picture, demonstrate that the squeezing increases with increased filter time, and provide an intuitive explanation for this. Since an infinite filtering time corresponds to a zero-bandwidth frequency filter at the cavity frequency, our results are consistent with the previous results.

A. PO output state
With K = 0, the system is Gaussian since the equation of motion (2) is only quadratic inĉ andĉ † . Hence the output field is completely characterized by its first-and second-order moments for which it is possible to find a closed set of equations. This means we can analytically solve for the state of the output field. For the selection of a wave packet mode, we choose a constant filter within a time interval T , a so-called boxcar filter [39][40][41]. This simplifies the analytical calculations.
As already mentioned, we will consider the steady-state emission from the resonator. Using the quantum regression theorem [42][43][44] we calculate the steady-state two-time correlations ĉ † (τ )ĉ(0) ss and ĉ(τ )ĉ(0) ss by assuming the steady-state condition ∂ t ss = 0 in Eq. (2). Correlations for the filtered output stateÂ f follow from integration: and From these correlators we can calculate the covariance matrix which completely determine the state of a Gaussian system. With these matrix elements, the Wigner function can be calculated as (9) and the Fock space representation of the state is determined following Refs. [45,46] (more details in Appendix A). The behavior of the filtered output field of the parametric oscillator as a function of the filtering time T for two different drive strengths β = 0.2 and β = 0.4 is shown in Fig. 1. Both T and β are given in units where the decay rate is set to γ = 1. In Figs. 1(a) and 1(b) we plot the first six Fock state populations ρ n . For both drive strengths, the single-photon state is the first state to be populated and is the dominant nonvacuum constituent of the complete state for T 1. Progressively, two-, three-, and higher photon number states become populated. At around T 2, the two-photon population overcomes the single-photon population. Increasing the filtering time, we observe the same dynamics for the other even photon number states (the four-photon population overcomes the three-photon one, and so on). As we go toward T → ∞ we see that the odd Fock state populations are suppressed and the state becomes a superposition of even Fock states. From the analytical solution (not shown), we can see that as we approach β th in the limit T → ∞ the populations of all the even states become identical. This tendency can be seen by comparing the asymptotic populations in Figs. 1(a) and1(b): the differences between the populations are smaller for the stronger drive.
In Figs. 1(c) and 1(d) we plot the squeezing defined by Eq. (5) as a function of T . Starting as vacuum noise for T 1, the output field becomes squeezed as soon as we start to generate photon pairs. It reaches the squeezing level of the intracavity field for T 1, i.e., around the natural bandwidth defined by the cavity decay rate, and then surpasses it for larger T . The maximum squeezing is attained when the state becomes a superposition of even Fock states.
A squeezed state which satisfies the minimum uncertainty allowed by the Heisenberg uncertainty principle is sometimes referred to as an ideal squeezed state [47]. To show that the PO output field becomes an ideal squeezed state in the limit T → ∞ we plot the determinant of the covariance matrix in Figs. 1(e) and 1(f). Because it is a real symmetric matrix, the covariance matrix can always be diagonalized. This diagonal matrix V contains the variances of the squeezed quadrature and the orthogonal one. Their product, i.e., the determinant of V , is the uncertainty product. The minimum uncertainty condition in our quadrature normalization corresponds to det V = 0.25. For values of T much smaller than the cavity characteristic decay time 1/γ , the minimum uncertainty condition is met because the probability of emission is very small, and thus, we are witnessing the minimum uncertainty of the vacuum state. As the squeezed states approach the minimum uncertainty condition for large T the populations of the odd Fock states go asymptotically to zero. In Appendix B we show analytically that the odd photon numbers vanish only for a squeezed state which is also a minimum uncertainty state.
Finally, in Fig. 2 we show the first 20 populations of the output field for two different widths of the boxcar filter, T = 4 and T = 500, for β = 0.2 as well as β = 0.4. We note that vacuum is always the leading contribution to the state. For T = 4 we can appreciate an enhanced two-photon population for both drive strengths [Figs. 2(a) and 2(b)], although for the stronger drive in Fig. 2(b) many higher number states are already also populated. For both drive strengths, as T becomes very large, in Figs. 2(c) and 2(d) we only observe even Fock states in agreement with Fig. 1. This distinct even-odd oscillation in the populations is a sign of nonclassicality, which will be elaborated further on in Sec. III B. The insets in Fig. 2 show the Wigner functions of the corresponding states, which become more and more squeezed for larger β and T .
These results can be intuitively understood by means of a very simple time-domain argument as follows: the two-photon drive places correlated pairs of photons inside the resonator. However, the photons leave the cavity one by one. This restricts the number of photon pairs inside the resonator and, therefore, the maximum amount of squeezing that can be sustained. On the other hand, we can access the photon pairs by monitoring the output field for a long time compared to the cavity lifetime 1/γ , as illustrated in Fig. 3. Thus, large two-photon correlations between states |n and |n − 2 can be obtained, and as described in the beginning of this section, therefore more squeezing is expected in the output.

B. Nonclassicality from the PO
There exist several indicators and measures of the nonclassical behavior of light, including sub-Poissonian statistics [48], antibunching [49], as well as the previously mentioned photon-number population oscillations and squeezing. As discussed in the introduction, it is also very common to define nonclassicality by the behavior of quantum phase space quasiprobability distributions such as the P function or Wigner function. Negativity of the Wigner function is a sufficient but not a necessary condition for nonclassicality. For instance, the Wigner functions of the PO cavity and output fields are always positive. This is a consequence of the states resulting from a linear, or Gaussian, system [50]. Nevertheless, quadrature squeezing can be related to a nonclassical P function [48,51]. Since the squeezed state exhibits popu- lation oscillations as the parametric drive creates excitations in pairs, a perhaps more useful [52] nonclassicality criterion for squeezed-like states is given by the Klyshko inequality [9] B n ≡ (n + 1)ρ n−1 ρ n+1 − nρ 2 n < 0.
The coefficient B n is defined in terms of the populations of three consecutive number states |n − 1 , |n , and |n + 1 (n 1). The Klyshko inequality sets a bound for the population ρ n in terms of those of its nearest neighbors. If B n < 0 for any n, the photon-number distribution of the corresponding field departs from a classical probability distribution; i.e., its P function is negative [53]. The nonclassical nature of states produced in parametric down-conversion has been established via the Klyshko inequality in the optical regime through direct use of photon counters [54]. Similarly, in the microwave regime the nonclassical nature of propagating squeezed states generated through a Josephson parametric amplifier was established by feeding the field into a 3D cavity and reading its population content via a dispersively coupled transmon qubit [55].
The odd Klyshko coefficients are always positive for our system, so in Fig. 4 we show the first three even Klyshko coefficients B 2 , B 4 , and B 6 for the filtered output of the PO as a function of the width T of the boxcar filter function, using two different drive strengths β = 0.2 and β = 0.4. It can be seen that for T > 1 the state is manifestly nonclassical as Eq. (10) is first satisfied for n = 2. As we have already seen in Figs. 1 and 2, for T 2, the two-photon state becomes the dominant nonvacuum contribution to the output field, and the inequality (10) for n = 2 implies that Therefore, when the population of the two-photon state overcomes this bound imposed by the populations of the singleand three-photon states the output field becomes nonclassical. The dominance of ρ 2 becomes stronger for larger values of T as the population of the odd-number states start to decrease. As T is increased, higher even-number states also become populated, which explains the successively appearing negative values of B 4 and B 6 .
Furthermore, we observe a qualitative difference between the behavior of the Klyshko coefficients for the two drive strengths. While the coefficient B 2 exhibits the largest negativity in both cases, for β = 0.2 it attains its minimum value in the limit T → ∞, when the odd-number states are fully suppressed. On the other hand, for β = 0.4, as we increase T we involve considerably more (even) number states, each with smaller populations [cf. Figs. 2(c) and 2(d)] as T → ∞. The Klyshko B 2 coefficient peaks along with the peak of the two-photon population, which is confirmed to occur just before T = 10 in Fig. 1(b).

IV. KERR PARAMETRIC OSCILLATOR (KPO)
The steady state of the cavity field for the case of a nonvanishing Kerr nonlinearity has also been exhaustively studied in the literature. In fact, despite being a nonlinear system, the steady state admits an analytical solution. It is a well-established result that the steady state of our system is characterized by a positive Wigner function. Specifically, with successively increased drive strength, the cavity steady-state field transitions from vacuum to a weakly squeezed state and finally into an incoherent superposition of coherent states, with quantum coherence washed out by interaction with the environment [29,[56][57][58][59].
While there has been comprehensive studies of the cavity field, to the best of our knowledge, a thorough characterization of the properties of the filtered output field is still missing in the literature. In the remainder of this paper we are going to study the filtered output field of the KPO, with special emphasis on its nonclassical properties using the above introduced Klyshko nonclassicality criterion as well as negativity of the Wigner function.
Characterizing the state of the propagating cavity output field is, in general, a difficult problem since it implies the calculation of multitime field correlations for different time orderings. For a Gaussian or linear system like the PO studied in Sec. III, the output field is completely characterized by its first-and second-order moments for which it is possible to find a closed set of equations. In the presence of the nonlinearity, this is no longer possible. Instead, we would get an infinite number of equations involving every possible order of the output field moments. For this reason we do not attempt to characterize the output field by calculating multitime correlations, but instead follow a different approach. In order to explore the features of the output field of the KPO we are going to make use of a technique recently introduced by Kiilerich and Mølmer [30,31]. We implemented it using QuTiP [60], and the code can be found in Ref. [61]. The numerical solutions were validated against the analytical solution for the PO. Alternatively, one could rely on stochastic methods to mimic a quantum tomography experiment, as done in Refs. [40,41], or numerically solve the full Schrödinger equation for the coupled cavity and output fields [24].
Before discussing the nonclassicality of the KPO output in depth, we will introduce the Poissonian regime of the KPO, in which both the cavity and the output are classical states with Poissonian photon-number statistics. We then discuss the twophoton population of the KPO output, which will be shown to have a strong connection to its nonclassical features.

A. KPO cavity Poissonian regime
The nonlinearity prevents the instability at the drive threshold β = β th . This means that the mean number of photons in the cavity n cav no longer diverges at this point, but grows steadily with β. When the photon number reaches n cav γ /8K, the cavity steady state is an incoherent superposition of coherent states [22]. Here the photon-number distribution is Poissonian, ρ n = n n cav exp(−n cav )/n!. Therefore, we will refer to this type of field configuration as the Poissonian regime.
The maximum of a Poissonian photon-number distribution is centered around the average number of photons present in the field. So while the average number of photons in the cavity grows along with the parametric drive strength β, when the Poissonian regime is reached, individual number state populations reach their peak and start decreasing in succession as the peak of the Poissonian distribution shifts to higher and higher number states. As such, we expect that the largest achievable population for the nth number state happens when the cavity field population has a Poissonian distribution with n cav = n. Then, it is possible to estimate the largest two-photon population in the KPO cavity, which is ρ * 2 0.27. In the next section we study how the KPO output field two-photon population depends on the filtering time and the nonlinearity strength. Similarly to the PO, we can get a two-photon population that dominates over the one-and threephoton populations for filtering times comparable to or longer than the cavity decay time.

B. Filtered KPO output field: Boundary of the Poissonian regime
Intuitively, the output two-photon population is expected to grow with β until the system reaches the Poissonian regime, similarly to the cavity state. Here we are going to show that, in contrast to the cavity field, for the output field it is possible to produce states which exhibit the same maximum two-photon state population ρ * 2 but with non-Poissonian photon number statistics. This is significant, because in Appendix C we show that when the cavity is in the Poissonian regime, so is the output field. Specifically, the output field is also a classical mixture of coherent states, with the average number of photons N f given by the cavity photon number n cav and a scale factor depending on the width of the filter function.
For a boxcar filter the mean number of photons in the output field is given by N f = γ n cav T . Thus, in the Poissonian regime, the largest two-photon population is achieved when N f = 2. Then, from the photon-number scaling relation, we can calculate the time T * at which N f = 2. This corresponds to T * = 2/(γ n cav ) [62]. Therefore, the largest two-photon population expected in the output of the KPO in the Poissonian regime is ρ * 2 0.27 for a filtering time T * . We need to be outside of the Poissonian regime to observe nonclassical effects. Below, we show the behavior of the photon populations in the filtered output field as we depart from this regime for different values of the nonlinearity strength. As stated in Sec. IV A, not being in the Poissonian regime implies a limited drive strength β, which in turn requires an increased filtering time T in order to collect a significant number of photon pairs. In Fig. 5(a) we show the output two-photon population as a function of β for a weak nonlinearity K = 0.1 and different values of T (in units of γ = 1). As expected, we observe that the two-photon population peaks at weaker drive strengths for larger values of T , and vice versa. Additionally, we note that as T is increased, the photon population statistics change. For the shortest filtering time T = 0.1, the drive strength β at which the maximum two-photon population is attained corresponds to the cavity field being in the Poissonian regime, and subsequently, the output is also Poissonian, as can be seen from the photon-number distribution in Fig. 5(b). Here the largest two-photon population observed corresponds to ρ * 2 0.27 (dotted line) and the filtering time at which it is achieved agrees with T * . For a slightly increased T it is possible to generate output states with a Poissonian-like photon-number distribution [cf. Fig. 5(c)] but with a smaller two-photon content than ρ * 2 , as for T = 0.4. The photon-number distribution of the output state slowly departs from the Poissonian behavior as we further increase T , as exemplified with T = 1 in Fig. 5(d). But crucially, the two-photon population is reduced more and more as T is increased and we move farther out of the Poissonian regime.
Interestingly, we do not need a strong nonlinearity to get dramatically different behavior. In Fig. 6 we show results for a moderate nonlinearity K = 0.5. Here, as we decrease the drive strength β and increase the filtering time T , we depart from the Poissonian regime similarly as with the weak nonlinearity. But in contrast, even when not in the Poissonian regime, the largest two-photon populations for K = 0.5 still correspond approximately to ρ * 2 . In fact, from our numerical simulations the largest two-photon output can even slightly exceed this value, as can be seen in Fig. 6(a).
For the longer filtering time in Fig. 6(d), population oscillations start to become evident. But it can be noted that opposed to the PO, for which we know that the odd-number states can be completely suppressed for a long filtering time as shown in Sec. III B, it is not possible with a nonzero Kerr nonlinearity (cf. Appendix B).
In the next section we are going to study the nonclassical features of the filtered output field of the KPO. As described in Sec. III B, nonclassicality is not only a consequence of the enhanced two-photon population, but by how large it is compared to its nearest neighbors, i.e., single-and three-photon states. For example, the photon distributions in Figs. 6(b), 6(c), and 6(d) have roughly the same two-photon populations, yet only Fig. 6(b) corresponds to a Poissonian distribution. In Fig. 6(c), the state begins to depart from the Poissonian regime and the two-photon population overcomes that of its nearest neighbors. Further away from the Poissonian regime this asymmetry is even larger, as in Fig. 6(d). We are going to quantify this using the Klyshko coefficient B 2 .

V. NONCLASSICALITY IN THE KPO OUTPUT
In this section, we will evaluate the nonclassicality of the KPO output field in terms of the Klyshko B 2 coefficient and Wigner negativity. Starting with the Klyshko criterion, we display the minimum B 2 for different nonlinearities K in Fig. 7. The Klyshko inequality B 2 < 0 is satisfied which establishes When the nonlinearity K is increased from zero, the Klyshko B 2 coefficient rapidly drops to its minimum for K = 0.7, corresponding to the "most nonclassical" state as determined by the Klyshko criterion. As K is further increased, the value of B 2 increases slightly before saturating. The minimum values were obtained by sweeping over a grid of β and T for each K.
the nonclassical nature of the KPO output field. Specifically, when K grows toward 0.5 and larger, the enhanced two-photon contribution leads to a very sharp population oscillation [an example of this is Fig. 6(d)] and, consequently, to a very large magnitude of the Klyshko coefficient. The largest magnitude is obtained for K = 0.7. A stronger nonlinearity does not increase the nonclassicality; in fact, it decreases slightly before it saturates for K 2. Further increasing K does not lead to fundamentally different behavior, since the cavity steady-state field is entirely determined by the ratio β/K in the strong nonlinearity regime [29,56,57,59]. For a weak, or even vanishing, nonlinearity the value of B 2 can remain <0, but the magnitude is severely diminished compared to what is attainable with a moderate or large K.
The rather large two-photon populations in the KPO output field may also result in more quantum coherence, namely, in larger density matrix elements ρ 02 and ρ 20 . This is because the magnitude of these matrix elements is bounded by the populations of the vacuum and two-photon states: |ρ 02 | √ ρ 0 ρ 2 (the equality is only achieved for a pure state). Quantum coherence typically translates into negative regions in the Wigner function. Unfortunately there is no simple relation between the ρ 02 coherence and Wigner negativity, as higherorder Fock states heavily influence the negativity.
In Fig. 8 we show the Wigner function and density matrix for the K = 0.5, β = 0.65 steady-state cavity field, and compare this with the output field for T = 2.5 and T = 5.0 boxcar filters. Surprisingly, even though the cavity steady-state field is always Wigner-positive [57,63], the KPO output field can obtain a negative Wigner function for certain values of the filter width T . Negativity in the output increases as a function of T until T = 5.0 where it peaks, and goes back down if T is further increased. Interestingly, the peak occurs at T = 5 not only for β = 0.65, but for all β 0.3. The ρ 02 coherences are similar for the two output fields, and they are only slightly larger than for the cavity field. But the main contribution to the Negativity appears in the Wigner function, and in the density matrix, the vacuum population is reduced and the two-photon population is increased compared to the cavity field. (e) Wigner function and (f) density matrix for the output field with a T = 5 boxcar filter, which gives the largest WLN (0.05). Here, the vacuum population is further reduced and the two-photon population is even more prominent. For clarity in the figure we truncate the density matrices at n = 3, but the full states occupy a larger Fock space. cavity-state coherence is from the vacuum population, which inhibits Wigner negativity.
As the source of quantum coherence is the two-photon drive and the output state is favored toward even-number states, it is no surprise that the KPO output Wigner function resembles a two-component kitten or cat state. Nevertheless, the resulting nonclassical states have a very low purity, roughly 0.617 for the most Wigner-negative state observed.
A possible measure of the negativity of the Wigner function is the integrated Wigner negativity [14], or alternatively the more recently introduced Wigner logarithmic negativity (WLN) [64]: It has the property W > 0 when the Wigner function W (x, p) has a negative part, and has the benefit of being an additive resource monotone [65]. The resource monotone can be defined using any logarithm base, but for our numerical calculations we use the natural logarithm. In Fig. 9(a) we show a map of the WLN as a function of both the two-photon drive strength β and the boxcar width T for K = 0.5. As can be seen, the maximum value of the negativity is achieved near T = 5. This holds for K 0.3. In the strong nonlinearity regime the maximum negativity occurs for β/K 1.
It is interesting to compare the behavior of the Klyshko coefficient B 2 and the WLN. This is shown in Fig. 9(b). As it can be seen, the nonclassical population oscillations around the two-photon state which result in negative values of B 2 directly translate into nonclassicality of the Wigner function. We show this correspondence for T = 5, but it holds for every value of T .
These results establish that whereas the cavity decay rate imposes a detection bandwidth for which the cavity state can be output with high fidelity, detection with a wave packet beyond this natural limit may reveal a completely different nature of the output field. In addition, while single-photon losses may destroy quantum coherence inside the cavity, it does not represent a loss mechanism for the output field. In the next section, we will briefly study the effects of the temporal profile of the wave packet.

Impact of the filter function on the Wigner negativity
So far, we have for convenience selected a temporal mode of the cavity output field with a boxcar filter, but any squareintegrable function can define a bosonic mode in accordance with Eq. (4). Since the filter response can significantly change the nature of the detected state [66], it is reasonable to suspect that the choice of filter function can have an impact on the observed Wigner negativity. In this section, we are going to target the temporal mode profile which gives the maximum WLN by means of numerical optimization.
Since there is literally an infinite number of possible temporal modes, the best filter could easily be left out if a selection of filter functions were tested manually. To ensure that the best filter function was found, even if it was not a well-behaved, smooth function, we performed an optimization of the numerical array that represents the filter using the "scipy.optimize" package [67]. The part of the filter that is zero until steady state has been reached was fixed and not included in the optimization. Besides the first and final points being zero, the initial filter was random (but properly normalized, i.e., i f 2 i = 1). These constraints on the first and final points as well as normalization were enforced during the optimization. Due to the randomness of the initial filter, the results of different optimization runs were not identical. But in general, the optimized filter obtained a Gaussian shape. A representative example is displayed in Fig. 10. In fact, the maximum WLN obtained with a Gaussian filter is twice the maximum of the boxcar. A comparison between the two is shown in Fig. 11 for T = 5 which gives the maximum WLN for the boxcar filter, and σ = 2.3 which gives the maximum WLN for the Gaussian filter.

VI. SUMMARY AND CONCLUSIONS
In this paper we have studied the nonclassicality of the output field of the steady-state Kerr parametric oscillator, in terms of the Wigner function and Klyshko coefficients. Our main result is that whereas the KPO cavity is Wigner-positive in the steady state, the output field can be Wigner-negative, depending on the properties of the selected field mode. In order to obtain the state of the output field we have defined bosonic modes in terms of wave packet functions, utilizing the new "input-output with quantum pulses" formalism introduced by Kiilerich and Mølmer [30], which allows us to obtain the density matrix of the output field.
We also revisited the well-studied linear parametric oscillator. The linear PO is instrumental for the generation of quadrature-squeezed states of light. It is also a paradigmatic example of the different properties exhibited by the cavity and output fields of a continuously driven setup. While the results for the KPO were obtained by numerical simulations, for this linear system we could study the properties of the filtered output field by reconstructing its density matrix analytically from two-time output field correlations. Here we also explored the so-called even-odd population oscillations as a function of the temporal width of the wave packet function. To the best of our knowledge, these oscillations have previously only been studied by direct photon detection, which is insensitive to the mode structure of the field.
We could then contrast the output of the PO to that of the KPO. The nonlinear Kerr parametric oscillator is also ubiquitous in quantum optics literature. But even though its cavity steady state has been analytically solved, to the best of our knowledge the output field has not been studied beyond its squeezing properties. We found that the presence of the nonlinearity leads to stronger population oscillations, which is expressed by the larger magnitude of the Klyshko coefficient. This is what gives rise to the Wigner negativity in the KPO output, as the magnitude of the Klyshko coefficient directly correlates to the integrated Wigner logarithmic negativity (WLN). Furthermore, by numerical optimization we have verified that the nonclassical properties of the output field are dependent on the chosen wave packet function, and that a Gaussian wave packet maximizes the WLN.
Typically, the important parameter responsible for driven nonlinear oscillators to reach quantum regimes is the ratio between the Kerr parameter and cavity decay rate; i.e., efficiency of quantum nonlinear effects requires a high nonlinearity with respect to dissipation [68]. In contrast, there is no need for a strong nonlinearity to observe Wigner negativity in the KPO steady-state output, as K/γ = 0.7 is sufficient to observe the largest Klyshko negativity and the largest Wigner negativity.
Our results could realistically be verified experimentally. Notably, in superconducting circuit experiments the nonlinearity strength can be tuned [69]. Quantum-state tomography of propagating microwave fields has been established using both phase-insensitive [1,70] and phase-sensitive amplification [71], and more recently via parity detection [72]. In addition, the Klyshko nonclassicality criterion is amenable to be tested for propagating fields following a recent proposal for a microwave number-resolved photon counter [73]. Finally, our approach could be extended to study the output field properties of recent higher-order squeezing realizations [74,75]. Also, it would be interesting to study the role of filtering in the output entanglement properties of multimode setups.

APPENDIX A: DENSITY MATRIX FROM THE COVARIANCE MATRIX FOR A GAUSSIAN STATE
A Gaussian state is defined by its covariance matrix V with matrix elements withx andp the position and momentum quadratures, respectively,x withb (b † ) the bosonic annihilation (creation) operation which might refer to a cavity or filtered propagating mode. Here we are assuming that x = p = 0 which is true for the models studied in the main text. Recall that in the steady state of (2) with Hamiltonian (1) we have ĉ ss = 0 and, consequently, Â f = 0 for the filtered output field. The density matrix elements in the Fock or number basis can be recovered from the covariance matrix V by means of the relation with d = det V and t = tr V the determinant and the trace of the covariance matrix, respectively [45]. The so-called multidimensional Hermite polynomials H R mn (x 1 , x 2 ) are defined in terms of a 2 × 2 matrix R with elements [46] The arguments x 1 and x 2 are related to the first moments of the field, which in our case are always zero. We get with H n (x) being the nth-order Hermite polynomial.

APPENDIX B: SUPPRESSION OF ODD FOCK STATE POPULATIONS
The Heisenberg uncertainty principle puts a lower bound on the minimum value of d = det V , where V is the covariance matrix of a quantum state, defined by Eqs. (A7)-(A9). Under the quadrature normalization used here we have d 1/4. A minimum uncertainty state is by definition a state for which d = 1/4. Examples of these are coherent states and the so-called ideal squeezed states [47]. In a squeezed state, noise (the variance) in one quadrature is reduced below the vacuum level at the expense of increased noise in the orthogonal quadrature. For an ideal squeezed state, the product of the quadrature variances equals the lower bound.
A quintessential example of an ideal squeezed state is the squeezed vacuum state, that is, the state that results from the action of the unitary squeezing operator with ξ = r exp(iθ ) (r, θ ∈ R) on the photon vacuum state |0 . For θ = 0, the variances satisfy V 11 = x 2 = 1 2 exp(−2r) ( B 2 ) and V 22 = p 2 = 1 2 exp(+2r), and, thus, d = 1/4 [47]. The suppression of the odd Fock state populations is a consequence of ideal squeezing. Indeed, the condition d = 1/4 results in R 12 = 0 in Eq. (A9). If this is the case, only the term k = 0 will contribute to the summation in Eq. (A10). For m = n the latter reduces to [H m (0)/m!] 2 . All of the odd-order Hermite polynomials are identical to zero at the origin (x = 0). Consequently, n|ρ|n = 0 in Eq. (A6) for odd n.

APPENDIX C: TIME FILTERING OF A COHERENT STATE
Let us assume that the steady state of the cavity field is a coherent state. A coherent state is completely defined by the first moment of the bosonic field operator, i.e., ĉ ss with higher-order moments factorizing in terms of it. Following the input-output relation, we have for the filtered output field moments where for simplicity we are assuming the filter function f to be real. Using the coherent field property that higher-order moments factorize, we get with similar relations holding for third-and higher-order correlations. For a boxcar and Gaussian wave packets, the above time integral depends on the width and the variance (T and σ ) of these functions, respectively. That is, the effect of the filter is just to rescale the cavity moments proportionally to the temporal width of the filter. Finally, for an incoherent superposition of coherent states the above arguments hold for every coherent state in the mixture.