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

Abstract

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.

1. Introduction

Networks of degenerate optical parametric oscillators (DOPOs), called coherent Ising machines (CIMs) [1,2,3,4,5,6,7,8,9,10,11,12,13], have been extensively studied from quantum optics and neural-network perspectives (for a recent review, see Ref. [14]). A DOPO network is constructed with a dissipative (Liouvillian) linear coupling rather than a conservative (Hamiltonian) coupling using either optical delay lines [2,5,6] or homodyne measurement feedback circuits [7,8,12]. The fundamental topics in quantum optics that can be studied with CIMs include Gaussian state entanglement [3,4], the Schrödinger cat state [10], and the entangled cat state [15]. Applications cover a broad spectrum, including spin-glass solvers [16,17], structure-based virtual screening for drug discovery [18], combinatorial optimization [12,19], compressed sensing [20], and fair sampling for deep machine learning [21].

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. Appendix B and Appendix C provide subsidiary information on the numerical and theoretical equations. Appendix D and Appendix E provide the detailed derivations of analytical results. Appendix F and Appendix G show the supplementary information about the entanglement criterion and the experimental NOPO [30].

2. Model

The density matrix equation of a single ${\chi }^{\left(2\right)}$ NOPO is written as

$$\begin{array}{c}\hfill \frac{\partial \widehat{\rho }}{\partial t}={\displaystyle \sum _{a=p,s,i}}{\gamma }_a([{\widehat{a}}_a,\widehat{\rho }{\widehat{a}}_a^{†}]+\mathrm{h}.\mathrm{c}.)+\epsilon [{\widehat{a}}_p^{†}-{\widehat{a}}_p,\widehat{\rho }]+\kappa [{\widehat{a}}_s^{†}{\widehat{a}}_i^{†}{\widehat{a}}_p{e}^{-i\Delta t}-{\widehat{a}}_p^{†}{\widehat{a}}_i{\widehat{a}}_s{e}^{+i\Delta t},\widehat{\rho }].\end{array}$$

(1)

Here, ${\widehat{a}}_p$, ${\widehat{a}}_s$, and ${\widehat{a}}_i$ are pump, signal, and idler modes, respectively. These three photon modes have cavity frequencies ${\omega }_a(a=p,s,i)$ and the linear dissipation rates ${\gamma }_a(a=p,s,i)$. $\Delta =({\omega }_p-{\omega }_s)-{\omega }_i$ is the detuning of the parametric interaction. $\epsilon$ is the strength of the coherent excitation of the pump mode. $\kappa$ is the strength of the ${\chi }^{\left(2\right)}$ parametric interaction. If the idler mode has a much larger dissipation than the other two photon modes (${\gamma }_i\gg {\gamma }_p,{\gamma }_s$), the quantum master equation can be written as

$$\begin{array}{ccc}\hfill {\left.\frac{\partial \widehat{\rho }}{\partial t}\right|}_{NOPO}& =& {\displaystyle \sum _{a=p,s}}{\gamma }_a([{\widehat{a}}_a,\widehat{\rho }{\widehat{a}}_a^{†}]+\mathrm{h}.\mathrm{c}.)+\epsilon [{\widehat{a}}_p^{†}-{\widehat{a}}_p,\widehat{\rho }]\hfill \\ & -& iK[{\widehat{a}}_p^{†}{\widehat{a}}_p{\widehat{a}}_s^{†}{\widehat{a}}_s,\widehat{\rho }]+G([{\widehat{a}}_s^{†}{\widehat{a}}_p,\widehat{\rho }{\widehat{a}}_p^{†}{\widehat{a}}_s]+\mathrm{h}.\mathrm{c}.).\hfill \end{array}$$

(2)

Here G is the coefficient of the effective Raman interaction in Refs. [43,44,45], and K is the coefficient of the non-degenerate Kerr effect [46]. If $\Delta =0$, the Kerr coefficient is $K=0$ and the signal mode has a parametric gain $G={\kappa }^2/{\gamma }_i$. We define the maximum parametric gain as $G_0:={\kappa }^2/{\gamma }_i$. For the detuned parametric interaction, the parametric gain is $G=\frac{{\kappa }^2{\gamma }_i}{{\gamma }_i^2+{\Delta }^2}$, and the Kerr coefficient is $K=\frac{{\kappa }^2\Delta }{{\gamma }_i^2+{\Delta }^2}$. We use the dimensionless detuning parameter $d=K^2/G^2$ to represent the normalized detuning. We also introduce the normalized excitation $p=\epsilon /{\epsilon }_{thr}$, where ${\epsilon }_{thr}:={\gamma }_p\sqrt{{\gamma }_s/G}$ is the strength of the excitation at the oscillation threshold.

Let us consider a CXM consisting of two NOPOs (NOPO1 with ${\widehat{a}}_{p1},{\widehat{a}}_{s1}$ and NOPO2 with ${\widehat{a}}_{p2},{\widehat{a}}_{s2}$) coupled by a dissipative Liouvillian ${\mathcal{L}}_c\widehat{\rho }$ for two signal modes [3,4,47,48,49],

$$\frac{\partial \widehat{\rho }}{\partial t}={\left.\frac{\partial \widehat{\rho }}{\partial t}\right|}_{NOPO1}+{\left.\frac{\partial \widehat{\rho }}{\partial t}\right|}_{NOPO2}+{\mathcal{L}}_c\widehat{\rho },$$

(3)

$${\mathcal{L}}_c\widehat{\rho }=J[{\widehat{a}}_{s1}-{\widehat{a}}_{s2},\widehat{\rho }({\widehat{a}}_{s1}^{†}-{\widehat{a}}_{s2}^{†})]+\mathrm{h}.\mathrm{c}..$$

(4)

We will evaluate the entanglement between the two signal modes using one of the entanglement criteria in Ref. [50],

$$HZ1=|\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s2}\rangle {|}^2-\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s1}{\widehat{a}}_{s2}^{†}{\widehat{a}}_{s2}\rangle .$$

(5)

If this value is larger than zero, the two signal modes are entangled. Above the threshold, we assume the ferromagnetic configuration $\langle {\widehat{a}}_{s1}\rangle =\langle {\widehat{a}}_{s2}\rangle$ is created by the coupling Liouvillian (4). Using the fluctuation of positive-P amplitudes (introduced in Appendix A and Appendix B), the normalized $HZ1$ is written as $\frac{HZ1}{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }\sim -2\langle \Delta {\alpha }_{s1}^{†}\Delta {\alpha }_{s1}\rangle -2\mathrm{Re}\langle \Delta {\alpha }_{s1}\Delta {\alpha }_{s2}\rangle$. Here, we assumed $\langle \Delta {\alpha }_{s1}^{†}\Delta {\alpha }_{s1}\rangle =\langle \Delta {\alpha }_{s2}^{†}\Delta {\alpha }_{s2}\rangle$ and that $\langle {\widehat{a}}_{s1}\rangle$ is real. Introducing the canonical coordinates and canonical momenta as $\Delta X_{si}=\frac{\Delta {\alpha }_{si}+\Delta {\alpha }_{si}^{†}}{\sqrt{2}}$ and $\Delta P_{si}=\frac{\Delta {\alpha }_{si}-\Delta {\alpha }_{si}^{†}}{\sqrt{2}i}(i=1,2)$, this normalized entanglement criterion is written as

$$\frac{HZ1}{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }=-\langle \Delta X_{s1}^2\rangle -\langle \Delta X_{s1}\Delta X_{s2}\rangle -\langle \Delta P_{s1}^2\rangle +\langle \Delta P_{s1}\Delta P_{s2}\rangle .$$

(6)

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 $|\psi \rangle =\frac{1}{\sqrt{2}}\left(\right|0,1\rangle +|1,0\rangle )$ [52]. This criterion seems useful to detect non-Gaussian entanglement in the below-threshold CXM, since the similar model [33] realizes a single photon state.

3. Large Pump Dissipation Limit

3.1. Quantum Master Equation after the Elimination of the Pump Mode

In this section, we consider the large pump mode dissipation limit, where $G,K,{\gamma }_p\gg {\gamma }_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:

$$\widehat{\rho }=\sum _{N_s,N_s^{\prime }}\int P_{N_s,N_s^{\prime }}({\alpha }_p,{\alpha }_p^{†})\frac{|{\alpha }_p\rangle \langle {\alpha }_p^{†*}|}{\langle {\alpha }_p^{†*}|{\alpha }_p\rangle }\otimes |N_s\rangle \langle N_s^{\prime }|d{\alpha }_{p}d{\alpha }_p^{†}.$$

(7)

Substituting this expansion into Equation (2), we can obtain a time development equation for $P_{N_s,N_s^{\prime }}({\alpha }_p,{\alpha }_p^{†})$.

$$\begin{array}{ccc}\hfill \frac{\partial P_{N_s,N_s^{\prime }}}{\partial t}& =& {\gamma }_s[2\sqrt{(1+N_s)(1+N_s^{\prime })}{P}_{N_s+1,N_s^{\prime }+1}-(N_s+N_s^{\prime })P_{N_s,N_s^{\prime }}]-\epsilon \frac{\partial P_{N_s,N_s^{\prime }}}{\partial {\alpha }_p}-\epsilon \frac{\partial P_{N_s,N_s^{\prime }}}{\partial {\alpha }_p^{†}}\hfill \\ & +& \frac{\partial }{\partial {\alpha }_p}({\gamma }_p+G(1+N_s)){\alpha }_p{P}_{N_s,N_s^{\prime }}+iK\frac{\partial }{\partial {\alpha }_p}{N}_s{\alpha }_p{P}_{N_s,N_s^{\prime }}\hfill \\ & +& \frac{\partial }{\partial {\alpha }_p^{†}}({\gamma }_p+G(1+N_s^{\prime })){\alpha }_p^{†}{P}_{N_s,N_s^{\prime }}-iK\frac{\partial }{\partial {\alpha }_p^{†}}{N}_s^{\prime }{\alpha }_p^{†}{P}_{N_s,N_s^{\prime }}\hfill \\ & -& iK{\alpha }_p^{†}{\alpha }_p(N_s-N_s^{\prime })P_{N_s,N_s^{\prime }}+G{\alpha }_p^{†}{\alpha }_p[2\sqrt{N_s{N}_s^{\prime }}{P}_{N_s-1,N_s^{\prime }-1}-(2+N_s+N_s^{\prime })P_{N_s,N_s^{\prime }}].\hfill \end{array}$$

(8)

Due to the linear dissipation of the signal mode (${\gamma }_s$), components with photon number indices $(N_s,N_s^{\prime })$ are excited by components with larger photon number indices $(N_s+1,N_s^{\prime }+1)$. The parametric gain ($G{\alpha }_p^{†}{\alpha }_p$) introduces a contribution from components with smaller photon number indices $(N_s-1,N_s^{\prime }-1)$. Other processes do not affect the photon number indices. Although the equation has drift terms for the pump amplitudes $({\alpha }_p,{\alpha }_p^{†})$, the signal photon number indices do not change when components are derived using ${\alpha }_p$ or ${\alpha }_p^{†}$. When ${\gamma }_p+G$ is sufficiently large, the complex-P amplitude (${\alpha }_p$) at photon number indices ($N_s$, $N_s^{\prime }$) rapidly converge to the steady-state value of the time-development equation $\frac{d{\alpha }_p}{dt}=-({\gamma }_p+G(1+N_s)+iK{N}_s){\alpha }_p+\epsilon$. We can eliminate the pump mode amplitudes by writing $P_{N_s,N_s^{\prime }}({\alpha }_p,{\alpha }_p^{†})={\rho }_{N_s,N_s^{\prime }}\delta \left({\alpha }_p-\frac{\epsilon }{{\gamma }_p+G(1+N_s)+iK{N}_s}\right)\delta \left({\alpha }_p^{†}-\frac{\epsilon }{{\gamma }_p+G(1+N_s^{\prime })-iK{N}_s^{\prime }}\right)$, and integrating the complex-P amplitudes by $\int d{\alpha }_{p}d{\alpha }_p^{†}$. The density matrix components of the signal mode ${\rho }_{N_s,N_s^{\prime }}=\langle N_s|\widehat{\rho }|N_s^{\prime }\rangle$ are,

$$\begin{array}{ccc}\hfill \frac{\partial {\rho }_{N_s,N_s^{\prime }}}{\partial t}& =& 2{\gamma }_s\sqrt{(N_s+1)(N_s^{\prime }+1)}{\rho }_{N_s+1,N_s^{\prime }+1}-{\gamma }_s(N_s+N_s^{\prime }){\rho }_{N_s,N_s^{\prime }}\hfill \\ & -& i{K}_e(N_s,N_s^{\prime })(N_s-N_s^{\prime }){\rho }_{N_s,N_s^{\prime }}\hfill \\ & +& 2G_e(N_s-1,N_s^{\prime }-1)\sqrt{N_s{N}_s^{\prime }}{\rho }_{N_s-1,N_s^{\prime }-1}-G_e(N_s,N_s^{\prime })(2+N_s+N_s^{\prime }){\rho }_{N_s,N_s^{\prime }}.\hfill \end{array}$$

(9)

Here $G_e(N,N^{\prime })=\frac{G{\epsilon }^2}{[{\gamma }_p+G(1+N)+iKN][{\gamma }_p+G(1+N^{\prime })-iK{N}^{\prime }]}$, and $K_e(N,N^{\prime })=\frac{K{\epsilon }^2}{[{\gamma }_p+G(1+N)+iKN][{\gamma }_p+G(1+N^{\prime })-iK{N}^{\prime }]}$. 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 ${\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}=\langle N_1,N_2|\widehat{\rho }|N_1^{\prime },N_2^{\prime }\rangle$ develop as,

$$\begin{array}{ccc}\hfill & & \frac{\partial {\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}}{\partial t}=({\gamma }_s+J)\left[2\sqrt{(N_1+1)(N_1^{\prime }+1)}{\rho }_{N_1+1,N_2,N_1^{\prime }+1,N_2^{\prime }}-(N_1+N_1^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}\right]\hfill \\ & +& 2G_e(N_1-1,N_1^{\prime }-1)\sqrt{N_1{N}_1^{\prime }}{\rho }_{N_1-1,N_2,N_1^{\prime }-1,N_2^{\prime }}-G_e(N_1,N_1^{\prime })(2+N_1+N_1^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}\hfill \\ & +& ({\gamma }_s+J)\left[2\sqrt{(N_2+1)(N_2^{\prime }+1)}{\rho }_{N_1,N_2+1,N_1^{\prime },N_2^{\prime }+1}-(N_2+N_2^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}\right]\hfill \\ & +& 2G_e(N_2-1,N_2^{\prime }-1)\sqrt{N_2{N}_2^{\prime }}{\rho }_{N_1,N_2-1,N_1^{\prime },N_2^{\prime }-1}-G_e(N_2,N_2^{\prime })(2+N_2+N_2^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}\hfill \\ & -& i{K}_e(N_1,N_1^{\prime })(N_1-N_1^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}-i{K}_e(N_2,N_2^{\prime })(N_2-N_2^{\prime }){\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}\hfill \\ & -& J[2\sqrt{(N_1+1)(N_2^{\prime }+1)}{\rho }_{N_1+1,N_2,N_1^{\prime },N_2^{\prime }+1}+2\sqrt{(N_2+1)(N_1^{\prime }+1)}{\rho }_{N_1,N_2+1,N_1^{\prime }+1,N_2^{\prime }}\hfill \\ & -& \sqrt{(N_1+1)N_2}{\rho }_{N_1+1,N_2-1,N_1^{\prime },N_2^{\prime }}-\sqrt{N_1^{\prime }(N_2^{\prime }+1)}{\rho }_{N_1,N_2,N_1^{\prime }-1,N_2^{\prime }+1}\hfill \\ & -& \sqrt{(N_2+1)N_1}{\rho }_{N_1-1,N_2+1,N_1^{\prime },N_2^{\prime }}-\sqrt{N_2^{\prime }(N_1^{\prime }+1)}{\rho }_{N_1,N_2,N_1^{\prime }+1,N_2^{\prime }-1}].\hfill \end{array}$$

(10)

3.2. 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\ll 1$). In this limit, $G_e(N_i,N_i^{\prime })$ and $K_e(N_i,N_i^{\prime })(i=1,2)$ are of order $O\left(p^2\right)$. Therefore, these contributions to ${\rho }_{N_1,N_2,N_1^{\prime },N_2^{\prime }}$ on the right hand side of Equation (10) are much smaller than those of ${\gamma }_s$ and J. However, the $G_e(N_i-1,N_i^{\prime }-1)(i=1,2)$ on the right hand side are not negligible. Although the $K_e(N_i,N_i^{\prime })(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^{\prime }-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, ${\rho }_{10,10}$, and ${\rho }_{10,01}$.

$$\frac{\partial {\rho }_{10,10}}{\partial t}=-2({\gamma }_s+J){\rho }_{10,10}+2J{\rho }_{10,01}+2G_e(0,0){\rho }_{00,00},$$

(11)

$$\frac{\partial {\rho }_{10,01}}{\partial t}=-2({\gamma }_s+J){\rho }_{10,01}+2J{\rho }_{10,10}.$$

(12)

From these equations, the mean signal photon number is ${\rho }_{10,10}/{\rho }_{00,00}$. In the steady-state, this becomes

$$\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s1}\rangle =\frac{{\gamma }_s+J}{{\gamma }_s+2J}\frac{{\gamma }_p^2}{{({\gamma }_p+G)}^2}{p}^2.$$

(13)

For a weak nonlinear gain saturation ($G\ll {\gamma }_p$), this value reaches one at $p\sim 1$. However, in general, $p\sim 1+G/{\gamma }_p$ is required for achieving $\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s1}\rangle =1$, where stimulated emission becomes dominant over spontaneous emission. ${\rho }_{10,01}/{\rho }_{00,00}$ provides an amplitude correlation function between the two signal modes,

$$\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s2}\rangle =\frac{J}{{\gamma }_s+2J}\frac{{\gamma }_p^2}{{({\gamma }_p+G)}^2}{p}^2.$$

(14)

Next, we calculate the values of order $O\left(p^4\right)$ from the equations of four components, ${\rho }_{20,20}$, ${\rho }_{11,11}$, ${\rho }_{20,02}$ and ${\rho }_{20,11}$.

$$\frac{\partial {\rho }_{20,20}}{\partial t}=-4({\gamma }_s+J){\rho }_{20,20}+2\sqrt{2}J\mathrm{Re}{\rho }_{20,11}+4G_e(1,1){\rho }_{10,10},$$

(15)

$$\frac{\partial {\rho }_{11,11}}{\partial t}=-4({\gamma }_s+J){\rho }_{11,11}+4\sqrt{2}J\mathrm{Re}{\rho }_{20,11}+4G_e(0,0){\rho }_{10,10},$$

(16)

$$\frac{\partial {\rho }_{20,02}}{\partial t}=-4({\gamma }_s+J){\rho }_{20,02}+2\sqrt{2}J\mathrm{Re}{\rho }_{20,11},$$

(17)

$$\frac{\partial {\rho }_{20,11}}{\partial t}=-4({\gamma }_s+J){\rho }_{20,11}+\sqrt{2}J({\rho }_{20,02}+{\rho }_{11,11}+{\rho }_{20,20})+2\sqrt{2}{G}_e(1,0){\rho }_{10,01}.$$

(18)

The second-order correlation function $g_s^{\left(2\right)}\left(0\right)=\frac{\langle {\widehat{a}}_{s1}^{†2}{\widehat{a}}_{s1}^2\rangle }{{\langle {\widehat{a}}_{s1}^{†}{\widehat{a}}_{s1}\rangle }^2}$ is obtained as $g_s^{\left(2\right)}\left(0\right)=2{\rho }_{20,20}{\rho }_{00,00}/{\rho }_{10,10}^2$. Its steady-state value is

$$g_s^{\left(2\right)}\left(0\right)=\frac{4{\gamma }_s({\gamma }_s+2J){({\gamma }_p+G)}^2+J^2[K^2+{(2{\gamma }_p+3G)}^2]}{2{({\gamma }_s+J)}^2[{({\gamma }_p+2G)}^2+K^2]}.$$

(19)

When the parametric coefficients are sufficiently small ($G,K\ll {\gamma }_p$), the signal modes are in blackbody radiation states with $g_s^{\left(2\right)}\left(0\right)=2$ [29]. In the case of a large parametric gain ($G,K$), however, the signal modes can have non-classical (photon anti-bunching) states with $g_s^{\left(2\right)}\left(0\right)<1$. In the single NOPO limit ($J\to 0$), $g_s^{\left(2\right)}\left(0\right)=\frac{2{({\gamma }_p+G)}^2}{{({\gamma }_p+2G)}^2+K^2}$ is obtained. This is identical to the ${\gamma }_s\to 0$ limit of $g_s^{\left(2\right)}\left(0\right)$ for a single NOPO, which is derived in Appendix A and shown in Equation (A15). This expression converges to a completely anti-bunching state $g_s^{\left(2\right)}\left(0\right)=0$ in the $K\to \infty$ limit [33]. In the ${\gamma }_s\ll J$ limit, $g_s^{\left(2\right)}\left(0\right)$ is $\frac{K^2+{(2{\gamma }_p+3G)}^2}{2[K^2+{({\gamma }_p+2G)}^2]}$, and it converges to $0.5$ in the large-K limit. This is larger than the minimum value for a single NOPO ($g_s^{\left(2\right)}\left(0\right)=0$). Hillery-Zubairy’s entanglement criterion (Equation (5)) for the two signal modes is obtained as $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2={\rho }_{10,01}^2/{\rho }_{10,10}^2-{\rho }_{11,11}{\rho }_{00,00}/{\rho }_{10,10}^2$, i.e.,

$$\frac{HZ1}{{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2}=\frac{J^2[K^2-(2{\gamma }_p^2+4{\gamma }_{p}G+G^2)]}{2{({\gamma }_s+J)}^2[{({\gamma }_p+2G)}^2+K^2]}-\frac{{\gamma }_s({\gamma }_s+2J)}{{({\gamma }_s+J)}^2}.$$

(20)

If this value is larger than zero, the two NOPOs are entangled far below the threshold. From the equation, $HZ1$ is always negative if $K\le G$. To achieve entanglement, the detuning of the parametric interaction must be larger than the idler linewidth ($d>1$). For below-threshold entanglement, a small ${\gamma }_p/G$ ratio and large J are preferred.

The following is an example of a parameter set for below-threshold entanglement, $J/{\gamma }_s=12$, $G=8{\gamma }_p$, and $d=5$. These parameters give $g_s^{\left(2\right)}\left(0\right)\sim 0.7361$ and $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2\sim 0.0074$. Therefore, the two NOPOs have anti-bunching states and are entangled. Figure 1a,b compares the analytically and numerically calculated $g_s^{\left(2\right)}\left(0\right)$ and $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2$. The numerical results are obtained by time developing the quantum master Equation (3) from the vacuum state $\widehat{\rho }=|0\rangle \langle 0|$ in accordance with the pump schedule $p\left(t\right)=0.01min(1,\sqrt{t{\gamma }_s/2})$. The values at $t{\gamma }_s=10$ are used as the approximate steady-state results for $p=0.01$, and are plotted for various ${\gamma }_p/{\gamma }_s$ ratio. When ${\gamma }_p/{\gamma }_s$ becomes larger than $J/{\gamma }_s=12$ by two orders of magnitude, $g_s^{\left(2\right)}\left(0\right)$ and the below-threshold entanglement criterion converge to the analytical results that assume infinite ${\gamma }_p/{\gamma }_s$.

3.3. Far-Above-Threshold Entanglement

Next, we analytically derive the above-threshold characteristics of the CXM by assuming a large pump dissipation ${\gamma }_p\gg {\gamma }_s$ and small parametric gain $G\ll {\gamma }_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 quasi-distribution function $P({\alpha }_p,{\alpha }_p^{†},{\alpha }_s,{\alpha }_s^{†})$ is derived in Appendix A. Here, ${\alpha }_p,{\alpha }_p^{†}({\alpha }_s,{\alpha }_s^{†})$ are positive-P amplitudes for the pump (signal) mode. Using the Ito rule, the equivalent c-number stochastic differential equations are expressed as:

$$\frac{d{\alpha }_p}{dt}=-({\gamma }_p+G){\alpha }_p+\epsilon -(G+iK){\alpha }_s^{†}{\alpha }_s{\alpha }_p-\sqrt{\frac{G+iK}{2}}{\alpha }_p{\xi }_C,$$

(21)

$$\frac{d{\alpha }_p^{†}}{dt}=-({\gamma }_p+G){\alpha }_p^{†}+\epsilon -(G-iK){\alpha }_s^{†}{\alpha }_s{\alpha }_p^{†}-\sqrt{\frac{G-iK}{2}}{\alpha }_p^{†}{\xi }_C^{†},$$

(22)

$$\frac{d{\alpha }_s}{dt}=-{\gamma }_s{\alpha }_s+(G-iK){\alpha }_p^{†}{\alpha }_p{\alpha }_s+\sqrt{\frac{G+iK}{2}}{\alpha }_s{\xi }_C^{*}+\sqrt{G}{\alpha }_p{\xi }_C^{‡},$$

(23)

$$\frac{d{\alpha }_s^{†}}{dt}=-{\gamma }_s{\alpha }_s^{†}+(G+iK){\alpha }_p^{†}{\alpha }_p{\alpha }_s^{†}+\sqrt{\frac{G-iK}{2}}{\alpha }_s^{†}{\xi }_C^{†*}+\sqrt{G}{\alpha }_p^{†}{\xi }_C^{‡*}.$$

(24)

Here ${\xi }_C$, ${\xi }_C^{†}$ and ${\xi }_C^{‡}$ are independent complex Gaussian random numbers satisfying correlation functions, such as $\langle {\xi }_C^{*}\left(t\right){\xi }_C\left(t^{\prime }\right)\rangle =2\delta (t-t^{\prime })$. ${\xi }_C^{‡}$ reflects random spontaneous emission of signal photons. ${\xi }_C$ and ${\xi }_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]:

$$\langle {\alpha }_p^{†}{\alpha }_p\rangle ={\gamma }_s/G.$$

(25)

The above-threshold signal photon number is obtained from ${\gamma }_s/G\sim {\left|\langle {\alpha }_p\rangle \right|}^2$ and

$$\langle {\alpha }_p\rangle \sim \frac{\epsilon }{{\gamma }_p+(G+iK)\langle {\alpha }_s^{†}{\alpha }_s\rangle },$$

(26)

which is the steady-state mean of Equation (21), $\frac{d\langle {\alpha }_p\rangle }{dt}\sim \epsilon -{\gamma }_p\langle {\alpha }_p\rangle -(G+iK)\langle {\alpha }_p\rangle \langle {\alpha }_s^{†}{\alpha }_s\rangle$. Here, we have assumed that $G\ll {\gamma }_p$ and $\langle {\alpha }_p{\alpha }_s^{†}{\alpha }_s\rangle \sim \langle {\alpha }_p\rangle \langle {\alpha }_s^{†}{\alpha }_s\rangle$. Using $d=K^2/G^2$ and $\epsilon ={\gamma }_{p}p\sqrt{{\gamma }_s/G}$, the steady-state signal photon number is obtained as

$$\langle {\alpha }_s^{†}{\alpha }_s\rangle =\frac{{\gamma }_p}{G}\frac{\pi -1}{1+d}.$$

(27)

Here, $\pi =\sqrt{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

$$\langle {\alpha }_p\rangle =\sqrt{\frac{{\gamma }_s}{G}}p\frac{1-i\sqrt{d}}{\pi -i\sqrt{d}}.$$

(28)

By introducing the phase factor

$$tan\varphi =\frac{\sqrt{d}(\pi -1)}{\pi +d},$$

(29)

the mean pump amplitude becomes $\langle {\alpha }_p\rangle =\sqrt{\frac{{\gamma }_s}{G}}{e}^{-i\varphi }$.

Now let us introduce the small amplitude fluctuations in the pump and signal modes, $\Delta {\alpha }_p$ and $\Delta {\alpha }_s$. The pump amplitude fluctuation is defined as the fluctuation after removing the phase factor denoted by $\varphi$.

$${\alpha }_p\sim e^{-i\varphi }\left(\sqrt{\frac{{\gamma }_s}{G}}+\Delta {\alpha }_p\right).$$

(30)

The mean signal amplitude rotates with the frequency $K\langle {\alpha }_p^{†}{\alpha }_p\rangle \sim \sqrt{d}{\gamma }_s$ due to Kerr-nonlinearity. The signal amplitude fluctuation is defined as the fluctuation after removing this time-dependent phase factor,

$${\alpha }_s\sim e^{-i\sqrt{d}{\gamma }_{s}t}\left(\sqrt{\frac{{\gamma }_p}{G}\frac{\pi -1}{1+d}}+\Delta {\alpha }_s\right).$$

(31)

These phase factors do not affect the products, ${\alpha }_p^{†}{\alpha }_p$ and ${\alpha }_s^{†}{\alpha }_s$, appearing in the positive-P equations. The equations for small amplitude fluctuations are

$$\frac{d\Delta {\alpha }_p}{dt}=-r{p}^2\frac{1+i\sqrt{d}}{\pi +i\sqrt{d}}\Delta {\alpha }_p-(1+i\sqrt{d})\sqrt{\frac{r(\pi -1)}{1+d}}(\Delta {\alpha }_s^{†}+\Delta {\alpha }_s)-\sqrt{\frac{1+i\sqrt{d}}{2}}{\xi }_C,$$

(32)

$$\frac{d\Delta {\alpha }_p^{†}}{dt}=-r{p}^2\frac{1-i\sqrt{d}}{\pi -i\sqrt{d}}\Delta {\alpha }_p^{†}-(1-i\sqrt{d})\sqrt{\frac{r(\pi -1)}{1+d}}(\Delta {\alpha }_s^{†}+\Delta {\alpha }_s)-\sqrt{\frac{1-i\sqrt{d}}{2}}{\xi }_C^{†},$$

(33)

$$\frac{d\Delta {\alpha }_s}{dt}=+(1-i\sqrt{d})\sqrt{\frac{r(\pi -1)}{1+d}}(\Delta {\alpha }_p+\Delta {\alpha }_p^{†})+\sqrt{\frac{1+i\sqrt{d}}{2}\frac{r(\pi -1)}{1+d}}{\xi }_C^{*}+{\xi }_C^{{‡}^{\prime }},$$

(34)

$$\frac{d\Delta {\alpha }_s^{†}}{dt}=+(1+i\sqrt{d})\sqrt{\frac{r(\pi -1)}{1+d}}(\Delta {\alpha }_p+\Delta {\alpha }_p^{†})+\sqrt{\frac{1-i\sqrt{d}}{2}\frac{r(\pi -1)}{1+d}}{\xi }_C^{†*}+{\xi }_C^{{‡}^{\prime }*}.$$

(35)

Here, $r={\gamma }_p/{\gamma }_s$ and ${\xi }_C^{{‡}^{\prime }}={\xi }_C^{‡}{e}^{-i\varphi }{e}^{i\sqrt{d}{\gamma }_{s}t}$. The time t is normalized so that $1/{\gamma }_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\gg 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}:=\frac{\langle \Delta {\widehat{n}}_s^2\rangle -\langle {\widehat{n}}_s\rangle }{\langle {\widehat{n}}_s\rangle }$, where ${\widehat{n}}_s={\widehat{a}}_s^{†}{\widehat{a}}_s$ and $\Delta {\widehat{n}}_s={\widehat{n}}_s-\langle {\widehat{n}}_s\rangle$. 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,

$$Q_{M,s}=-\frac{4+j}{4(2+j)}+\frac{1+d}{4(\pi -1)}-\frac{d}{4\pi }+\frac{(1+j)(\pi +d)}{(2+j)[(2+j){\pi }^2-2\pi +dj]}$$

(36)

is obtained. In the large excitation limit ($\pi \to \infty$), $Q_{M,s}$ converges to a negative value $-\frac{4+j}{4(2+j)}$. Therefore, an amplitude squeezed state is obtained even in a small-G NOPO. In the single NOPO limit ($j\to 0$), $Q_{M,s}=-\frac{1}{2}+\frac{1+d}{2(\pi -1)}-\frac{d}{2\pi }$. This value converges to $Q_{M,s}\to -\frac{1}{2}$ for a large pump excitation, which is identical to the value obtained in Ref. [35]. In terms of the order $O\left(d\right)$, Mandel’s Q parameter is written as $Q_{M,s}=-\frac{1}{2}+\frac{1}{2(p-1)}-\frac{d}{4p}$. Therefore, for the same p, detuning in the parametric interaction increases the amplitude squeezing. In the strong dissipative coupling limit, i.e., $j\to \infty$, 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

$$\frac{HZ1}{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }=\frac{j-2}{4j}-\frac{\pi +d}{4\pi (\pi -1)}-\frac{d(\pi -1)}{j({\pi }^2+d)}\left[1+\frac{(d+1)\pi }{(j+1){\pi }^2-\pi +dj}+\frac{\pi ({\pi }^2-2\pi -d)}{(j+2){\pi }^2-2\pi +dj}\right].$$

(37)

The detailed derivation is provided in Appendix D. Here, in the large excitation limit ($\pi \to \infty$), $HZ1$ converges as $HZ1/\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \to \frac{j^2-4(1+d)}{4j(j+2)}$. The coupling coefficient j must satisfy $j>2\sqrt{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/{\gamma }_s=12$ and $d=5$ and plotted as the function of normalized excitation p. The analytical results assume that ${\gamma }_p/{\gamma }_s\to \infty$ and $G/{\gamma }_s\to 0$. The numerical results were obtained by integrating the positive-Pc-number SDEs in Appendix B, for ${\gamma }_p/{\gamma }_s=50$ and $G/{\gamma }_s={10}^{-7}$. To obtain steady-state statistics from the positive-P calculation, we calculated a single trajectory, starting from ${\alpha }_{pi}={\alpha }_{pi}^{†}={\alpha }_{si}={\alpha }_{si}^{†}=0(i=1,2)$ with excitation depending on time as $p\left(t\right)=pmin(1,\sqrt{t{\gamma }_s/{10}^5})$, and took the time average from $t{\gamma }_s={10}^5$ to $t{\gamma }_s={10}^6$. The theoretical Mandel’s $Q_{M,s}$ parameter from Equation (36) becomes negative at $p\sim 1.7$. The numerical $Q_{M,s}$ values were slightly larger than the theoretical values due to the finite ${\gamma }_p/{\gamma }_s$ ratio. When $d=0$, Equation (36) always gave $Q_{M,s}=0$ at $p=2$. Therefore, the non-degenerate Kerr effect causes amplitude squeezing for smaller p. Moreover, the theoretical $HZ1$ from Equation (37) becomes positive for $p\sim 2.1$. The numerical values are slightly smaller than the analytical results due to the finite ${\gamma }_p/{\gamma }_s$ ratio. Above the threshold, the use of a small ${\gamma }_p/{\gamma }_s$ ratio helps to decrease the degrees of both amplitude squeezing and entanglement. This differs from the below threshold case, where using small ${\gamma }_p/{\gamma }_s$ ratio improves the anti-bunching and entanglement.

3.4. Numerical Results

Here, we give the numerical simulation results as a function of normalized excitation p. The mean signal photon number $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle$ and $HZ1$ entanglement criterion are presented in Figure 2a,b for ${\gamma }_p/{\gamma }_s=50$, $G/{\gamma }_s=400$, $d=5$, and $J/{\gamma }_s=12$. For a small excitation p, we numerically calculated the quantum master equation (QME) in Equation (3) with a photon number space expansion. When the maximum pump photon number is $M_p-1$ and maximum signal photon number is $M_s-1$ in the expansion, the size of the density matrix and the calculation time are of order $O\left(M_p^4{M}_s^4\right)$. This means that the calculation time increases rapidly for large p, where the mean signal photon number $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle$ increases and we have to prepare a sufficiently larger $M_s$ satisfying $M_s\gg \langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle$. The maximum value of p calculated by QME was $p=10$, where the mean signal photon number is $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim 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\left(M_p^2{M}_s^2\right)$. In the WFMC calculation, time development started from the vacuum state with the excitation $p\left(t\right)=pmin(1,\sqrt{t{\gamma }_s/{10}^2})$. A time average was taken from $t{\gamma }_s={10}^2$ to $t{\gamma }_s={10}^3-2.5\times {10}^4$, depending on the values of normalized excitation p. The calculation for small p with $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle <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\sim 178$, where the mean signal photon number was $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim 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\left(M_s^4\right)$, so we could use it to calculate even the above-threshold characteristics. We calculated the characteristics at $p\sim 316$ with $M_s=46$, where $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim 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\left(p^2\right)$ 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\sqrt{1+d}$ is satisfied). However, the small stimulated emission below the threshold gives a negative correction to $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2$. $HZ1$ turned negative and reached a minimum value near the threshold ($p\sim 19$), where the signal photon number $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle$ exceeds one. Far above the threshold $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^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/{\gamma }_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.

4. 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/{\gamma }_s$ to the same order of $G/{\gamma }_s$. In Figure 1a,b, it can be seen that a large linewidth ratio (${\gamma }_p/{\gamma }_s$) reduces the anti-bunching and entanglement. The theory for a large pump mode dissipation assumes the coupling coefficient J is much smaller than ${\gamma }_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 (${\gamma }_p,{\gamma }_s,G,K$).

4.1. 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 $\frac{\partial P}{\partial t}{|}_{NOPO,\epsilon \to 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,

$$\begin{array}{ccc}\hfill \frac{\partial P}{\partial t}& =& {\left.\frac{\partial P}{\partial t}\right|}_{NOPO1,\epsilon \to 0}+J\frac{\partial }{\partial {\alpha }_{s1}}({\alpha }_{s1}-{\alpha }_{s2})P+J\frac{\partial }{\partial {\alpha }_{s1}^{†}}({\alpha }_{s1}^{†}-{\alpha }_{s2}^{†})P\hfill \\ & +& {\left.\frac{\partial P}{\partial t}\right|}_{NOPO2,\epsilon \to 0}+J\frac{\partial }{\partial {\alpha }_{s2}}({\alpha }_{s2}-{\alpha }_{s1})P+J\frac{\partial }{\partial {\alpha }_{s2}^{†}}({\alpha }_{s2}^{†}-{\alpha }_{s1}^{†})P.\hfill \end{array}$$

(38)

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 $\langle {\alpha }_p^{†i}{\alpha }_p^j\rangle ={\left(\frac{{\gamma }_{p}p}{{\gamma }_p+G}\sqrt{\frac{{\gamma }_s}{G}}\right)}^{i+j}$, which is identical to those of a single NOPO. However, as $J\to \infty$, the signal photon number becomes one half of that of a single NOPO $\langle {\alpha }_{s1}^{†}{\alpha }_{s1}\rangle \sim \frac{{\gamma }_p^2{p}^2}{2{({\gamma }_p+G)}^2}$. This is the same as the $J\to \infty$ limit of Equation (13) obtained for large pump mode dissipation, and is also identical to the same limit of the two signal amplitudes’ correlation $\langle {\alpha }_{s1}^{†}{\alpha }_{s2}\rangle$ in Equation (14). We assume that such amplitude products connected by J, i.e., $\langle {\alpha }_{p1}^{†i}{\alpha }_{p1}^j{\alpha }_{s1}^{†k}{\alpha }_{s1}^l\rangle$, to any replacement of ${\alpha }_{s1}^{†}\to {\alpha }_{s2}^{†}$ and ${\alpha }_{s1}\to {\alpha }_{s2}$, have the same value in the zeroth order of $1/J$.

The details of the derivation of $g_s^{\left(2\right)}\left(0\right)$ are shown in Appendix E. In the large-J limit, the below-threshold second-order correlation function $g_s^{\left(2\right)}\left(0\right)$ can be written as,

$$g_s^{\left(2\right)}\left(0\right)=4\frac{{\gamma }_s+2({\gamma }_p+G)\mathrm{Re}\frac{{\gamma }_p+2{\gamma }_s+G}{2{\gamma }_p+4{\gamma }_s+3G+iK}}{2{\gamma }_p+2{\gamma }_s+3G}.$$

(39)

In the limit ${\gamma }_s\to 0$, this value converges to $g_s^{\left(2\right)}\left(0\right)=\frac{8{({\gamma }_p+G)}^2}{{(2{\gamma }_p+3G)}^2+K^2}$. This is different from the $J\to \infty$ limit of Equation (19). If we start from large dissipative coupling $J\to \infty$, the anti-bunching state $g_s^{\left(2\right)}\left(0\right)\to 0$ is obtained with large non-degenerate Kerr coefficient K. From Equation (39), the below-threshold entanglement criterion obeys as $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2=1-g_s^{\left(2\right)}\left(0\right)$. The normalized $HZ1$ with no detuning in the parametric interaction ($K=0$) is

$$\frac{HZ1}{{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2}=\frac{G^2-2G(2{\gamma }_p+5{\gamma }_s)-4({\gamma }_p+{\gamma }_s)({\gamma }_p+2{\gamma }_s)}{(2{\gamma }_p+2{\gamma }_s+3G)(2{\gamma }_p+4{\gamma }_s+3G)}.$$

(40)

This equation shows that the far-below-threshold entanglement criterion is satisfied even with $K=0$. When ${\gamma }_p={\gamma }_s$, it is satisfied for $G/{\gamma }_s>15.6$.

Figure 3a compares the normalized entanglement criterion $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2$ obtained from the quantum master Equation (3) and the theory in the $J\to \infty$ limit ($1-g_s^{\left(2\right)}\left(0\right)$ with Equation (39)). We set the linewidth ratio ${\gamma }_p/{\gamma }_s$ to 4 and the maximum parametric gain $G_0/{\gamma }_s={\kappa }^2/\left({\gamma }_i{\gamma }_s\right)$ to 40 and investigate the impact of normalized detuning d in the parametric interaction. With detuning, the parametric gain becomes $G=\frac{G_0}{1+d}$ and the non-degenerate Kerr effect is $K=\frac{G_0\sqrt{d}}{1+d}$. 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\to \infty$. Even in the $J\to \infty$ limit, normalized $HZ1$ slightly increases for a small detuning d. As d increases, K decreases as $\propto \frac{1}{\sqrt{d}}$ and gets closer to ${\gamma }_s$, where below-threshold entanglement is impossible. When $J/{\gamma }_s=360$, entanglement criterion is satisfied for $d=0$. However, when $J/{\gamma }_s=120$, the entanglement criterion is not satisfied with zero detuning, although $HZ1$ becomes positive with a non-zero detuning parameter.

4.2. Far-Above-Threshold Entanglement

Next, we derive the far-above-threshold entanglement in the limit $J\to \infty$. First of all, we will consider a single NOPO. From the cSDEs in Equations (32) and (33), we obtain the equations of the pump amplitude fluctuations, $\Delta X_p=\frac{\Delta {\alpha }_p+\Delta {\alpha }_p^{†}}{\sqrt{2}}$ and $\Delta P_p=\frac{\Delta {\alpha }_p-\Delta {\alpha }_p^{†}}{\sqrt{2}i}$. The equations in this representation are shown in Appendix C. The mean fluctuation products are

$$\begin{array}{c}\hfill \frac{d\langle \Delta X_p^2\rangle }{dt}=-2rpcos\varphi \langle \Delta X_p^2\rangle +2rpsin\varphi \langle \Delta X_p\Delta P_p\rangle -4\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta X_p\Delta X_s\rangle ,\end{array}$$

(41)

$$\begin{array}{c}\hfill \frac{d\langle \Delta P_p^2\rangle }{dt}=-2rpcos\varphi \langle \Delta P_p^2\rangle -2rpsin\varphi \langle \Delta X_p\Delta P_p\rangle -4\sqrt{\frac{rd(\pi -1)}{1+d}}\langle \Delta P_p\Delta X_s\rangle ,\end{array}$$

(42)

$$\begin{array}{ccc}\hfill \frac{d\langle \Delta X_p\Delta P_p\rangle }{dt}& =& -2rpcos\varphi \langle \Delta X_p\Delta P_p\rangle -rpsin\varphi \langle \Delta X_p^2\rangle +rpsin\varphi \langle \Delta P_p^2\rangle \hfill \\ & -& 2\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta P_p\Delta X_s\rangle -2\sqrt{\frac{rd(\pi -1)}{1+d}}\langle \Delta X_p\Delta X_s\rangle .\hfill \end{array}$$

(43)

These pump amplitude fluctuations are excited by two fluctuation correlations, $\langle \Delta X_p\Delta X_s\rangle$, and $\langle \Delta P_p\Delta X_s\rangle$, where $\Delta X_s$ is the amplitude fluctuation of the signal mode. We introduce the normalized amplitude correlations,

$$B=-2\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta X_p\Delta X_s\rangle ,$$

(44)

$$A=-2\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta P_p\Delta X_s\rangle .$$

(45)

The steady-state pump amplitude fluctuations in terms of A and B are

$$\langle \Delta X_p^2\rangle =\frac{1+d}{2r(\pi +d)}(B+\sqrt{d}A)+\frac{\pi B-\sqrt{d}A}{2r{p}^2},$$

(46)

$$\langle \Delta P_p^2\rangle =\frac{1+d}{2r(\pi +d)}(B+\sqrt{d}A)-\frac{\pi B-\sqrt{d}A}{2r{p}^2},$$

(47)

$$2\langle \Delta X_p\Delta P_p\rangle =\frac{\pi A+\sqrt{d}B}{r{p}^2}.$$

(48)

Next, we obtain the time development equations of $\langle \Delta X_p\Delta X_s\rangle$, $\langle \Delta P_p\Delta X_s\rangle$ and $\langle \Delta X_s^2\rangle$.

$$\begin{array}{ccc}\hfill \frac{d\langle \Delta X_p\Delta X_s\rangle }{dt}& =& -rpcos\varphi \langle \Delta X_p\Delta X_s\rangle +rpsin\varphi \langle \Delta P_p\Delta X_s\rangle \hfill \\ & +& 2\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta X_p^2\rangle -\sqrt{\frac{r(\pi -1)}{1+d}}(2\langle \Delta X_s^2\rangle +1),\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill \frac{d\langle \Delta P_p\Delta X_s\rangle }{dt}& =& -rpcos\varphi \langle \Delta P_p\Delta X_s\rangle -rpsin\varphi \langle \Delta X_p\Delta X_s\rangle \hfill \\ & +& 2\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta X_p\Delta P_p\rangle -\sqrt{\frac{rd(\pi -1)}{1+d}}(2\langle \Delta X_s^2\rangle +1),\hfill \end{array}$$

(50)

$$\frac{d\langle \Delta X_s^2\rangle }{dt}=4\sqrt{\frac{r(\pi -1)}{1+d}}\langle \Delta X_p\Delta X_s\rangle +2.$$

(51)

$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,

$$\langle \Delta X_s^2\rangle =-\frac{1}{4}+\frac{1+d}{4(\pi -1)}-\frac{A\sqrt{d}}{4}+\langle \Delta X_p^2\rangle .$$

(52)

Here,

$$A=\frac{\sqrt{d}}{\pi }\frac{\frac{r}{2(\pi -1)}+\frac{1}{\pi +d}+\frac{\pi -1}{{\pi }^2+d}}{\frac{r}{2(\pi -1)}+\frac{1}{\pi +d}-\frac{\pi -1}{{\pi }^2+d}}.$$

(53)

This fluctuation theory for a single NOPO is sufficient for calculating the normalized $HZ1$ in a CXM in the $j\gg 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

$$\begin{array}{c}\hfill \frac{d\langle \Delta P_{s1}^2\rangle }{dt}=-4\sqrt{\frac{rd(\pi -1)}{1+d}}\langle \Delta X_{p1}\Delta P_{s1}\rangle -2(j+\delta )\langle \Delta P_{s1}^2\rangle +2j\langle \Delta P_{s1}\Delta P_{s2}\rangle +2,\end{array}$$

(54)

$$\begin{array}{c}\hfill \frac{d\langle \Delta P_{s1}\Delta P_{s2}\rangle }{dt}=-4\sqrt{\frac{rd(\pi -1)}{1+d}}\langle \Delta X_{p1}\Delta P_{s2}\rangle -2(j+\delta )\langle \Delta P_{s1}\Delta P_{s2}\rangle +2j\langle \Delta P_{s1}^2\rangle .\end{array}$$

(55)

The steady-state difference of these two fluctuation values is

$$\begin{array}{c}\hfill \langle \Delta P_{s1}^2\rangle -\langle \Delta P_{s1}\Delta P_{s2}\rangle =\frac{1}{2j}-\frac{1}{j}\sqrt{\frac{rd(\pi -1)}{1+d}}(\langle \Delta X_{p1}\Delta P_{s1}\rangle -\langle \Delta X_{p1}\Delta P_{s2}\rangle ),\end{array}$$

(56)

which converges to zero in the large j limit. The sum $\langle \Delta X_{s1}^2\rangle +\langle \Delta X_{s1}\Delta X_{s2}\rangle$ appearing in the $HZ1$ is the same as $\langle \Delta X_s^2\rangle$ of a single NOPO shown in Equation (52). Finally, the normalized $HZ1$ in the $J\to \infty$ limit is

$$\begin{array}{c}\hfill \frac{HZ1}{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }=\frac{1}{4}-\frac{1+d}{4(\pi -1)}+\frac{\sqrt{d}A}{4}-\frac{(1+d)(1+\sqrt{d}A)}{2r(\pi +d)}-\frac{\pi -\sqrt{d}A}{2r{p}^2}.\end{array}$$

(57)

Mandel’s Q parameter for the signal mode $Q_{M,s}=2\langle \Delta X_{s1}^2\rangle$ is the same as $-HZ1/\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim \langle \Delta X_{s1}^2\rangle +\langle \Delta X_{s1}\Delta X_{s2}\rangle$ in the $J\to \infty$ limit. If we take the $r\to \infty$ limit in Equation (57), the normalized entanglement criterion converges to $\frac{HZ1}{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }=\frac{1}{4}-\frac{\pi +d}{4\pi (\pi -1)}$, which is the same as the $j\to \infty$ 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/\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle$ is negative: $-\frac{1}{r\pi }-\frac{2d(\pi -1)}{r\pi ({\pi }^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 ${\gamma }_p/{\gamma }_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/{\gamma }_s=40$. Here, to check the far-above threshold entanglement, we performed the positive-P calculation for $G/{\gamma }_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\to \infty$ 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\left(d\right)$ correction makes it easier for an above-threshold NOPO to have a non-classical state.

4.3. Numerical Results

Here, we present the numerical simulation results as a function of normalized excitation p, for ${\gamma }_p/{\gamma }_s=4,G_0/{\gamma }_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{\gamma }_s={10}^2$ to $t{\gamma }_s={10}^4-{10}^5$ depending on the excitation p. For the smallest p (∼5.6), time average was taken from $t{\gamma }_s={10}^2$ to $t{\gamma }_s={10}^5$ because of the slow convergence. For the maximum $p\sim 56$ for WFMC, where the mean signal photon number is $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim 5.0$, we used a photon number space with $M_p=5$ and $M_s=21$ and took the time average from $t{\gamma }_s={10}^2$ to $t{\gamma }_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/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2$. Because of this correction, the entanglement criterion is satisfied, even at the threshold. Next, we present the numerical simulation results for ${\gamma }_p/{\gamma }_s=4,G_0/{\gamma }_s=40,d=1,j=120$. As shown in Figure 4c,d, for the maximum p value (∼32) the mean signal photon number was $\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle \sim 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.

5. 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 ${\chi }^{\left(2\right)}$-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 ${\gamma }_s$. Moreover, the normalized entanglement criterion $HZ1/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^2$ is decreased by a small stimulated emission, while in the second limit, the same value increased. The experimentally required $G/{\gamma }_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/{\langle {\widehat{a}}_s^{†}{\widehat{a}}_s\rangle }^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, $\kappa /{\gamma }_s\sim {10}^{-2}$, has been experimentally confirmed in a second-order nonlinear (${\chi }^{\left(2\right)}$) 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 $\frac{G}{{\gamma }_s}\sim \frac{\hslash c^2{g}_R{Q}_{cs}}{{\epsilon }_r{V}_{cav}}\sim 9\times {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, ${\epsilon }_r$ is the relative dielectric coefficient of the material, and $V_{cav}$ is the volume of the cavity (we used $g_R\sim 57\mathrm{cm}/\mathrm{GW}$ [65], ${\epsilon }_r=12$ and $V_{cav}\sim 0.5{\mathrm{m}}^3$ [62] for silicon photonic crystal nanocavity). The cavity quality factor must be as large as $Q_{cs}\sim {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 (${\chi }^{\left(3\right)}$) nonlinearity [30] can have similar characteristics in terms of below-threshold anti-bunching or above-threshold amplitude squeezing (Appendix G). The CXMs of ${\chi }^{\left(3\right)}$ NOPOs seem to realize entanglement due to these non-classical characteristics.