Entanglement and Photon Anti-Bunching in Coupled Non-Degenerate Parametric Oscillators

We analytically and numerically show that the Hillery-Zubairy’s entanglement criterion is satisfied both below and above the threshold of coupled non-degenerate optical parametric oscillators (NOPOs) with strong nonlinear gain saturation and dissipative linear coupling. We investigated two cases: for large pump mode dissipation, below-threshold entanglement is possible only when the parametric interaction has an enough detuning among the signal, idler, and pump photon modes. On the other hand, for a large dissipative coupling, below-threshold entanglement is possible even when there is no detuning in the parametric interaction. In both cases, a non-Gaussian state entanglement criterion is satisfied even at the threshold. Recent progress in nano-photonic devices might make it possible to experimentally demonstrate this phase transition in a coherent XY machine with quantum correlations.

The fundamental quantum resources of CIMs come from the phase-sensitive amplification/deamplification of two quadrature amplitudes in the DOPO [22]. Quantum correlation among DOPOs, evaluated by entanglement and quantum discord, reaches a maximum at the DOPO threshold and decreases below and above the threshold [23]. This quantum resource comes at the cost of limiting the spin-like degrees of freedom for computation and simulation. The in-phase quadrature amplitude (a canonical coordinate of a harmonic oscillator) takes either one of the bi-stable values above the threshold. We can relax this binary constraint using phase-insensitive oscillators such as lasers [24][25][26] or exciton-polariton condensates [27,28], in which the optical or polaritonic field can take a continuous phase above the threshold, so that an XY Hamiltonian (or Kuramoto model) can be naturally implemented instead of the Ising Hamiltonian. We call such a network as a coherent XY machine (CXM) in this paper. However, the bosonic fields in lasers and polariton condensates are thermal statistical mixture states below the threshold, while they approach coherent states far above the threshold [29]. The lack of non-classical states in lasers and polariton condensates seems to exclude the possibility of creating quantum correlations in a CXM based on such classical oscillators.
Non-degenerate optical parametric oscillators (NOPOs) have recently been used to construct a CXM [30]. One of the experimental advantages of an NOPO-based CXM is that we can easily transform a CIM to a CXM by introducing the frequency non-degeneracy between the signal and idler waves without changing the basic structure of the machine. An NOPO is a phase insensitive oscillator with a continuous phase degree of freedom, but its quantum statistical features are unique and distinct from standard lasers and polariton condensates. An NOPO below the threshold has a stronger gain saturation than a standard laser, so that a sub-Poissonian light or even a single photon state may be generated if the system parameters are chosen appropriately. This non-classical behavior below the threshold is analogous to those of a strongly coupled atom-cavity system [31,32] and coherently excited Raman three-level system with a large coupling and detuning [33]. On the other hand, an NOPO above the threshold can produce an amplitude squeezed state with a reduced amplitude fluctuation due to the strong gain saturation and diverging phase fluctuation due to a random walk diffusion process. This non-classical behavior above the threshold is analogous to those of pump-noise-suppressed lasers [34][35][36][37][38][39]. When such non-classical states of light are mixed by dissipative (Liouvillian) coupling, it is expected that the quantum correlations will form among the two NOPOs, just as in the case of two DOPOs in a CIM [3,4,23]. This would be another advantage of NOPO-based CXMs.
In this paper, we investigate the formation of entanglement in a CXM consisting of two NOPOs with a large parametric gain and show that entanglement is achieved below, above, and at the threshold. We use both analytical and numerical methods, although the analytical one is used to calculate the entanglement characteristics only far below and far above the threshold. Our analytical study on above-threshold CXMs, started with c-number stochastic differential equations (cSDEs) in the positive-P representation [40,41]. The numerical simulations were performed using both the quantum master equation (QME) in the photon number representation, and the wave function Monte Carlo (WFMC) method [42]. The paper is organized as follows. We introduce the model of CXM consisting of two NOPOs in Section 2. In Section 3, we investigate the case of a large dissipation in the pump mode and find that a large detuning of the parametric interaction is required to satisfy the entanglement criterion below the threshold. Next, in Section 4, we consider the case of a large dissipative coupling. We find that, below the threshold, entanglement is obtained even in the absence of detuning in the parametric interaction. Section 5 summarizes the paper. Appendix A shows the Fokker-Planck equation of the positive-P quasi-distribution function and derives the below-threshold characteristics of a single NOPO. Appendixs B and C provide subsidiary information on the numerical and theoretical equations. Appendixs D and E provide the detailed derivations of analytical results. Appendixs F and G show the supplementary information about the entanglement criterion and the experimental NOPO [30].

Model
The density matrix equation of a single χ (2) NOPO is written as Here,â p ,â s , andâ i are pump, signal, and idler modes, respectively. These three photon modes have cavity frequencies ω a (a = p, s, i) and the linear dissipation rates γ a (a = p, s, i). ∆ = (ω p − ω s ) − ω i is the detuning of the parametric interaction. ε is the strength of the coherent excitation of the pump mode. κ is the strength of the χ (2) parametric interaction. If the idler mode has a much larger dissipation than the other two photon modes (γ i γ p , γ s ), the quantum master equation can be written as Here G is the coefficient of the effective Raman interaction in Ref. [43][44][45], and K is the coefficient of the non-degenerate Kerr effect [46]. If ∆ = 0, the Kerr coefficient is K = 0 and the signal mode has a parametric gain G = κ 2 /γ i . We define the maximum parametric gain as G 0 := κ 2 /γ i . For the detuned parametric interaction, the parametric gain is G = κ 2 γ i γ 2 i +∆ 2 , and the Kerr coefficient is We use the dimensionless detuning parameter d = K 2 /G 2 to represent the normalized detuning. We also introduce the normalized excitation p = ε/ε thr , where ε thr := γ p γ s /G is the strength of the excitation at the oscillation threshold.
Let us consider a CXM consisting of two NOPOs (NOPO1 withâ p1 ,â s1 and NOPO2 withâ p2 ,â s2 ) coupled by a dissipative Liouvillian L cρ for two signal modes [3,4,[47][48][49], We will evaluate the entanglement between the two signal modes using one of the entanglement criteria in Ref. [50], If this value is larger than zero, the two signal modes are entangled. Above the threshold, we assume the ferromagnetic configuration â s1 = â s2 is created by the coupling Liouvillian (4).
Using the fluctuation of positive-P amplitudes (introduced in Appendixes A and B), the normalized HZ1 is written as HZ1 Here, we assumed ∆α † s1 ∆α s1 = ∆α † s2 ∆α s2 and that â s1 is real. Introducing the canonical coordinates and canonical momenta as (i = 1, 2), this normalized entanglement criterion is written as This criterion is equivalent to the Duan's sufficient condition for entanglement (Theorem 1 of Ref. [51]). Moreover, HZ1 criterion can detect an entanglement in a non-Gaussian state, for example, entanglement of a superposition of single photon states |ψ = 1 √ 2 (|0, 1 + |1, 0 ) [52]. This criterion seems useful to detect non-Gaussian entanglement in the belowthreshold CXM, since the similar model [33] realizes a single photon state.

Quantum Master Equation after the Elimination of the Pump Mode
In this section, we consider the large pump mode dissipation limit, where G, K, γ p γ s , J. In this limit, the linear loss of the pump mode, due to the dissipation into the reservoir and spontaneous emission into the signal mode, is much larger than the linear loss of the signal mode. We will use the expansion of the density matrix with the complex-P representation [53,54] for the pump mode and photon number state representation for the signal mode: Substituting this expansion into Equation (2), we can obtain a time development equation Due to the linear dissipation of the signal mode (γ s ), components with photon number indices (N s , N s ) are excited by components with larger photon number indices (N s + 1, N s + 1). The parametric gain (Gα † p α p ) introduces a contribution from components with smaller photon number indices (N s − 1, N s − 1). Other processes do not affect the photon number indices. Although the equation has drift terms for the pump amplitudes (α p , α † p ), the signal photon number indices do not change when components are derived using α p or α † p . When γ p + G is sufficiently large, the complex-P amplitude (α p ) at photon number indices (N s , N s ) rapidly converge to the steady-state value of the time-development equation We can eliminate the pump mode amplitudes by writing , and integrating the complex-P amplitudes by dα p dα † p . The density matrix components of the signal mode ρ N s ,N s = N s |ρ|N s are, Here G e (N, N ) = . The denominators of these terms have higher dependence on the signal photon number than in Scully-Lamb's theory [55]. In regard to the CXM consisting of two NOPOs, we can omit the subscript for the signal mode s, after eliminating the pump mode. The density matrix components of the two signal modes ρ N 1 ,N 2 ,N 1 ,N 2 = N 1 , N 2 |ρ|N 1 , N 2 develop as,

Far-Below-Threshold Entanglement
From the above Equation (10), we will analytically derive the photon anti-bunching and entanglement characteristics far below the threshold (p 1). In this limit, G e (N i , N i ) and K e (N i , N i )(i = 1, 2) are of order O(p 2 ). Therefore, these contributions to ρ N 1 ,N 2 ,N 1 ,N 2 on the right hand side of Equation (10) are much smaller than those of γ s and J. However, the G e (N i − 1, N i − 1)(i = 1, 2) on the right hand side are not negligible. Although the K e (N i , N i )(i = 1, 2) on the right hand side do not contribute for small p, the non-degenerate Kerr coefficient K in the denominators of G e (N i − 1, N i − 1)(i = 1, 2) contributes to the characteristics far below the threshold. The signal photon number of the CXM is obtained from the equations of the two density-matrix components, ρ 10,10 , and ρ 10,01 . ∂ρ 10,10 ∂t = −2(γ s + J)ρ 10,10 + 2Jρ 10,01 + 2G e (0, 0)ρ 00,00 , ∂ρ 10,01 ∂t = −2(γ s + J)ρ 10,01 + 2Jρ 10,10 .

Far-Above-Threshold Entanglement
Next, we analytically derive the above-threshold characteristics of the CXM by assuming a large pump dissipation γ p γ s and small parametric gain G γ s . We derive an analytical expression for the above-threshold fluctuations from the equations of positive-P representation [40,41]. For a single NOPO, the Fokker-Planck equation of the positive-P Here, α p , α † p (α s , α † s ) are positive-P amplitudes for the pump (signal) mode. Using the Ito rule, the equivalent c-number stochastic differential equations are expressed as: Here ξ C , ξ † C and ξ ‡ C are independent complex Gaussian random numbers satisfying correlation functions, such as ξ * C reflects random spontaneous emission of signal photons. ξ C and ξ † C contribute to the correlation between the pump and signal amplitudes. From Equations (23) and (24), the above-threshold pump photon number is identical to that of an NOPO without a Kerr term K [39]: The above-threshold signal photon number is obtained from γ s /G ∼ | α p | 2 and which is the steady-state mean of Equation (21), Here, we have assumed that G γ p and α p α † s α s ∼ α p α † s α s . Using d = K 2 /G 2 and ε = γ p p γ s /G, the steady-state signal photon number is obtained as Here, π = p 2 (1 + d) − d is the normalized excitation modified by the detuning parameter d = K 2 /G 2 in the parametric interaction. The steady-state pump amplitude is written as By introducing the phase factor the mean pump amplitude becomes α p = γ s G e −iφ . Now let us introduce the small amplitude fluctuations in the pump and signal modes, ∆α p and ∆α s . The pump amplitude fluctuation is defined as the fluctuation after removing the phase factor denoted by φ.
The mean signal amplitude rotates with the frequency K α † p α p ∼ √ dγ s due to Kerrnonlinearity. The signal amplitude fluctuation is defined as the fluctuation after removing this time-dependent phase factor, These phase factors do not affect the products, α † p α p and α † s α s , appearing in the positive-P equations. The equations for small amplitude fluctuations are Here, r = γ p /γ s and ξ ‡ The time t is normalized so that 1/γ s is the unit time. The equivalent equations written with canonical coordinates and momenta are shown in Appendix C.
From these equations in the limit of r 1, we can calculate the steady-state value of Mandel's Q parameter [56] for a signal mode. It is defined as Q M,s := ∆n 2 s − n s n s , wherê n s =â † sâs and ∆n s =n s − n s . When Q M,s is smaller than zero, the signal mode is in an amplitude squeezed state. From the calculation shown in the Appendix D, is obtained. In the large excitation limit (π → ∞), Q M,s converges to a negative value − 4+j 4(2+j) . Therefore, an amplitude squeezed state is obtained even in a small-G NOPO. In the single NOPO limit (j → 0), Q M,s = − 1 2 + 1+d 2(π−1) − d 2π . This value converges to Q M,s → − 1 2 for a large pump excitation, which is identical to the value obtained in Ref. [35].
. Therefore, for the same p, detuning in the parametric interaction increases the amplitude squeezing. In the strong dissipative coupling limit, i.e., j → ∞, the Mandel's Q parameter is halved from that of the single NOPO limit.
The steady-state value of HZ1 entanglement criterion is also obtained for large r, as The detailed derivation is provided in Appendix D. Here, in the large excitation limit (π → ∞), HZ1 converges as HZ1/ â † sâs → j 2 −4(1+d) 4j(j+2) . The coupling coefficient j must satisfy j > 2 √ 1 + d for entanglement far above the threshold. Therefore, detuning d is not preferred for satisfying above-threshold entanglement criterion, although d > 1 must be satisfied for below-threshold entanglement.
The analytical and numerical Mandel's Q parameter and HZ1 entanglement criterion are compared in Figure 1c,d. Both results are obtained for J/γ s = 12 and d = 5 and plotted as the function of normalized excitation p. The analytical results assume that γ p /γ s → ∞ and G/γ s → 0. The numerical results were obtained by integrating the positive-P c-number SDEs in Appendix B, for γ p /γ s = 50 and G/γ s = 10 −7 . To obtain steady-state statistics from the positive-P calculation, we calculated a single trajectory, with excitation depending on time as p(t) = p min(1, tγ s /10 5 ), and took the time average from tγ s = 10 5 to tγ s = 10 6 . The theoretical Mandel's Q M,s parameter from Equation (36) becomes negative at p ∼ 1.7. The numerical Q M,s values were slightly larger than the theoretical values due to the finite γ p /γ s ratio. When d = 0, Equation (36) always gave Q M,s = 0 at p = 2. Therefore, the nondegenerate Kerr effect causes amplitude squeezing for smaller p. Moreover, the theoretical HZ1 from Equation (37) becomes positive for p ∼ 2.1. The numerical values are slightly smaller than the analytical results due to the finite γ p /γ s ratio. Above the threshold, the use of a small γ p /γ s ratio helps to decrease the degrees of both amplitude squeezing and entanglement. This differs from the below threshold case, where using small γ p /γ s ratio improves the anti-bunching and entanglement.

Numerical Results
Here, we give the numerical simulation results as a function of normalized excitation p. The mean signal photon number â † sâs and HZ1 entanglement criterion are presented in Figure 2a,b for γ p /γ s = 50, G/γ s = 400, d = 5, and J/γ s = 12. For a small excitation p, we numerically calculated the quantum master equation (QME) in Equation (3)  . This means that the calculation time increases rapidly for large p, where the mean signal photon number â † sâs increases and we have to prepare a sufficiently larger M s satisfying M s â † sâs . The maximum value of p calculated by QME was p = 10, where the mean signal photon number is â † sâs ∼ 0.42, and we have used M p = 4 and M s = 7. The method of the QME calculation is the same as for Figure 1a,b.
For a larger excitation p, we used Mølmer's wave-function Monte Carlo (WFMC) calculation [42], whose calculation time is O(M 2 p M 2 s ). In the WFMC calculation, time development started from the vacuum state with the excitation p(t) = p min(1, tγ s /10 2 ). A time average was taken from tγ s = 10 2 to tγ s = 10 3 − 2.5 × 10 4 , depending on the values of normalized excitation p. The calculation for small p with â † sâs < 1 required more samples (longer period for time averaging) due to the slow convergence, although the size of the Hilbert space is smaller than that of the case of large p. The maximum p value calculated by WFMC was p ∼ 178, where the mean signal photon number was â † sâs ∼ 10, where we set M p = 4 and M s = 36. The blue dashed line shows the results of the density matrix calculation with the pump mode eliminated in Equation (10). The calculation time of this method is O(M 4 s ), so we could use it to calculate even the abovethreshold characteristics. We calculated the characteristics at p ∼ 316 with M s = 46, where â † sâs ∼ 17. The numerical simulation indicated that the mean signal photon number does not show a rapid increase at the threshold. This behavior is known as a thresholdless lasing [57] and is obtained for a large nonlinear saturation coefficient. For an NOPO, different from a Scully-Lamb laser [55], the dependence of the signal photon number on p becomes smaller above the threshold. The below-threshold signal photon number is proportional to p 2 (Equation (13)), and the above-threshold signal photon number is proportional to p (Equation (27)). As a consequence of this behavior, a small stimulated emission below the threshold helps to reduce the signal photon number below the O(p 2 ) line. The Hillery-Zubairy's entanglement criterion is satisfied far-below and far-above the threshold. This is an expected from the analytical results (K is much larger than G below the threshold, while above the threshold j > 2 √ 1 + d is satisfied). However, the small stimulated emission below the threshold gives a negative correction to HZ1/ â † sâs a minimum value near the threshold (p ∼ 19), where the signal photon number â † sâs exceeds one. Far above the threshold HZ1/ â † sâs 2 converges to the theoretical values (purple dash line) obtained from dividing Equation (37) by Equation (27). Figure 2c,d presents the numerical simulation results for larger coupling coefficient J/γ s = 120 (other parameters are the same as in Figure 2a,b). Equation (20) leads us to expected that a larger coupling coefficient increases the normalized HZ1 value far below the threshold. Although a stimulated emission contributes negatively to the entanglement criterion, in a similar way to Figure 2b, HZ1 always has a positive value even at the threshold. We thus obtained entanglement in below, above and at the threshold of the CXM, starting from a large pump mode dissipation.

Large Dissipative Coupling Limit
In the previous section, we showed that achieving entanglement at the threshold is not straightforward. The correction of the small stimulated emission to normalized HZ1 is negative, and the entanglement criterion at the threshold can only be satisfied by increasing J/γ s to the same order of G/γ s . In Figure 1a,b, it can be seen that a large linewidth ratio (γ p /γ s ) reduces the anti-bunching and entanglement. The theory for a large pump mode dissipation assumes the coupling coefficient J is much smaller than γ p , G, K, which keeps J small and prevents quantum correlation between the two signal modes. Here, we consider another limit where J is sufficiently larger than the other parameters (γ p , γ s , G, K).

Far-Below-Threshold Entanglement
First, we will consider the quantum statistics far below the threshold. Appendix A describes the procedure for deriving quantum statistics for a single NOPO far below the threshold. Briefly, this procedure starts from the truncated Fokker-Planck equation (Equation (A3)) in the positive-P representation, after the stimulated emission and Kerr non-linearity terms removed (which we write as ∂P ∂t | NOPO,ε→0 ). We integrate the Fokker-Planck equation to obtain the time development equations of mean amplitude products and use them to derive the steady-state photon number and correlation functions. For a CXM, we start from the following Fokker-Planck equation, Some of the steady-state results obtained from this Fokker-Planck equation are similar to, or even identical to the results for a single NOPO in Appendix A. Even with the large coupling of a signal mode, the pump amplitudes satisfy α †i p α j p = ( γ p p γ p +G γ s G ) i+j , which is identical to those of a single NOPO. However, as J → ∞, the signal photon number becomes one half of that of a single NOPO α † s1 α s1 ∼ γ 2 p p 2 2(γ p +G) 2 . This is the same as the J → ∞ limit of Equation (13) obtained for large pump mode dissipation, and is also identical to the same limit of the two signal amplitudes' correlation α † s1 α s2 in Equation (14). We assume that such amplitude products connected by J, i.e., α †i p1 α j p1 α †k s1 α l s1 , to any replacement of α † s1 → α † s2 and α s1 → α s2 , have the same value in the zeroth order of 1/J. The details of the derivation of g (2) s (0) are shown in Appendix E. In the large-J limit, the below-threshold second-order correlation function g (2) s (0) can be written as, In the limit γ s → 0, this value converges to g (2) s (0) = 8(γ p +G) 2 (2γ p +3G) 2 +K 2 . This is different from the J → ∞ limit of Equation (19). If we start from large dissipative coupling J → ∞, the anti-bunching state g (2) s (0) → 0 is obtained with large non-degenerate Kerr coefficient K. From Equation (39), the below-threshold entanglement criterion obeys as HZ1/ â † sâs 2 = 1 − g (2) s (0). The normalized HZ1 with no detuning in the parametric interaction (K = 0) is This equation shows that the far-below-threshold entanglement criterion is satisfied even with K = 0. When γ p = γ s , it is satisfied for G/γ s > 15.6. Figure 3a compares the normalized entanglement criterion HZ1/ â † sâs 2 obtained from the quantum master Equation (3) and the theory in the J → ∞ limit (1 − g s (0) with Equation (39)). We set the linewidth ratio γ p /γ s to 4 and the maximum parametric gain G 0 /γ s = κ 2 /(γ i γ s ) to 40 and investigate the impact of normalized detuning d in the parametric interaction. With detuning, the parametric gain becomes G = G 0 1+d and the non-degenerate The numerical methods are the same as in Figure 1b. The numerical results with the largest j (j = 3240) have almost the same values as those in the theory assuming J → ∞. Even in the J → ∞ limit, normalized HZ1 slightly increases for a small detuning d. As d increases, K decreases as ∝ 1 √ d and gets closer to γ s , where below-threshold entanglement is impossible. When J/γ s = 360, entanglement criterion is satisfied for d = 0. However, when J/γ s = 120, the entanglement criterion is not satisfied with zero detuning, although HZ1 becomes positive with a non-zero detuning parameter.

Far-Above-Threshold Entanglement
Next, we derive the far-above-threshold entanglement in the limit J → ∞. First of all, we will consider a single NOPO. From the cSDEs in Equations (32) and (33) These pump amplitude fluctuations are excited by two fluctuation correlations, ∆X p ∆X s , and ∆P p ∆X s , where ∆X s is the amplitude fluctuation of the signal mode. We introduce the normalized amplitude correlations, The steady-state pump amplitude fluctuations in terms of A and B are Next, we obtain the time development equations of ∆X p ∆X s , ∆P p ∆X s and ∆X 2 s .
B = 1 is obtained from the steady state of Equation (51). Substituting Equations (46) and (48) into Equations (49) and (50), we obtain the steady-state signal amplitude fluctuation, Here, This fluctuation theory for a single NOPO is sufficient for calculating the normalized HZ1 in a CXM in the j 1 limit. HZ1 above the threshold obeys Equation (6) and includes a contributions from canonical momenta. However, for canonical momenta of signal modes, the fluctuation products of the CXM obey The steady-state difference of these two fluctuation values is which converges to zero in the large j limit. The sum ∆X 2 s1 + ∆X s1 ∆X s2 appearing in the HZ1 is the same as ∆X 2 s of a single NOPO shown in Equation (52). Finally, the normalized HZ1 in the J → ∞ limit is Mandel's Q parameter for the signal mode Q M,s = 2 ∆X 2 s1 is the same as −HZ1/ â † sâs ∼ ∆X 2 s1 + ∆X s1 ∆X s2 in the J → ∞ limit. If we take the r → ∞ limit in Equation (57), the normalized entanglement criterion converges to HZ1 â † sâs = 1 4 − π+d 4π(π−1) , which is the same as the j → ∞ limit of the theory with a large pump dissipation (Equation (37)). This independence of the order of taking limits does not apply below the threshold. The O(1/r) correction to HZ1/ â † sâs is negative: − 1 rπ − 2d(π−1) rπ(π 2 +d) . Therefore, above-threshold entanglement is more easily obtained with large r (see Figure 1d).
We performed a numerical simulation using a positive-P equations in Appendix B. From Figure 3a with γ p /γ s = 4, we chose two sets of parameters: (d = 0, j = 360) and (d = 1, j = 120). Below-threshold entanglement was achieved with these parameters, when we set a large parametric gain G 0 /γ s = 40. Here, to check the far-above threshold entanglement, we performed the positive-P calculation for G/γ s = 10 −7 , with the same method as Figure 1d. The results are shown in Figure 3b. The numerical results fit the analytical ones (Equation (57)), which assume the J → ∞ limit. For d = 1, j = 120, the numerical HZ1 values are slightly smaller than the theoretical values, because the correction of O(1/j) reduces HZ1, as shown in Equation (37). Nevertheless, the values of normalized HZ1 for d = 1 are larger than for d = 0, because, as discussed in relation to Figure 1c, the O(d) correction makes it easier for an above-threshold NOPO to have a non-classical state.

Numerical Results
Here, we present the numerical simulation results as a function of normalized excitation p, for γ p /γ s = 4, G 0 /γ s = 40, d = 0, and j = 360. The mean signal photon number and normalized HZ1 are shown in Figure 4a,b. The numerical results are obtained by solving the QME (3) for small-p, and by performing a WFMC calculation [42] for large-p. The Methods are the same as those used to calculate the results shown in Figure 2, but for WFMC, the time average was taken from tγ s = 10 2 to tγ s = 10 4 − 10 5 depending on the excitation p. For the smallest p (∼5.6), time average was taken from tγ s = 10 2 to tγ s = 10 5 because of the slow convergence. For the maximum p ∼ 56 for WFMC, where the mean signal photon number is â † sâs ∼ 5.0, we used a photon number space with M p = 5 and M s = 21 and took the time average from tγ s = 10 2 to tγ s = 10 4 . The purple dashed lines far-above the threshold are plots of Equation (27), or Equation (57) divided by Equation (27). The numerical results converge to the theoretical values far-above the threshold. The entanglement criterion is satisfied below and above the threshold, although the nonlinear Kerr effect is absent (K = 0). In contrast to Figure 2b,d, a small stimulated emission gave a positive correction to the normalized entanglement criterion HZ1/ â † sâs 2 . Because of this correction, the entanglement criterion is satisfied, even at the threshold. Next, we present the numerical simulation results for γ p /γ s = 4, G 0 /γ s = 40, d = 1, j = 120. As shown in Figure 4c,d, for the maximum p value (∼32) the mean signal photon number was â † sâs ∼ 4.3. We used M p = 6 and M s = 21 for calculating this value. As with Figure 4b, stimulated emission gives a positive correction to the normalized HZ1 value, which also results in entanglement at the threshold. As expected from Figure 3b, the peak value of the normalized HZ1 is slightly larger than in Figure 4b.

Summary
We showed that Hillery-Zubairy's entanglement criterion is satisfied in coherent XY machines, below, above, and even at the threshold of a CXM consisting of two highly non-Gaussian χ (2) -NOPOs. We investigated two limits: (1) the pump mode has much larger dissipation than the signal mode, and (2) the dissipative coupling coefficient is much larger than the other parameters. In the first limit, below-threshold entanglement is possible only when the parametric coupling is detuned. In the second limit, below-threshold entanglement was obtained even when the parametric coupling is not detuned. In the first limit, although detuning of the parametric interaction is necessary to achieve below-threshold entanglement, it prevents above-threshold entanglement if J is comparable to γ s . Moreover, the normalized entanglement criterion HZ1/ â † sâs 2 is decreased by a small stimulated emission, while in the second limit, the same value increased. The experimentally required G/γ s for at-threshold entanglement is smaller for the second limit than the first. From these considerations, the second case with large dissipative coupling seems to be a more effective way to achieve entanglement at the threshold. For a more detailed study, other entanglement criteria should be discussed (In Appendix F, we discussed the Simon's necessary and sufficient criterion for Gaussian state entanglement [58]). For obtaining entanglement at the threshold, the small stimulated emission correction to below-threshold HZ1/ â † sâs 2 , would be important. By further increasing G and K from the values we used in this paper, quantum state production in quantum spin model [59] would be possible in a CXM.
A large parametric interaction, κ/γ s ∼ 10 −2 , has been experimentally confirmed in a second-order nonlinear (χ (2) ) rib-waveguide-based microring resonator [60]. The theoretical model studied in this paper could be realized by stimulated Raman scattering of traveling-wave modes in a silicon rib-waveguide [61], or standing-wave modes in silicon photonic crystal nanocavities [62,63], coupled via low-Q cavity mode [64]. For traveling wave model, the pulse period must be longer than the phonon lifetime to avoid unintentional correlation. In the standing-wave model, the parametric gain normalized by linear dissipation is calculated as G γ s ∼¯h c 2 g R Q cs ε r V cav ∼ 9 × 10 −10 Q cs [62]. Here, c is the speed of light in vaccum, Q cs is the quality factor of the signal cavity mode, g R is the coefficient of Raman amplification, ε r is the relative dielectric coefficient of the material, and V cav is the volume of the cavity (we used g R ∼ 57 cm/GW [65], ε r = 12 and V cav ∼ 0.5 µm 3 [62] for silicon photonic crystal nanocavity). The cavity quality factor must be as large as Q cs ∼ 10 10 to achieve below-threshold entanglement. This is only an order of magnitude larger than the numerically achieved Q factor [66]. The hyper-parametric oscillation which enabled the CXM with third-order (χ (3) ) nonlinearity [30] can have similar characteristics in terms of below-threshold anti-bunching or above-threshold amplitude squeezing (Appendix G). The CXMs of χ (3) NOPOs seem to realize entanglement due to these non-classical characteristics.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Fokker-Planck Equation and Far-Below Threshold Characteristics of Single NOPO
Here, we derive the below-threshold characteristics of a single NOPO using the positive-P Fokker-Planck equation. The positive-P distribution function P(α, α † ) is defined for a bosonic modeâ as [40,41], For a single NOPO obeying Equation (2), the motion of the positive-P amplitudes, for the pump mode (α p , α † p ) and for the signal mode (α s , α † s ), depends on the Fokker-Planck equation of the quasi-distribution function P(α p , α † p , α s , α † s ).

∂P ∂t
Next, we derive the far below-threshold characteristics (ε → 0). Assuming small signal photon number (α † s α s 1), we remove terms representing stimulated emission from the pump mode to the signal mode, and terms representing Kerr nonlinearity. After these terms have been removed, the Fokker-Planck equation of From this truncated Fokker-Planck equation, the mean of pump amplitude α p (= α p Pd 2 α p d 2 α † p d 2 α s d 2 α † s ) develops as The second-order correlation function of the signal mode of a below-threshold single NOPO is written as g (2) s (0) = 2 γ s + (γ p + G)Re γ p +2γ s +G γ p +2γ s +2G+iK This value is 2 for small G, K [29], and it converges to zero in the large-K limit [33]. Even when K = 0 (no detuning in the parametric interaction), a single NOPO can have an anti-bunching state (g (2) s (0) < 1) for large G. However, the minimum value is 0.5, instead of zero as obtained from the detuned parametric interaction.

Appendix D. Details of Above-Threshold Theory in the Large Pump Mode Dissipation Limit
In this appendix, we derive the above-threshold theory of CXM in the large pump mode dissipation limit. Assuming large r( 1), we can eliminate the pump mode and obtain the equations of signal-mode amplitude fluctuations from Equations (32)-(35), as Here, ξ R and ξ † R are independent real Gaussian noise satisfying ξ R (t)ξ R (t ) = δ(t − t ) and ξ † R (t)ξ † R (t ) = δ(t − t ). Using these equations, we will consider an CXM consisting of two NOPOs. Defining the coefficients, A : π(π−1) p 2 (1+d) , B := 2π p 2 , and D := − 2(π−1)(π+i , the fluctuations of the canonical coordinates obey, where j = J/γ s . Assuming a Gaussian state, the Mandel's Q parameter is obtained as (36)). Next, we calculate the above-threshold values of Hillery-Zubairy's entanglement criterion following Equation (6). The first two terms related to ∆X s are obtained as The mean fluctuation of ∆P s diverges above the threshold of the NOPO. We can use a small lasing linewidth parameter δ 1 to produce a finite fluctuation, The HZ1 contains the difference of these values, which is independent of small δ, Next, we consider the mean products of the ∆X s and ∆P s amplitudes, appearing in Equation (A30): Since ∆X 2 s1 − ∆X s1 ∆X s2 = ReD+B 4(ReA+j) , the steady-state difference of these values is Finally, the normalized HZ1 criterion is obtained as Equation (37).

Appendix E. Details of Below-Threshold Theory in the Large Dissipative Coupling Limit
In this Appendix, we derive the second order correlation function g (2) s (0) for belowthreshold CXM in the large dissipative coupling limit, using the assumption that values coupled by J have the same value in the large J limit. α †2 s1 α 2 s1 and the three values coupled to it via J obey, There is no coupling coefficient J in the time-development equation of . Assuming that this weighted average is identical to α †2 s1 α 2 s1 in the zeroth order of 1/J, we obtain From α † s1 α s1 = G 2γ s α † p1 α p1 , the steady-state value of the second-order correlation function can be written as g is the correlation function between the pump photon number and signal photon number.
Next, we calculate the correlation g X between the pump and signal photon numbers from the three following equations: Moreover, there is no coupling coefficient J in the equation of . Assuming that this equation is identical to α † p1 α p1 α † s1 α s1 in the zeroth order of inverse coupling coefficient 1/J, the mean photon number product is, From this equation, the steady-state pump-signal photon number correlation can be written as, is the correlation between the pump amplitude and the signal photon number.
To obtain the steady-state correlation function g X between the pump amplitude and the signal photon number, we derived the time development equations of four mean amplitude products, The average of these four mean amplitude products obeys the equation without the coupling coefficient J. If we assume this average is identical to α p1 α † s1 α s1 in the zeroth order of 1/J, we obtain From this equation, the steady-state correlation g X between the pump amplitude and the signal photon number is The second order correlation function g (2) s (0) of the far below-threshold large-J CXM is obtained in the form of Equation (39).

Appendix F. Simon's Criterion for Above-Threshold Entanglement
For entanglement above the threshold, we can use Simon's criterion [58], which is necessary and sufficient condition for entanglement when the states are Gaussian. When the above-threshold states of CXMs are treated as the Gaussian states with infinite fluctuations of canonical momenta, we can use this criterion. This criterion is described by a covariance matrix of two signal modes σ ij = R i R j − R i R j + δ ij , where − → R = √ 2[X s1 , P s1 , X s2 , P s2 ] (X si and P si (i = 1, 2) are defined for positive-P amplitudes). If the minimum eigenvalue of iΩ(ΛσΛ T ) (which we callν − ) is smaller than 1, for a partial-transposition matrix Λ = diag(1, 1, 1, −1), and Ω = 0 1 −1 0 ⊕ 0 1 −1 0 , the two signal modes are entangled [67]. From the equations in Appendix C (for d∆P s /dt, we add −δ∆P s to the right hand side to avoid divergence), we obtain the steady-state values of 20 fluctuation products, ∆X 2 p1 , ∆P 2 p1 , ∆X 2 s1 , ∆P 2 s1 , ∆X p1 ∆P p1 , ∆X p1 ∆X s1 , ∆X p1 ∆P s1 , ∆P p1 ∆X s1 , ∆P p1 ∆P s1 , ∆X s1 ∆P s1 , and these between different NOPOs. From six values ∆X 2 s1 , ∆P 2 s1 , ∆X s1 ∆P s1 , ∆X s1 ∆X s2 , ∆P s1 ∆P s2 , and ∆X s1 ∆P s2 , we calculated the minimum symplectic eigenvalueν − of a partially transposed covariance matrix. First, we study the limit of large pump mode dissipation. For γ p /γ s = 10 3 , d = 5, j = 12, and δ = 10 −10 , we plot theν − − 1 as a function of excitation p. The comparison with HZ1 criterion is shown in Figure A1a. The entanglement criterion is satisfied for smaller p(∼1.99) than the HZ1 criterion (p ∼ 2.11). In the special case of d = 0, the Simon's criterion is calculated asν 2 − = 1 − (j−1)p−2j 2j(p−1) . Therefore, we need j > 1 for entanglement in the p → ∞ limit, instead of larger j(>2) for HZ1. When γ p /γ s = 10 3 , p = 100 and d = 100, we show the j dependency ofν − − 1 and the normalized HZ1 in Figure A1b. The required j for above-threshold entanglement was also by a factor of 2 smaller for Simon's criterion. Next, we show the results for large dissipative coupling. For γ p /γ s = 4, d = 1, j = 10 3 , and δ = 10 −10 , we plot theν − − 1 in Figure A1c. The Simon's entanglement criterion is satisfied at the same p(∼2.46) as the HZ1 criterion. From the recursion relation of ρ N s , the distribution of signal photon number ρ N s is obtained as For G = γ p , and G/γ s = 100, we show the comparison between numerical results of QME (calculated in the same methods as the main text) and the analytical results obtained by Equations (A52) and (A53). In Figure A2, two methods produced the same mean signal photon number, and the similar behavior of second order correlation function g (2) s (0) and Mandel's Q parameter Q M,s . Both methods show that the signal mode has a photon antibunching state below the threshold and an amplitude squeezed state above the threshold.