Entanglement and quantum discord in optically coupled coherent Ising machines

TL;DR

This paper investigates the presence and behavior of entanglement and quantum discord in dissipatively coupled coherent Ising machines, using analytical and numerical methods across different coupling regimes and system sizes.

Contribution

It provides the first comprehensive analysis of entanglement and quantum discord in CIMs, including analytical solutions and numerical simulations for various coupling strengths and system configurations.

Findings

Entanglement disappears when background loss exceeds Ising coupling.

Quantum discord remains finite even with large linear loss.

Entanglement diminishes with distance in a 1D lattice, but discord persists.

Abstract

We present analytical and numerical simulation results for squeezing, entanglement, and quantum discord in a dissipatively coupled coherent Ising machine (CIM). Both analytical solutions and numerical simulation results, which are obtained with positive-P, truncated-Wigner and truncated-Husimi representations for the density operator, predict the presence of entanglement and quantum discord, below and above the threshold of CIM. The entanglement criteria and quantum discord are evaluated as a function of the dissipative coupling strength relative to the background loss. For coupled two DOPOs, while entanglement disappears as the background loss exceeds the Ising coupling strength, the quantum discord remains finite even with a large linear loss. For one-dimensional lattice of DOPOs, while entanglement disappears for DOPO pair with large distance, quantum discord remains finite.

Equations

$$\frac{\partial \hat{\rho }}{\partial t}=\mathcal{L}\hat{\rho }=\epsilon [{\hat{a}}_p^{†}-{\hat{a}}_p,\hat{\rho }]+\frac{\kappa }{2}[{\hat{a}}_s^{†2}{\hat{a}}_p-{\hat{a}}_p^{†}{\hat{a}}_s^2,\hat{\rho }]$$

$$+({\gamma }_p[{\hat{a}}_p,\hat{\rho }{\hat{a}}_p^{†}]+{\gamma }_s[{\hat{a}}_s,\hat{\rho }{\hat{a}}_s^{†}]+h.c.).$$

$$\frac{d{\hat{a}}_p}{dt}$$

$$\frac{d{\hat{a}}_s}{dt}$$

$${\hat{a}}_p=\frac{\epsilon }{{\gamma }_p}-\frac{\kappa }{2{\gamma }_p}{\hat{a}}_s^2+\sqrt{\frac{1}{{\gamma }_p}}{\hat{\xi }}_1.$$

$$\mathcal{L}_{DOPO}\hat{\rho}=\frac{S}{2}[\hat{a}^{\dagger 2}-\hat{a}^{2},\hat{\rho}]+\Bigl{(}\gamma_{s}[\hat{a},\hat{\rho}\hat{a}^{\dagger}]+\frac{B}{2}[\hat{a}^{2},\hat{\rho}\hat{a}^{\dagger 2}]+{\rm h.c.}\Bigr{)}.$$

$$\frac{d\hat{a}}{dt}=-{\gamma }_s\hat{a}+S{\hat{a}}^{†}-B{\hat{a}}^{†}{\hat{a}}^2+\sqrt{2B}{\hat{a}}^{†}{\hat{\xi }}_1+\sqrt{{\gamma }_s}{\hat{\xi }}_2.$$

$${\mathcal{L}}_C\hat{\rho }=J[{\hat{a}}_1-{\hat{a}}_2,\hat{\rho }({\hat{a}}_1^{†}-{\hat{a}}_2^{†})]+h.c..$$

$${\hat{H}}_C=\hslash \chi ({\hat{a}}_1^{†}{\hat{a}}_C+{\hat{a}}_C^{†}{\hat{a}}_1)-\hslash \chi ({\hat{a}}_2^{†}{\hat{a}}_C+{\hat{a}}_C^{†}{\hat{a}}_2).$$

$$\frac{d}{dt}\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_1^{†}\end{array}\right]$$

$$\frac{d}{dt}\left[\begin{array}{c}{\alpha }_2\\ {\alpha }_2^{†}\end{array}\right]$$

$$\frac{d{\alpha }_1}{dt}{|}_C$$

$$\frac{d{\alpha }_2}{dt}{|}_C$$

$${\sigma }_{ij}=\frac{1}{2}⟨[{\hat{R}}_i,{\hat{R}}_j{]}_{+}⟩-⟨{\hat{R}}_i⟩⟨{\hat{R}}_j⟩.$$

$$\sigma =\left[\begin{array}{cc}\underline{\alpha }& \underline{\gamma }\\ \underline{\gamma }& \underline{\alpha }\end{array}\right].$$

$$\left[\begin{array}{cccc}1+j& -p& -j& 0\\ -p& 1+j& 0& -j\\ -j& 0& 1+j& -p\\ 0& -j& -p& 1+j\end{array}\right]\left[\begin{array}{c}⟨\Delta {\alpha }_1^2⟩\\ ⟨\Delta {\alpha }_1^{†}\Delta {\alpha }_1⟩\\ ⟨\Delta {\alpha }_1\Delta {\alpha }_2⟩\\ ⟨\Delta {\alpha }_1^{†}\Delta {\alpha }_2⟩\end{array}\right]=\frac{p}{2}\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].$$

$$\underline{\alpha}={\rm diag}\Bigl{[}1+\frac{p(1-p+j)}{(1-p)(1-p+2j)},1-\frac{p(1+p+j)}{(1+p)(1+p+2j)}\Bigr{]},$$

$$\underline{\gamma}={\rm diag}\Bigl{[}\frac{pj}{(1-p)(1-p+2j)},-\frac{pj}{(1+p)(1+p+2j)}\Bigr{]}.$$

$$\left[\begin{array}{cccc}2p-1+j& -1& -j& 0\\ -1& 2p-1+j& 0& -j\\ -j& 0& 2p-1+j& -1\\ 0& -j& -1& 2p-1+j\end{array}\right]\left[\begin{array}{c}⟨\Delta {\alpha }_1^2⟩\\ ⟨\Delta {\alpha }_1^{†}\Delta {\alpha }_1⟩\\ ⟨\Delta {\alpha }_1\Delta {\alpha }_2⟩\\ ⟨\Delta {\alpha }_1^{†}\Delta {\alpha }_2⟩\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right].$$

$$\underline{\alpha}={\rm diag}\Bigl{[}1+\frac{2p-2+j}{4(p-1)(p-1+j)},1-\frac{2p+j}{4p(p+j)}\Bigr{]},$$

$$\underline{\gamma}={\rm diag}\Bigl{[}\frac{j}{4(p-1)(p-1+j)},-\frac{j}{4p(p+j)}\Bigr{]}.$$

$$\mathfrak{D}=⟨\Delta {\hat{u}}^2⟩+⟨\Delta {\hat{v}}^2⟩,$$

$$\frac{\mathfrak{D}}{2}=\frac{(a_1-c_1)+(a_2+c_2)}{2}.$$

$$\frac{\mathfrak{D}}{2}=1-\frac{p(j-p)}{(1+p)(1-p+2j)}.$$

$$\frac{\mathfrak{D}}{2}=1-\frac{j-1}{4p(p-1+j)}.$$

$${\tilde{\nu }}_{-}^2=1-\frac{p(2j-p)}{(1+p)(1-p+2j)}.$$

$${\tilde{\nu }}_{-}^2=1-\frac{2j-1}{4p(p-1+j)}.$$

$$\mathcal{D}(\sigma )=f(\sqrt{a_1{a}_2})+\underset{\underline{{\sigma }_0}}{\inf }f(\sqrt{\det \underline{\epsilon }})-f({\nu }_{+})-f({\nu }_{-}).$$

$$\underset{\underline{{\sigma }_0}}{\inf }\det \underline{\epsilon }=\frac{a_2}{a_1}(a_1^2-c_1^2).$$

$${\mathcal{L}}_C\hat{\rho }$$

$$\underline{\alpha}=\Bigl{[}1+\frac{p}{1-p+j},1-\frac{p}{1+p+j}\Bigr{]},$$

$$\underline{\alpha}=\Bigl{[}1+\frac{1}{2(p-1)+j},1-\frac{1}{2p+j}\Bigr{]},$$

Full text

I Introduction

A coherent Ising machine (CIM) with delay-line coupled degenerate optical parametric oscillators (DOPOs) has been theoretically studied Wang13 ; Takata15 ; Maruo16 ; Hamerly16 and experimentally demonstratedMarandi14 ; Takata16 ; Inagaki16 . In contrast to a CIM with measurement-feedback-coupled DOPOsMcMahon16 ; Inagaki16b ; Shoji17 ; Yamamura17 , a CIM with optically coupled DOPOs has been shown to satisfy an entanglement criterion Takata15 ; Maruo16 . The theoretical model of optically coupled CIM employs the linear Liouvillian coupling of squeezed states. The entanglement generation from squeezed states and linear Hamiltonian coupling has been considered in the transient caseFurusawa98 ; Furusawa07 ; Paris99 and in the steady state caseOlsen06 . The study of dissipative coupling induced entanglement generation started from atomic systems in the conditional caseCabrillo99 ; Plenio99 and in the steady state caseSchneider02 . The coherent Ising machine realizes dissipative coupling induced steady-state entanglement. The previous numerical results of CIM showed that it satisfies the sufficient condition for entanglement criterionDuan00 , not only below the threshold but also above it Takata15 ; Maruo16 . The degree of entanglement above the threshold depends on various system parameters: for some cases the entanglement criterion is violated before the threshold is reachedTakata15 . The quantum discord Ollivier01 ; Henderson01 was also calculated Takata15 , but the difference in behavior with entanglement has not been clarified. To obtain the condition for entanglement and to understand a different behavior from the quantum discord in a CIM, a comprehensive theoretical decription which passes the numerical test for a wide range of system parameters is required.

The previous theoretical investigation of CIM Takata15 ; Maruo16 depended on numerical simulation, particularly in the phase space method. A rigorous numerical simulation of cavity quantum electrodynamics (C-QED) system can be performed when the field density operator is expanded by Fock states with discrete spectraDalibard92 ; Mu92 ; Molmer93 . The Fock state approach is effective for various nonlinear quantum optical effects Gisin92 if a system is in a small-photon-number regime. In CIMs, however, the system is in a large-photon-number regime where the Fock state approach is poorly suited. Open bosonic quantum systems, such as lasers and optical parametric oscillators (OPOs) are alternatively described by Heisenberg-Langevin equations in the Heisenberg picture Lax66 or by $c$-number SDEs in the Schrödinger picture Lax69 . For lasers with a vastly increasing number of photons, the field density operator can be expanded by using diagonal coherent state expansion Glauber63 ; Lax67 , and the Fokker-Planck equation for the Glauber’s $P$ function is derived. By subsequently using the Ito rule, the $c$-number stochastic differential equation (SDE) can be derived. For a DOPO, however, the Fokker-Planck equation of the $P$ function has a negative diffusion coefficient, and this approach fails.

A CIM with dissipatively coupled DOPOs was first studied by using $c$-number Langevin equations Wang13 , which constitute a $c$-number counterpart to the $q$-number Heisenberg-Langevin equationsLax66 . Later, an equivalence between the $c$-number Heisenberg-Langevin equations and SDEs in Wigner representation was established with the truncation of Fokker-Planck equation in the Wigner representationMaruo16 . The phase space treatments of DOPOs McNeil78 ; Drummond79 ; Drummond81 ; Milburn81 ; Wolinsky88 ; Drummond89 ; Kinsler91 were, on the other hand, performed with positive-$P$ Drummond80 ; Drummond81b or complex-$P$Drummond81b representations. The SDEs in the positive-$P$ representation have been introduced into the numerical simulation of a CIM Takata15 ; Maruo16 . These phase-space methods overcome the difficulties in Glauber’s diagonal $P$ representation and can deal with a varying photon number over many orders of magnitude. Various quantum features can be computed from such SDEs if a sufficient number of trajectories are ensemble averaged. Positive-$P$ and truncated Wigner approaches produce similar values for entanglement criterionMaruo16 , if DOPOs are operating in a so-called weak noise limitWolinsky88 .

On the analytical side, quantum statistical properties of a single DOPO at a steady state have been rigorously obtained by the integration of the Fokker-Planck equation in the complex-$P$ representationDrummond81 ; Milburn81 . In coupled nonlinear quantum optical systems, some recent theoretical studies have used the mean field approximation for the coupling part between subsystemsSheshadri93 ; Boite13 ; Savona17 . Such treatment, motivated by condensed matter theorySheshadri93 , is useful in searching for a macroscopic or global order, but neglects quantum correlation, including entanglement, between subsystems. To elucidate the quantum correlation among constituent DOPOs in a CIM, we consider a weak noise limit Wolinsky88 of a DOPO, instead of the general solution of a DOPO Drummond81 ; Milburn81 . In such a weak noise case, we can analytically calculate entanglement and quantum discord. To check the validity of analytical results, we compute the squeezing, entanglement and quantum discord numerically by using the Positive-$P$, truncated-Wigner and truncated-Husimi SDEs. The numerical simulation results agree completely with independently derived analytical solutions.

The paper is organized as follows. In Sec.II, the quantum master equation for a single DOPO is presented. In Sec.III, we introduce the theoretical model of two dissipatively coupled DOPOs and present the analytical and numerical results on the degrees of squeezing/anti-squeezing, entanglement, and quantum discord. Analytical results are compared with numerical results. Section IV examines a one-dimensional (1D) lattice of DOPOs with only nearest-neighbor coupling. Section V summarizes the main results. Appendix A summarizes the numerical simulation methods based on the positive-$P$, truncated-Wigner and truncated-Husimi SDEs for a single DOPO. Appendix B summarizes the analytical method for a single DOPO. Appendix C shows the traveling wave model for 1D nearest-neighbor-coupled DOPOs.

II Master equation of a DOPO

In this section, we present the theoretical model for a single DOPO. The system of a single DOPO consists of the pump-mode operator $\hat{a}_{p}$ and the signal-mode operator $\hat{a}_{s}$. The quantum master equation of the system is represented as follows:

[TABLE]

Here, $[\hat{A},\hat{B}]=\hat{A}\hat{B}-\hat{B}\hat{A}$. This equation has Hamiltonian and Liouvillian parts. The Hamiltonian part has two terms $\hat{H}_{1}+\hat{H}_{2}$, where $\hat{H}_{1}=i\hbar\varepsilon(\hat{a}_{p}^{\dagger}-\hat{a}_{p})$ is the coherent excitation of the pump mode by an external injection field $\varepsilon$ and $\hat{H}_{2}=i\hbar\frac{\kappa}{2}(\hat{a}_{s}^{\dagger 2}\hat{a}_{p}-\hat{a}_{p}^{\dagger}\hat{a}_{s}^{2})$ is the nonlinear parametric coupling between the signal and pump modes via the second-order nonlinear coupling constant $\kappa$. The Liouvillian parts consist of the dissipation of the pump mode, represented by the half width at half maximum (HWHM) $\gamma_{p}$, and that of the signal mode, represented by $\gamma_{s}$. In the Heisenberg picture, the corresponding Heisenberg-Langevin equations Lax66 are as followsWang13 .

[TABLE]

where $\langle\hat{\xi}_{i}^{\dagger}(t)\hat{\xi}_{j}(t^{\prime})\rangle=0$, and $\langle\hat{\xi}_{i}(t)\hat{\xi}_{j}^{\dagger}(t^{\prime})\rangle=2\delta_{ij}\delta(t-t^{\prime})$. The injected pump field at the oscillation threshold is $\varepsilon_{thr}=\gamma_{p}\gamma_{s}/\kappa$. Assuming that $\gamma_{p}$ is sufficiently large compared to $\gamma_{s}$, the pump mode can be eliminated via

[TABLE]

After eliminating the pump mode adiabatically by substituting Eq.(4) into $\hat{H}_{2}$ in the quantum master equation (II) and averaging over the pump mode and the noise operator $\hat{\xi}_{1}$, we can obtain the quantum master equation for only the signal mode (below, we omit the signal mode subscript $s$ ) Drummond81 :

[TABLE]

Here, the parameter $S=\gamma_{s}\varepsilon/\varepsilon_{thr}$ represents the squeezing/anti-squeezing rate due to the parametric interaction between the pump mode and the signal mode, and $B=\kappa^{2}/(2\gamma_{p})$ represents the degenerate two photon absorption that describes the saturation of the parametric gain. This master equation is equivalent to the following Heisenberg-Langevin equationWang13 :

[TABLE]

The $c$-number counterpart of this equation is equivalent to the SDE in the truncated Wigner representationMaruo16 . The three phase space methods (positive-$P$, truncated Wigner, and truncated Husimi) for dealing with the quantum master equation [Eq.(5)] are presented in Appendix A.

III Two coupled DOPOs

III.1 Model

We now consider a simplest CIM consisting of two DOPOs with ferromagnetic dissipative coupling, as shown in Fig.1. For the two signal modes represented by $\hat{a}_{1}$ and $\hat{a}_{2}$, the Liouvillian of an entire system is the sum of a single DOPO part and coupling part: $\frac{\partial\hat{\rho}}{\partial t}=\sum_{r}\mathcal{L}_{DOPO}^{(r)}\hat{\rho}+\mathcal{L}_{C}\hat{\rho}$. Here, $\mathcal{L}_{DOPO}^{(r)}$ operates on only $r$-th DOPO and $\mathcal{L}_{C}\hat{\rho}$ operates on both. We consider the following Liouvillian which represents the dissipative ferromagnetic couplingTakata15 :

[TABLE]

This Liouvillian is derived when the two signal modes $\hat{a}_{1}$ and $\hat{a}_{2}$ with eigenfunctions $E_{s1}$ and $E_{s2}$, respectively, are coupled via the standing-wave cavity mode $\hat{a}_{C}$ with eigenfunction $E_{C}$, which has a node at the center of the cavity and has the oppositely displaced anti-nodes at the interface of $\hat{a}_{1}$ and $\hat{a}_{2}$:

[TABLE]

Here, the coupling mode $\hat{a}_{C}$ has a dissipation rate (HWHM) $\gamma_{C}$. When the two signal modes $\hat{a}_{1}$ and $\hat{a}_{2}$ are excited with the same phase (ferromagnetic order), the dissipation induced by $\gamma_{C}$ does not occur. When $\gamma_{C}$ is large, mode $\hat{a}_{C}$ can be adiabatically eliminated via the Lindblad procedure, and then equation (7) is obtained with $J\sim\chi^{2}/\gamma_{C}$.

The positive-$P$ SDE for two coupled DOPOs is obtained, as the extension of that for a single DOPO, shown in Appendix A, as

[TABLE]

[TABLE]

Here $\alpha_{r}$ and $\alpha_{r}^{\dagger}(r=1,2)$ are independent $c$-numbers satisfying $\langle\alpha_{r}\rangle=\langle\alpha_{r}^{\dagger}\rangle^{*}$ and $\xi_{Rr}$ and $\xi_{Rr}^{\dagger}$ are independent real number random variables satisfying $\langle\xi_{Rr}(t)\xi_{Rr^{\prime}}(t^{\prime})\rangle=\delta_{rr^{\prime}}\delta(t-t^{\prime})$ and $\langle\xi_{Rr}^{\dagger}(t)\xi_{Rr^{\prime}}^{\dagger}(t^{\prime})\rangle=\delta_{rr^{\prime}}\delta(t-t^{\prime})$. When the system has the ferromagnetic order $\alpha_{1}\sim\alpha_{2}$, each cavity has no additional loss due to coupling. Therefore, the oscillation occurs with the same excitation as for a single DOPO i.e., at normalized pump rate $p=S/\gamma_{s}=1$. The dissipative coupling Liouvillian (7) does not introduce new noise terms into the positive-$P$ SDE Takata15 . In the Wigner and Husimi SDEs, however, the dissipative coupling introduces new noise terms. When the coupling parts of the SDE are written as $\frac{d\alpha_{1}}{dt}|_{C}$ and $\frac{d\alpha_{2}}{dt}|_{C}$, they are represented as:

[TABLE]

Here, $A=1$ for the Husimi representation and $A=\frac{1}{2}$ for the Wigner representation. The complex noise source $\xi_{C}$, satisfying $\langle\xi^{*}_{C}(t)\xi_{C}(t^{\prime})\rangle=2\delta(t-t^{\prime})$, is common for $\alpha_{1}$ and $\alpha_{2}$, since it comes from the vacuum noise of the coupling mode $\hat{a}_{C}$.

III.2 Covariance matrix

We consider the steady-state covariance matrix for two coupled DOPOs by applying the fluctuation analysis in the positive-$P$ representation, that for a single DOPO is presented in Appendix B. If two modes have the same excitations, then they have symmetry under the exchange of modes. We assume that the fluctuation products of the positive-$P$ amplitudes satisfy $\langle\Delta\alpha_{1}^{2}\rangle=\langle\Delta\alpha_{2}^{2}\rangle$ and $\langle\Delta\alpha_{1}^{\dagger}\Delta\alpha_{1}\rangle=\langle\Delta\alpha_{2}^{\dagger}\Delta\alpha_{2}\rangle$. The covariance matrix $\sigma$ is defined for the vectors $\overrightarrow{\hat{R}}=\sqrt{2}[\hat{X}_{1},\hat{P}_{1},\hat{X}_{2},\hat{P}_{2}]$,

[TABLE]

Here, $[\hat{A},\hat{B}]_{+}=\hat{A}\hat{B}+\hat{B}\hat{A}$. For two DOPOs with dissipative coupling, this matrix can be written as follows:

[TABLE]

Here, $\underline{\alpha}={\rm diag}(a_{1},a_{2})$ and $\underline{\gamma}={\rm diag}(c_{1},c_{2})$. $a_{1}=2\langle\Delta\hat{X}_{1}^{2}\rangle=2\langle\Delta\hat{X}_{2}^{2}\rangle$, and $a_{2}=2\langle\Delta\hat{P}_{1}^{2}\rangle=2\langle\Delta\hat{P}_{2}^{2}\rangle$ represent the single-mode variances of $\hat{X}$ and $\hat{P}$, respectively. On the other hand, $c_{1}=2\langle\Delta\hat{X}_{1}\Delta\hat{X}_{2}\rangle$, and $c_{2}=2\langle\Delta\hat{P}_{1}\Delta\hat{P}_{2}\rangle$ represent the intermode correlation of $\hat{X}$ and $\hat{P}$, respectively.

Below the oscillation threshold ($p<1$), where $\langle\alpha_{1}\rangle=\langle\alpha_{2}\rangle=0$, the steady-state amplitude fluctuations satisfy the following:

[TABLE]

Here, $j:=J/\gamma_{s}$ is the normalized coupling coefficient. $\langle\Delta\alpha_{1}^{2}\rangle$ and $\langle\Delta\alpha_{1}^{\dagger}\Delta\alpha_{1}\rangle$ represent the averaged fluctuation products for the same DOPO, while $\langle\Delta\alpha_{1}\Delta\alpha_{2}\rangle$ and $\langle\Delta\alpha^{\dagger}_{1}\Delta\alpha_{2}\rangle$ represent the averaged fluctuation products for different DOPOs. For the system of two DOPOs below the threshold, all the nonzero components of the covariance matrix are obtained as follows.

[TABLE]

[TABLE]

Above the threshold, the mean amplitudes are given by $\langle\alpha_{1}\rangle=\langle\alpha_{2}\rangle=\sqrt{\frac{\gamma_{s}}{B}(p-1)}$. We consider the fluctuation around this mean value using the linearization technique Chaturvedi77 . This procedure in a single DOPO is shown in Appendix B. Under a steady-state condition, the fluctuation amplitude products are represented as follows:

[TABLE]

These amplitude fluctuation products produce the following nonzero components of covariance matrix:

[TABLE]

[TABLE]

III.3 Sufficient condition for entanglement

The sufficient criterion for entanglement is given in Ref.Duan00 , using

[TABLE]

where $\hat{u}=\hat{X}_{1}-\hat{X}_{2}$, $\hat{v}=\hat{P}_{1}+\hat{P}_{2}$, $\hat{X}_{r}=\frac{\hat{a}_{r}+\hat{a}_{r}^{\dagger}}{\sqrt{2}}(r=1,2)$ and $\hat{P}_{r}=\frac{\hat{a}_{r}-\hat{a}_{r}^{\dagger}}{\sqrt{2}i}(r=1,2)$. If $\mathfrak{D}/2<1$, then the system’s density matrix is not separated into the product states of the individual DOPOs : $\hat{\rho}\neq\hat{\rho}_{1}\otimes\hat{\rho}_{2}$. This criterion can be calculated using the components of covariance matrix as

[TABLE]

From Eqs.(16)(17), below threshold $\mathfrak{D}/2$ is represented as follows:

[TABLE]

Then the sufficient condition for entanglement is satisfied if $j>p$ is satisfied. This means that when $j<1$, the entanglement criterion is violated before the threshold is reached as seen in Ref.Takata15 . Above the threshold, the entanglement criterion is represented as:

[TABLE]

We can thus see that the entanglement criterion is satisfied even above the threshold, when the dissipative coupling rate is larger than a linear loss rate $j=J/\gamma_{s}>1$.

III.4 Necessary and sufficient entanglement criterion

We now consider the necessary and sufficient entanglement criterion following Simon’s theorySimon00 ; Vidal02 ; Serafini03 ; Adesso04 . For symmetric two boson gaussian system, this is easily shown to be equivalent to ”Theorem 2” in Ref.Duan00 . In the description in Ref.Adesso04 , the necessary and sufficinet entanglement criterion is described by the symplectic eigenvalues of covariance matrix $\sigma_{ij}$. The covariance matrix can be diagonalized with a symplectic transformation $\Sigma$ as $\Sigma\sigma\Sigma^{T}={\rm diag}(\nu_{+},\nu_{+},\nu_{-},\nu_{-})$, where $\nu_{+}>\nu_{-}$. The symplectic eigenvalues are related to the uncertainty relation, and a physical state must satisfy $\nu_{-}\geq 1$ Adesso04 . These eigenvalues $\nu_{+}$ and $\nu_{-}$ are obtained as the eigenvalues of $i\Omega\sigma$ where $\Omega=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}\oplus\begin{bmatrix}0&1\\ -1&0\end{bmatrix}$. Simon’s entanglement criterion is defined on the partially transposed covariance matrix $\bar{\sigma}=\Lambda\sigma\Lambda^{T}$, where $\Lambda={\rm diag}(1,1,1,-1)$. The symplectic eigenvalues of $\bar{\sigma}$ are eigenvalues of $i\Omega\bar{\sigma}$ and satisfy $\Sigma^{\prime}\bar{\sigma}\Sigma^{\prime T}={\rm diag}(\tilde{\nu}_{+},\tilde{\nu}_{+},\tilde{\nu}_{-},\tilde{\nu}_{-})$ ($\tilde{\nu}_{+}>\tilde{\nu}_{-}$), for a symplectic transformation $\Sigma^{\prime}$. The separable state must satisfy $\tilde{\nu}_{-}\geq 1$Adesso04 .

The above procedures produce the symplectic eigenvalues of $\sigma$. The two independent symplectic eigenvalues of $\sigma$ are calculated as eigenvalues $\lambda$ of $i\Omega\sigma$, and are reduced to $\lambda^{2}=(a_{1}+c_{1})(a_{2}+c_{2}),(a_{1}-c_{1})(a_{2}-c_{2})$. Substituting Eqs.(16)(17) into $\lambda^{2}$, we obtain $\nu_{+}^{2}=1+\frac{p^{2}}{1-p^{2}}$ and $\nu_{-}^{2}=1+\frac{p^{2}}{(1+p+2j)(1-p+2j)}$. To obtain necessary and sufficient entanglement criterion, two independent symplectic eigenvalues of $\bar{\sigma}$ are calculated as eigenvalues $\lambda^{\prime}$ of $i\Omega\bar{\sigma}$, which are reduced to $\lambda^{\prime 2}=(a_{1}+c_{1})(a_{2}-c_{2}),(a_{1}-c_{1})(a_{2}+c_{2})$. The smallest symplectic eigenvalue of $\bar{\sigma}$ is

[TABLE]

When $\tilde{\nu}_{-}<1$, the two coupled DOPOs are inseparable. This entanglement condition for below the threshold is satisfied when $j>p/2$.

We next consider the above-threshold characteristics described by the covariance matrix $\sigma$. The symplectic eigenvalues of $\sigma$ are $\nu_{+}^{2}=1+\frac{1}{4p(p-1)}$ and $\nu_{-}^{2}=1+\frac{1}{4(p+j)(p-1+j)}$. The smallest symplectic eigenvalue of $\bar{\sigma}$ is represented as

[TABLE]

Therefore, the necessary and sufficient entanglement criterion for above the threshold is satisfied for $j>1/2$. The required $j$ to satisfy the entanglement criterion is thus smaller than that for sufficient criterion obtained in the previous section ($j>1$). The symplectic eigenvalue of partially transposed covariance matrix $\bar{\sigma}$ represents the degree of entanglement and the minimum value of $\tilde{\nu}_{-}$ is $\tilde{\nu}_{-}\rightarrow\frac{1}{\sqrt{2}}$ when $p=1$ and $j\rightarrow\infty$. This is equivalent to that obtained for one-dimensional model in Ref.Yanagimoto19 .

III.5 Quantum discord

Next, we consider the quantum discord of two dissipatively coupled DOPOs. When the total system density operator $\hat{\rho}$ is composed of the first DOPO with $\hat{\rho}_{1}$ and the second DOPO with $\hat{\rho}_{2}$, the von Neumann entropy of the total system is defined by $S(\hat{\rho})=-{\rm Tr}\hat{\rho}\log\hat{\rho}$. The total entropy also seems to be written as the summation $S(\hat{\rho}_{1})+S(\hat{\rho}_{2|1})$, where $S(\hat{\rho}_{1})$ is the reduced entropy of the first DOPO, while $S(\hat{\rho}_{2|1})$ is the conditional (or residual) entropy of the second DOPO which is given by the subtraction of mutual information from $S(\hat{\rho}_{2})$. In classical information theory, these two formula for the total entropy are identical, but when the systems have quantumness they generally differOllivier01 ; Henderson01 . The quantum discord is defined as the difference between these two expressions for entropies. Non-zero quantum discord is identified as a metric for quantifying the quantum correlation of mixed states and as a useful quantum resource in a specific quantum computational model Datta07 ; Lanyon08 .

For a system of two DOPOs, using the theory of two-mode gaussian states Giorda10 ; Adesso10 , quantum discord for the covariance matrix $\sigma$ is written as follows:

[TABLE]

Here $f(x)=\frac{x+1}{2}\log\frac{x+1}{2}-\frac{x-1}{2}\log\frac{x-1}{2}$, and $\underline{\varepsilon}=\underline{\alpha}-\underline{\gamma}(\underline{\sigma_{0}}+\underline{\alpha})^{-1}\underline{\gamma}^{T}$ where $\underline{\sigma_{0}}$ relates to the measurement which removes the mutual informationEisert02 ; Giedke02 . The first two terms on the right hand side of Eq.(27) represent the entropies $S(\hat{\rho}_{1})$ and $S(\hat{\rho}_{2|1})$. The last two terms represent the entropy $S(\hat{\rho})$ consisting of two DOPOsHolevo01 . In the calculation of $S(\hat{\rho}_{2|1})$, we consider $\theta=0$ case of Ref.Adesso10 where $\underline{\sigma_{0}}={\rm diag}(\lambda,\lambda^{-1})$($\lambda>0$). When $(a_{2}c_{1}^{2}-a_{1}c_{2}^{2}(a_{1}^{2}-c_{1}^{2}))(a_{2}c_{1}^{2}(a_{2}^{2}-c_{2}^{2})-a_{1}c_{2}^{2})<0$, $\frac{d\det\underline{\varepsilon}(\lambda)}{d\lambda}=0$ can be obtained for positive $\lambda$. When $(a_{2}c_{1}^{2}-a_{1}c_{2}^{2}(a_{1}^{2}-c_{1}^{2}))(a_{2}c_{1}^{2}(a_{2}^{2}-c_{2}^{2})-a_{1}c_{2}^{2})\geq 0$, $\frac{d\det\underline{\varepsilon}(\lambda)}{d\lambda}=0$ is satisfied for only negative $\lambda$, and the minimum value of the determinant of the matrix $\underline{\varepsilon}$ is obtained at $\lambda\rightarrow 0$ Adesso10 as

[TABLE]

In our two-DOPO system, the residual entropy is minimized for $\lambda\rightarrow 0$. This is shown from the analytical result of $D:=a_{2}c_{1}^{2}(a_{2}^{2}-c_{2}^{2})-a_{1}c_{2}^{2}(a_{1}^{2}-c_{1}^{2})$, which is smaller than both $a_{2}c_{1}^{2}-a_{1}c_{2}^{2}(a_{1}^{2}-c_{1}^{2})$, and $a_{2}c_{1}^{2}(a_{2}^{2}-c_{2}^{2})-a_{1}c_{2}^{2}$, because of $a_{1}^{2}-c_{1}^{2}>1$ and $a_{2}^{2}-c_{2}^{2}<1$. Both $D=\frac{2p^{3}j^{2}(1+j)(1+2j)}{(1+p)^{2}(1-p)^{2}(1-p+2j)^{2}(1+p+2j)^{2}}$ below the threshold and $D=\frac{j^{2}(2p-1)(2p-1+j)(2p-1+2j)}{128p^{2}(p-1)^{2}(p+j-1)^{2}(p+j)^{2}}$ above the threshold have positive values. Therefore $(a_{2}c_{1}^{2}-a_{1}c_{2}^{2}(a_{1}^{2}-c_{1}^{2}))(a_{2}c_{1}^{2}(a_{2}^{2}-c_{2}^{2})-a_{1}c_{2}^{2})\geq 0$ is always satisfied. The quantum discord is always nonzero in an optically coupled CIM, no matter how large the dissipation $\gamma_{s}$ is compared to the mutual coupling $J$, because of the nonzero $c_{1}$ and $c_{2}$ values. The quantum discord at the threshold of two DOPOs is calculated with the approximation $f(x)\rightarrow 1+\log(x/2)(x\rightarrow\infty)$. In the limit $p\rightarrow 1$ and $j\rightarrow\infty$, we have the maximum quantum discord $\mathcal{D}\sim\frac{1}{2}\log\frac{3}{4}+f(\sqrt{\frac{3}{2}})\sim 0.22$.

III.6 Numerical simulation

To test the validity of analytical results for a ferromagnetically coupled two-DOPO system, we performed a numerical simulation based on the positive-$P$, truncated Wigner, and truncated Husimi SDEs. The simulation was performed with a small two-photon absorption loss $B/\gamma_{s}=10^{-4}$ and coupling coefficient $j=7/3$. The fluctuations were calculated with the time average for a single trajectory, assuming ergodicity. The total simulation time was a period of $10^{6}/\gamma_{s}$. In the first time period of $t_{f}=10^{4}/\gamma_{s}$, the time average was not taken, and the pump excitation varies on time via $p(t)=p\sqrt{t/t_{f}}$. After the first period of $t_{f}$, the excitation was held constant to $p$, and the time average was taken. The time step was $\Delta t=10^{-3}/\gamma_{s}$, $\Delta t=5\times 10^{-5}/\gamma_{s}$, and $\Delta t=2.5\times 10^{-5}/\gamma_{s}$, respectively for positive-$P$, truncated Wigner and truncated Husimi SDEs.

We first present the results for $\langle\Delta\hat{X}^{2}\rangle=a_{1}/2$ and $\langle\Delta\hat{P}^{2}\rangle=a_{2}/2$ for the two-DOPO system in Fig.2(a). The circles with different colors represent the numerical results from positive-$P$, truncated-Wigner (T-Wigner), and truncated-Husimi (T-Husimi) calculations. The filled circles represent $\langle\Delta\hat{X}^{2}\rangle$, while the open circles represent $\langle\Delta\hat{P}^{2}\rangle$. The analytical results are shown by black solid lines. The fluctuation of vacuum is shown by a gray dashed line. Far below the threshold $(p\ll 1)$ and far above the threshold $(p\gg 1)$, the DOPO field should be in a vacuum state and a coherent state, respectively, so that the variance $\langle\Delta\hat{X}^{2}\rangle$ and $\langle\Delta\hat{P}^{2}\rangle$ are expected to approach an asymptotic value $\langle\Delta\hat{X}^{2}\rangle=\langle\Delta\hat{P}^{2}\rangle=\frac{1}{2}$. On the other hand, the degree of squeezing and anti-squeezing should be maximum at the threshold $(p=1)$. The numerical simulations reproduce these theoretical predictions well and thus verified the analytical results presented in the previous sections. The correlations of canonical coordinates are shown in Fig.2(b). Mean fluctuation product of canonical coordinates $\langle\Delta\hat{X}_{1}\Delta\hat{X}_{2}\rangle$ has a positive correlation and has a singular increase at the threshold, similarly to the $\langle\Delta\hat{X}^{2}\rangle$. The fluctuation product of canonical momenta $\langle\Delta\hat{P}_{1}\Delta\hat{P}_{2}\rangle$ has a negative correlation and remains finite at the threshold.

Duan-Giedke-Cirac-Zoller sufficient criterion for entanglement ($\mathfrak{D}/2$) is shown in Fig.2(c). When this value is smaller than one, the system is inseparable and entangled. We see that the entanglement is maximum at the threshold, but exists even above the threshold. Note that the entanglement condition above the threshold ($j>1$) is satisfied in this numerical simulation. The theoretical values shown by black line were calculated from Eqs.(23) and (24), and these analytical results match the numerical simulation results. The Duan-Simon’s necessary and sufficient entanglement criterion is shown in Fig.2(d), which also shows a maximum entanglement at the threshold. Next, Fig.2(e) shows the quantum discord as a function of $p$, calculated from the covariance matrix $\sigma_{ij}$ Takata15 . The analytical results were obtained from Eqs.(27) and (28) and are shown by a black line. As shown by filled circles, the quantum discord takes a maximum value at the threshold and remains finite even at the pump rate far below and far above it. Figure 2(f) shows the analytical quantum discord of dissipatively coupled DOPOs at three different pump rates as a function of normalized dissipative coupling $j=J/\gamma_{s}$. A dashed vertical line represents the necessary and sufficient entanglement criterion for above threshold ($j=1/2$). For a large loss limit ($j\ll 1$), the entanglement disappears, but the quantum discord survives particularly around the threshold: $p\sim 1$.

III.7 Mean field approximation

For comparison, we present here a theoretical model for mean-field coupled CIMs. A mean-field coupling theory in nonlinear quantum optical system has been developed in Bose-Hubbard like models Boite13 ; Savona17 . In such an approach, the Hamiltonian coupling between different sites, $\hat{H}_{C}=\hbar t\sum_{\langle ij\rangle}(\hat{a}_{i}^{\dagger}\hat{a}_{j}+\hat{a}_{j}^{\dagger}\hat{a}_{i})$, is approximated by the mean-field via $\hat{a}_{i}^{\dagger}\hat{a}_{j}\rightarrow\hat{a}_{i}^{\dagger}\langle\alpha_{j}\rangle+\langle\alpha_{i}\rangle^{*}\hat{a}_{j}$. Here, $t$ is a hopping element. Similarly, we can consider the mean-field coupling in two dissipatively coupled DOPOs, where the Liouvillian coupling [Eq.(7)] is replaced by $[\hat{a}_{1}-\hat{a}_{2},\hat{\rho}(\hat{a}_{1}^{\dagger}-\hat{a}_{2}^{\dagger})]\rightarrow[\hat{a}_{1}-\langle\alpha_{2}\rangle,\hat{\rho}(\hat{a}_{1}^{\dagger}-\langle\alpha_{2}\rangle^{*})]+[\hat{a}_{2}-\langle\alpha_{1}\rangle,\hat{\rho}(\hat{a}_{2}^{\dagger}-\langle\alpha_{1}\rangle^{*})]$. Assuming that $\langle\alpha_{r}\rangle(r=1,2)$ are real, $\mathcal{L}_{C}\hat{\rho}$ in Eq.(7) is replaced as

[TABLE]

Such a model leads to a covariance matrix Eq.(14) replaced by

[TABLE]

below the threshold and by

[TABLE]

above the threshold. $c_{1}=c_{2}=0$ is satisfied both below and above the threshold in the mean field approximation. Far below the threshold ($p\ll 1$) and far above the threshold ($p\gg 1$), these covariance matrices are identical to those without mean field approximation [Eqs.(16)(17)(19) and (20)], in the small coupling limit ($j\ll 1$). However, the difference emerges near the threshold ($p=1$). In the mean field approximation, $\langle\Delta\hat{X}_{1}^{2}\rangle$ is finite and $\langle\Delta\hat{X}_{1}\Delta\hat{X}_{2}\rangle$ is zero at the threshold, although they diverge without mean field approximation. The mean field coupled model has zero quantum discord, for $\inf_{\underline{\sigma_{0}}}\det\underline{\varepsilon}=\nu_{+}^{2}=\nu_{-}^{2}=a_{1}a_{2}$.

Mean field characteristics is calculated by the analytical results of the single DOPODrummond81 ; Milburn81 and self-consistent loopBoite13 . The mean amplitude is

[TABLE]

where $c=\sqrt{\frac{S}{B}}$, $x=\frac{\gamma_{s}+J}{B}$, and $e=\frac{J\langle\hat{a}\rangle}{\sqrt{SB}}$. The generalized moment is obtained by

[TABLE]

where $k^{\prime}=k+m$ and $k^{\prime\prime}=k+n$. We present the characteristics of mean-field coupled coherent Ising machines via self-consistent equation. In Fig.3(a), the bifurcation diagram is presented. When we consider the initial mean amplitude $\langle\hat{a}\rangle>0$ (shown with filled circles in Fig.3(a)), from the self-consistent loop [Eq.(32)], the mean amplitude is expected to converge onto $\langle\hat{a}\rangle=0$ below the threshold, and $\langle\hat{a}\rangle\sim\sqrt{\frac{\gamma_{s}}{B}(p-1)}$ above the threshold, shown with open circles in Fig.3(a). With the self-consistent loop, non-gaussian correction is precisely considered. We considered the results with $B/\gamma_{s}=0.019$ and $j=7/3$, starting from $\langle\hat{a}\rangle=1$. In Fig.3(b), we present the mean amplitude depending on the count of iteration in the self-consistent loop. For $p\leq 1$, the mean amplitude converges into zero, while for $p>1$ the mean amplitude converges into finite values. In Fig.3(c), the mean amplitudes after $2000$ iterations were presented as a function of $p$. The mean amplitudes above the threshold are described by $\sqrt{\frac{\gamma_{s}}{B}(p-1)}$ shown by black line. The fluctuations are shown as circles in Fig.3(d). These fluctuations follow the black lines given by Eqs.(30)(31).

Finally, we present the fluctuations by finite-particle simulation using positive-$P$ SDE. When mean field is constructed from $N_{p}$ positive-$P$ amplitudes, below the threshold steady-state fluctuations are

[TABLE]

Above the threshold, steady state fluctuations of $\hat{X}$ are represented by Eqs.(III.7)(35) with substitution $p\rightarrow 1$ in the numerator and $1-p\rightarrow 2(p-1)$ in the denominator, and those of $\hat{P}$ are represented by Eqs.(III.7)(37) with substitution $p\rightarrow 1$ in the numerator and $1+p\rightarrow 2p$ in the denominator. These converge into the values of gaussian theory [Eqs.(30)(31)] in the $N_{p}\rightarrow\infty$ limit. Fig.3(e)(f) show the positive-$P$ numerical results for $B/\gamma_{s}=10^{-4}$ and $j=7/3$, using $N_{p}=5$ and $9\times 10^{4}$ particles, respectively. Numerical simulation was performed in the similar way as Fig.2, where we considered $t_{f}=10^{3}/\gamma_{s}$ and total time development with period $2\times 10^{3}/\gamma_{s}$. For $N_{p}=5$, numerical results are well fitted by Eqs.(III.7)(III.7), considering the correction of finiteness of particle number $N_{p}$. For $N_{p}=9\times 10^{4}$, numerical results are fitted well by Eqs.(30)(31). At $p=1$, the divergence of $\langle\Delta\hat{X}^{2}\rangle$ remaining in Eqs.(III.7)(III.7) cannot be observed in Fig.3(f) due to large $N_{p}$ and finite $\gamma_{s}/B=10^{4}$.

IV 1D lattice of $N$ DOPOs

IV.1 Model

We consider a CIM consisting of $N$-DOPOs ($\hat{a}_{r}(r=1,\cdots,N)$) which is the ferromagnetic 1D lattice with nearest neighbor coupling. The dissipative coupling Liouvillian is obtained from the extension of Eq.(7) as follows:

[TABLE]

Here, we consider a periodic system such that $\hat{a}_{N+1}$ is identical to $\hat{a}_{1}$. Such a 1D lattice can be constructed with $N$ dissipative standing-wave modes $\hat{a}_{C,r}(r=1,\cdots,N)$ with half-width $\gamma_{C}$ :

[TABLE]

Here, the $r$-th DOPO is coupled to the $(r-1)$-th DOPO via the dissipative mode $\hat{a}_{C,r-1}$, and is also coupled to the $(r+1)$-th DOPO via the dissipative mode $\hat{a}_{C,r}$. When all dissipative modes $\hat{a}_{C,r}$ have the same loss rate $\gamma_{C}$ with $\gamma_{C}\gg\gamma_{s}$, the coupling coefficient of the Liouvillian follows $J/2\sim\chi^{2}/\gamma_{C}$. In the Wigner and Husimi SDE, from Eq.(38), the coupling part in the SDE for the $r$-th DOPO is represented as follows:

[TABLE]

Here, $A=1$ for the Husimi SDE and $A=\frac{1}{2}$ for the Wigner SDE. The noise sources satisfy $\langle\xi_{C,r}^{*}(t)\xi_{C,r^{\prime}}(t^{\prime})\rangle=2\delta_{r,r^{\prime}}\delta(t-t^{\prime})$. This model was obtained as the extension of the standing wave model in Fig.1. The theory of the traveling pulse model of a CIM Marandi14 in the Wigner representation is shown in Appendix C.

IV.2 Extended Duan-Giedke-Cirac-Zoller criterion

First, we take the analytical approach using the positive-$P$ SDE. For the $r$-th signal mode $\hat{a}_{r}(r=1,\cdots,N)$, the positive-$P$ SDE is written as:

[TABLE]

We introduce the Fourier-transformed fluctuation amplitude $\Delta\tilde{\alpha}_{k}$:

[TABLE]

Here, $\theta_{k}=\frac{2\pi}{N}k$. Below the threshold, we obtain the following SDE for the Fourier components:

[TABLE]

The Fourier-transformed random noise sources are defined as

[TABLE]

and satisfy $\langle\tilde{\xi}_{R,k}(t)\tilde{\xi}_{R,-k^{\prime}}(t^{\prime})\rangle=\delta_{k,k^{\prime}}\delta(t-t^{\prime})$. We obtain the steady-state fluctuation of Fourier components as follows.

[TABLE]

We now calculate the extended Duan-Giedke-Cirac-Zoller criterion for even number $N$-DOPOs Maruo16 , denoted as $\mathfrak{D}^{\prime}=\langle\Delta\hat{u}^{2}\rangle+\langle\Delta\hat{v}^{2}\rangle$, where $\hat{u}=\sum_{r=1}^{N}(-1)^{r+1}\hat{X}_{r}$ and $\hat{v}=\sum_{r=1}^{N}\hat{P}_{r}$. When $\mathfrak{D}^{\prime}/N<1$, the entire system is not separated into the product state, i.e. $\hat{\rho}\neq\otimes_{r=1}^{N}\hat{\rho}_{r}$. The extended criterion can be written with Fourier components as

[TABLE]

We can see that this expression is reduced to $\frac{\mathfrak{D}^{\prime}}{N}=1-\frac{p(j-p)}{(1+p)(1-p+2j)}$. Note that this $\frac{\mathfrak{D}^{\prime}}{N}$ contains the case of $N=2$ as a special case (see Eq.(23)).

Above the threshold, we obtain the following SDE for the Fourier components:

[TABLE]

Here, $S-B\langle\alpha\rangle^{2}\sim\gamma_{s}$ is satisfied above the threshold. A similar equation was obtained in the dissipative Bose Hubbard model with the nonlinear Kerr effectBoite14 . From Eq.(IV.2), we obtain the steady-state fluctuation correlations as follows.

[TABLE]

The extended Duan-Giedke-Cirac-Zoller criterion $\mathfrak{D}^{\prime}/N$ is given by Eq.(IV.2), and the result is reduced to $\frac{\mathfrak{D}^{\prime}}{N}=1-\frac{j-1}{4p(p-1+j)}$, which is identical to $\mathfrak{D}/2$ [Eq.(24)] for two DOPOs. Therefore, $j=1$ is required to satisfy the extended entanglement criterion above the threshold.

To understand these results, we introduce an intuitive picture based on the loss spectrum model. After linearizaton, we consider the loss matrix $\underline{M_{k}}$ for fluctuation defined as

[TABLE]

where $\overrightarrow{A_{k}}:=\begin{bmatrix}\Delta\tilde{\alpha}_{k}\\ \Delta\tilde{\alpha}_{-k}^{\dagger}\end{bmatrix}$. The matrix $\underline{M_{k}}$ in a 1D DOPO system with nearest-neighbor coupling is described by

[TABLE]

below the threshold, and by

[TABLE]

above the threshold. When the loss matrix is diagonalized, the $\hat{X}$-like ($\Delta\tilde{\alpha}_{k}+\Delta\tilde{\alpha}_{-k}^{\dagger}$) mode has a normalized loss represented by $\Gamma_{X}(k)/\gamma_{s}=1-p+j(1-\cos\theta_{k})$ below threshold and $\Gamma_{X}(k)/\gamma_{s}=2(p-1)+j(1-\cos\theta_{k})$ above threshold. Likewise, the $\hat{P}$-like ($\Delta\tilde{\alpha_{k}}-\Delta\tilde{\alpha}_{-k}^{\dagger}$) mode has a normalized loss spectrum of $\Gamma_{P}(k)/\gamma_{s}=1+p+j(1-\cos\theta_{k})$ below threshold and $\Gamma_{P}(k)/\gamma_{s}=2p+j(1-\cos\theta_{k})$ above threshold. Figure 4(a) shows these spectra with $j=7/3$ for three excitation strengths $p=0.3$, $p=1$, and $p=3$. Here, the filled circles represent $\hat{X}$-like modes, while the open circles represent $\hat{P}$-like modes. As the fluctuation term does not depend on $\theta_{k}$ and the mode type (i.e., $\hat{X}$-like or $\hat{P}$-like ), the fluctuation magnitude is represented by the loss function. When the loss is smaller, the fluctuation is larger in proportion to $\Gamma(k)^{-1}$. We can see that when $p=1$, the $\hat{X}$-like mode becomes lossless at $k=0$ and lasing occurs at $k=0$. The extended entanglement criterion for the entire lattice is represented by fluctuations at two symmetric points in wavenumber space [Eq.(IV.2)]. The extended entanglement criterion ($\mathfrak{D}^{\prime}/N$) is satisfied when the loss of the $\hat{X}$ -like band at $k=\frac{N}{2}$ is larger than that of the $\hat{P}$-like band at $k=0$. This is always the case far below the threshold. Above the threshold, however, the criterion is satisfied only when $j>1$. In Fig.4(b), we present the wavenumber dependent loss for $j=2/3$. In this case, the extended entanglement criterion is satisfied for $p=0.3$, but not satisfied for $p=1$ and $p=3$.

We performed a numerical simulation to confirm the analytical results for the 1D DOPO ring with nearest-neighbor coupling. The simulation was performed in the same way as Fig.2(a-e), with a small two-photon absorption loss $B/\gamma_{s}=10^{-4}$. We considered $N=32$ DOPOs with $j=7/3$ or $j=2/3$. Figs.5(a)(b) show the extended entanglement criterion for an entire lattice. As seen in the loss spectrum of Fig.4, the extended entanglement criterion is always satisfied for $j=7/3$. However, for $j=2/3$ it is satisfied only for $p<2/3$.

IV.3 Entanglement criterion between nearest-neighbor DOPOs

Instead of the entire lattice, we can consider the entanglement for two DOPOs, after tracing out the remaining $(N-2)$ DOPOs. Duan-Giedke-Cirac-Zoller sufficient entanglement criterion for a pair with site distance $r$ is calculated as follows: $\frac{\mathfrak{D}(r)}{2}=1+\frac{2}{N}\sum_{k}(\langle\Delta\tilde{\alpha}^{\dagger}_{k}\Delta\tilde{\alpha}_{k}\rangle-\langle\Delta\tilde{\alpha}_{-k}\Delta\tilde{\alpha}_{k}\rangle\cos(r\theta_{k}))$. We first consider the $r=1$ case (nearest neighbor pair). In the limit $N\rightarrow\infty$, we compute the above summation by replacing it with the integration $\frac{1}{N}\sum_{k=0}^{N-1}f(\theta_{k})\rightarrow\int_{-\pi}^{+\pi}f(\theta)\frac{d\theta}{2\pi}$. After the integration, we obtain

[TABLE]

below the threshold, and

[TABLE]

above it. We can see that $\mathfrak{D}(1)/2<1$ is always satisfied far below the threshold ($p\ll 1$), when $j>3$ satisfied at the threshold ($p\rightarrow 1$), and when $j>2$ satisfied far above the threshold ($p\gg 1$). When $2<j<3$, Duan-Giedke-Cirac-Zoller sufficient entanglement criterion is not satisfied at $p\sim 1$, even though it is satisfied well below and well above the threshold. Figure 5(c) shows the Duan-Giedke-Cirac-Zoller sufficient entanglement criterion for a nearest-neighbor pair in $N=32$ DOPOs with $j=7/3$. The analytical results [Eqs.(54) and (55)] agree very well with the numerical results. It is seen that the entanglement criterion ($\frac{\mathfrak{D}(r)}{2}<1$) is not satisfied at the threshold but is satisfied below and above it.

IV.4 Correlation as a function of distance $r$

Here, we consider the covariance matrix [Eq.(14)] between $\hat{a}_{1}$ and $\hat{a}_{1+r}$, calculated with $\overrightarrow{R}(r)=\sqrt{2}[\hat{X}_{1},\hat{P}_{1},\hat{X}_{1+r},\hat{P}_{1+r}]$. In the large $N$ limit, the small-$k$ part of the loss spectrum is important in calculating the spatial correlation function of the canonical components, $c_{1}=2\langle\Delta\hat{X}_{1}\Delta\hat{X}_{1+r}\rangle$ and $c_{2}=2\langle\Delta\hat{P}_{1}\Delta\hat{P}_{1+r}\rangle$. These are calculated from

[TABLE]

In the limit $N\rightarrow\infty$, these expressions are reduced to the integration forms:

[TABLE]

after the approximation $\cos\theta\sim 1-\frac{\theta^{2}}{2}$. Here, $\psi(p)=pH(1-p)+H(p-1)$ and $H(p)$ is the Heaviside’s step function. $m_{\pm}^{2}=\frac{2(1\pm p)}{j}$ below the threshold and $m_{\pm}^{2}=\frac{2(2p-1\pm 1)}{j}$ above the threshold. Assuming small $m_{\pm}$, the covariance matrix can be calculated with complex integral as

[TABLE]

We see that the correlation length $r_{c}$ of the $\hat{X}$ component satisfies $r_{c}\propto|1-p|^{-\frac{1}{2}}$. The correlation of $\hat{X}$ becomes long range at the threshold. This result resembles the correlation length $r_{c}\propto|1-\frac{T}{T_{c}}|^{-\frac{1}{2}}$ in Landau’s equilibrium phase transition theory Landau80 , where $T/T_{c}$ is a system temperature normalized by the critical temperature. The correspondence between an equilibrium phase transition and a non-equilibrium phase transition has been noted in Ref.DeGiorgio70 . For $\hat{P}$, on the other hand, the maximum correlation length is obtained at $p\rightarrow 0$: $r_{c}=\sqrt{j/2}$. Even for $p\rightarrow 0$, however, this correlation is short ranged. Duan-Giedke-Cirac-Zoller sufficient condition for entanglement for a pair with distance $r$ ($\mathfrak{D}(r)$) can be calculated from correlations shown above as

[TABLE]

Particularly, at the threshold ($p=1$), $\mathfrak{D}(r)$ is proportional to $r$ due to small $m_{-}$. With such dependence on $r$, the entanglement disappears for DOPO pair with large distance at the threshold.

Considering both nearest neighbor value and the long-range characteristics, we obtain the equation,

[TABLE]

[TABLE]

The theoretical Duan-Giedke-Cirac-Zoller sufficient condition for entanglement between DOPOs with site distance $r$ is shown in Fig.5(d), where $j=7/3$. Lines represent the values with complex integral and exact nearest-neighbor values[Eq.(63)(64)], and filled circles are calculated from positive-$P$ SDEs with $N=32$. Duan-Giedke-Cirac-Zoller entanglement criterion is satisfied only for $r=1$ even with $p\neq 1$. In Fig.5(e), necessary and sufficient condition for entanglement is presented. Entanglement criterion is satisfied even for nearest-neighbor pair for $p=1$. However, for pairs with large $r$, entanglement criterion is not satisfied. In Fig.5(d)(e), numerical results obtained by positive-$P$ calculation are shown with filled circles, which match the analytical results. At the threshold $p\rightarrow 1$, difference between the theory and numerical results emerges for large $r$ presumably because analytical results do not assume periodicity of one-dimensional lattice.

We next consider the quantum discord between two DOPOs with site distance $r$. The theoretical quantum discord between DOPOs with site distance $r$ is shown with lines in Fig.5(f), where $j=7/3$. Lines represent the values calculated with Eqs.(60)(61). Filled circles represent the numerical results obtained by positive-$P$ simulation. Even for large distance $r$, non-zero quantum discord can be formed at the threshold. This is because $c_{1}$ has a long-range correlation at the threshold. Even when $c_{2}$ is negligibly small for large $r$, the non-zero quantum discord can be obtained via finite $c_{1}$ value.

V Summary

We have presented analytical results on the degrees of squeezing/anti-squeezing, entanglement and quantum discord for a network of DOPOs both below and above threshold. We confirmed the validity of the analytical results through numerical simulations based on the positive-$P$, truncated Wigner, and truncated Husimi SDEs. The comparison establishes well the validity of the phase-space methods based on these SDEs. The $c$-number Langevin equations are essentially identical to the truncated Wigner SDE and valid for DOPOs below and above threshold. In a high-loss CIM based on optical delay-line coupling, entanglement between DOPO pulses disappears before reaching the thresholdTakata15 . It does, however, have quantum discord, indicating the existence of quantum correlations among the quantum fluctuations of DOPO pulses below and above threshold. Such quantum correlations disappear if we use a mean-field coupling that does not account for the fluctuations in the coupling fields. In the one-dimensional lattice of DOPOs, the entanglement criterion is satisfied only for short ranged pair of DOPOs. The quantum discord of far distant pair, on the other hand, has finite values particularly around the threshold, due to the long range correlation of canonical coordinates.

The analytical methods and insights obtained here for a 1D lattice of DOPOs with nearest-neighbor coupling seem to be applicable to other systems: DOPOs with measurement feedback couplingShoji17 , more complicated DOPO lattices with frustrationHamerly16 ; Leleu17 or topological nontriviality of the loss spectrum, and optical networks with other nonlinear components. For example, in the coherent XY machineTakeda17 with dissipative coupling, the entanglement has not been considered. Some lasers, for example Raman lasers, are known to have photon number squeezed state (sub-Poissonian state) only far above the thresholdYamamoto86 ; Ritsch92 . As the dissipatively coupled squeezed states can satisfy entanglement criterionTakata15 ; Maruo16 , the coherent XY machines could also satisfy entanglement criterion above the threshold, due to the dissipative coupling of the photon number squeezed states.

The linearization treatment developed in this paper is effective only for weak noise limit. They will not be applied with sufficient accuracy for large quantum noise case, where non-gaussianity Yamamura17 ; Drummond89 becomes important.

Acknowledgements.

The authors thank Z.Wang, T.Leleu, K.Takata, P.D.Drummond, Y.Huang, H.Mabuchi and E.Ng for their critical discussions. This work was supported partially by the ImPACT program of the cabinet office of the government of Japan.

Appendix A Phase-space method

We expand the density operator in the quantum master equation (5) by using the positive-$P$, Wigner, and Husimi quasi-distribution functions. First, we discuss Glauber’s diagonal $P$ distribution function, which is represented as $\hat{\rho}=\int P(\alpha)|\alpha\rangle\langle\alpha|d^{2}\alpha$Glauber63 . For a DOPO, the Fokker-Planck equation of the $P$ function is derived as :

[TABLE]

The diagonal $P$ function is not positive, however, and it does not have a corresponding $c$-number SDE because of the difficulty in adding noise term to $d\langle\alpha^{2}\rangle/dt$ without adding it to $d\langle|\alpha^{2}|\rangle/dt$. The positive-$P$ representationDrummond80 is a modification to Glauber’s $P$ representation and is known to have a similar Fokker-Planck equation. We can write the distribution function in positive-$P$ representation as Drummond80

[TABLE]

Here, $\alpha$ and $\alpha^{\dagger}$ are two independent complex numbers satisfying $\langle\alpha\rangle^{*}=\langle\alpha^{\dagger}\rangle$. Using Eq.(66) for Eq. (5), we obtain the quantum-mechanical Fokker-Planck equation for $P(\alpha,\alpha^{\dagger})$. This Fokker-Planck equation for the positive-$P$ representation is the same as that for Glauber’s diagonal $P$ distribution, except for the replacement $\alpha^{*}\rightarrow\alpha^{\dagger}$. Applying the Ito rule to the resulting Fokker-Planck equation, we can obtain the following $c$-number SDE for the signal modeShoji17 :

[TABLE]

[TABLE]

Here, $\xi_{R}$ and $\xi_{R}^{\dagger}$ are independent real-number random noise sources satisfying $\langle\xi_{R}(t)\xi_{R}(t^{\prime})\rangle=\delta(t-t^{\prime})$ and $\langle\xi_{R}^{\dagger}(t)\xi^{\dagger}_{R}(t^{\prime})\rangle=\delta(t-t^{\prime})$. This SDE under positive-$P$ representation is defined with no truncation and is equivalent to the quantum master equation (5). In numerical simulation, however, these equations have some difficulty for large $B/\gamma_{s}$. It is known that when $B$ is large, positive-$P$ calculation becomes unstable Gilchrist97 ; Deuar02 . This is because a negative photon number is allowed in positive-$P$ calculation, which causes the photon number to diverge to negative infinity because of two-photon absorption. In this paper, we thus consider the case of small $B$ and a sufficiently slow adiabatic excitation schedule. In our simulations, the divergence problem was not observed.

As noted above, the SDEs with positive-$P$ representation are equivalent to the Liouville equation (5). On the other hand, the Wigner and the Husimi functions for Eq.(5) have third and higher order derivatives in their Fokker-Planck equations. As these higher order derivatives are proportional to $B/\gamma_{s}$, we can neglect higher-order derivatives when $B/\gamma_{s}$ is small, and we obtain truncated SDEs. In return for the possible incorrectness due to truncation, the truncated Wigner and Husimi SDEs represent the signal boson operator $\hat{a}$ by only one $c$-number amplitude. The Wigner distribution function is written with a density matrix as the followingCahill69 ; Mishkin69 :

[TABLE]

The Fokker-Planck equation can be obtained from the time derivative of this equation. The truncated Wigner SDE is then obtained from the Fokker-Planck equation after neglecting third-order derivatives and terms resulting from the Weyl ordering in the drift and diffusion terms with the assumption of $B\ll\gamma_{s}$ :

[TABLE]

We can see that this Wigner SDE resembles the $c$-number counterpart of the Heisenberg-Langevin equation after eliminating the pump mode Wang13 , but Eq.(70) is the $c$-number SDE in the Schrödinger picture rather than the $q$-number equation in the Heisenberg picture.

The Husimi distribution function is then written as $Q(\alpha)=\frac{1}{\pi}\langle\alpha|\hat{\rho}|\alpha\rangle$. The Fokker-Planck equation of the Husimi distribution function for degenerate two-photon absorption is presented in Drummond81b , and it has third- and fourth- order derivatives. Here, we consider the truncated Husimi Fokker-Planck equation of Eq.(5) only up to second-order derivatives. As this approximation assumes small $B/\gamma_{s}$, we also neglect small $B/\gamma_{s}$ dependent terms resulting from the anti-normal ordering of drift and diffusion terms. The truncated Husimi SDE is obtained through the Fokker-Planck equation as the following:

[TABLE]

Here, $\xi_{C}$ is a complex-number noise source with $\langle\xi_{C}^{*}(t)\xi_{C}(t^{\prime})\rangle=2\delta(t-t^{\prime})$, and $\xi_{Ri}$ is a real-number noise source with $\langle\xi_{Ri}(t)\xi_{Rj}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$. We can use Eqs.(67) and (68) for the first principle numerical simulation of a CIM, and Eqs.(70) and (71) for the approximated simulation of a CIM with $B\ll\gamma_{s}$.

Appendix B Analytical method for single DOPO

Here, we derive the analytical solutions for the quantum noise of a CIM. We use the positive-$P$ representation, while noting that the truncated Wigner and truncated Husimi representations can produce identical results in the weak noise limit ($B/\gamma_{s}\ll 1$). Below the threshold, we neglect the two-photon absorption terms with $B$ in Eqs.(67) and (68). We can then obtain the following positive-$P$ SDEs:

[TABLE]

[TABLE]

The term with $B$ in the square root of Eqs.(67) and (68) can contribute even below the threshold when $B/\gamma_{s}$ is large, as $d\langle\alpha^{2}\rangle/dt=-(2\gamma_{s}+B)\langle\alpha^{2}\rangle+\cdots$. When $B$ is large, the saturation term $-B\alpha^{\dagger}\alpha^{2}$ in Eq.(67) is also not negligible even below the threshold, because the rate of two-photon absorption is enhanced by the intensity correlation function $g^{(2)}(0)=\frac{\langle\hat{a}^{\dagger 2}\hat{a}^{2}\rangle}{\langle\hat{a}^{\dagger}\hat{a}\rangle^{2}}$ Shen67 , and a squeezed vacuum state in a DOPO has large $g^{(2)}(0)$ proportional to $\frac{1}{2\langle\hat{a}^{\dagger}\hat{a}\rangle}$. We can neglect the saturation term, however, if $B\ll\gamma_{s}$. Below the threshold, the mean amplitude of the system is $\langle\alpha\rangle=0$. Using the normalized excitation $p=S/\gamma_{s}$, which is equal to $\varepsilon/\varepsilon_{thr}$, we can obtain the following steady-state equation:

[TABLE]

We assume that $\langle\Delta\alpha^{2}\rangle$ is real, because from Eqs.(72) and (73), $\alpha$ and $\alpha^{\dagger}$ are real: the time development starts from the vacuum state $\alpha=\alpha^{\dagger}=0$, and both the coefficients and random noises are real. We define the canonical coordinate and momentum as $\hat{X}=\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}}$ and $\hat{P}=\frac{\hat{a}-\hat{a}^{\dagger}}{\sqrt{2}i}$, respectively. In the positive-$P$ representation, the fluctuations are calculated as $\langle\Delta\hat{X}^{2}\rangle=\langle\Delta\alpha^{\dagger}\Delta\alpha\rangle+\langle\Delta\alpha^{2}\rangle+\frac{1}{2}$ and $\langle\Delta\hat{P}^{2}\rangle=\langle\Delta\alpha^{\dagger}\Delta\alpha\rangle-\langle\Delta\alpha^{2}\rangle+\frac{1}{2}$. In a DOPO below the threshold under the steady-state condition, these canonical components have the following fluctuations: $\langle\Delta\hat{X}^{2}\rangle=\frac{1}{2}+\frac{p}{2(1-p)}$ and $\langle\Delta\hat{P}^{2}\rangle=\frac{1}{2}-\frac{p}{2(1+p)}$. We can see that $\langle\Delta\hat{X}^{2}\rangle$ has a singularity at the threshold. On the other hand, $\langle\Delta\hat{P}^{2}\rangle$ has a minimized value $\langle\Delta\hat{P}^{2}\rangle=\frac{1}{4}$, which is one half the vacuum fluctuation, at the thresholdMilburn81 . Similar results can be obtained from other representations; for example, in the truncated Wigner representation, we set $B\rightarrow 0$ in Eq.(70). From $\begin{bmatrix}1&-p\\ -p&1\end{bmatrix}\begin{bmatrix}\langle\Delta\alpha^{2}\rangle\\ \langle|\Delta\alpha|^{2}\rangle\end{bmatrix}=\frac{1}{2}\begin{bmatrix}0\\ 1\end{bmatrix}$, we obtain the same mean-squared fluctuations via $\langle\Delta\hat{X}^{2}\rangle=\langle|\Delta\alpha|^{2}\rangle+\langle\Delta\alpha^{2}\rangle$ and $\langle\Delta\hat{P}^{2}\rangle=\langle|\Delta\alpha|^{2}\rangle-\langle\Delta\alpha^{2}\rangle$.

Next, we consider the fluctuations of a DOPO above the threshold. We separate the signal amplitudes into the mean amplitudes $\langle\alpha\rangle=\langle\alpha^{\dagger}\rangle=\sqrt{\frac{\gamma_{s}}{B}(p-1)}$ and small fluctuations ( $\Delta\alpha=\alpha-\langle\alpha\rangle$ and $\Delta\alpha^{\dagger}=\alpha^{\dagger}-\langle\alpha^{\dagger}\rangle$ ). Such linearization treatment was previously applied in quantum nonlinear optics for degenerate two-photon absorption Chaturvedi77 and an above-threshold DOPOMcNeil78 ; Drummond79 . We obtain the positive-$P$ SDE for the fluctuation part as follows.

[TABLE]

[TABLE]

Using the normalized excitation $p=\varepsilon/\varepsilon_{thr}$, we then obtain the following equation under the steady-state condition:

[TABLE]

The squared fluctuations of the canonical coordinate and momentum are obtained as $\langle\Delta\hat{X}^{2}\rangle=\frac{1}{2}+\frac{1}{4(p-1)}$ and $\langle\Delta\hat{P}^{2}\rangle=\frac{1}{2}-\frac{1}{4p}$, respectively. Above the threshold, the canonical coordinate and momentum monotonically reach the vacuum fluctuation. Although squeezing exists above the threshold, the squeezing becomes weaker for large $p$ Milburn81 . We can see that Eqs.(74) and (77) are reduced to an identical equation at the threshold (in the limit $p\rightarrow 1$). For the canonical momentum $\hat{P}$, the same mean-squared fluctuations are obtained at $p=1$ from both analytical theories below and above the threshold. For the canonical coordinate $\hat{X}$, the mean-squared fluctuation above the threshold is obtained from that below the threshold with substitution $p\rightarrow 1$ in the numerator and $1-p\rightarrow 2(p-1)$ in the denominator. For the canonical momentum $\hat{P}$, the mean-squared fluctuation above the threshold is obtained from that below the threshold with substitution $p\rightarrow 1$ in the numerator and $1+p\rightarrow 2p$ in the denominator.

Appendix C Traveling-pulse model of 1D lattice

We also consider the traveling-pulse model Marandi14 of the nearest-neighbor-coupled 1D lattice of DOPOs, shown in Fig.6. This model uses traveling optical pulses in a single ring cavity with a nonlinear crystal, in which the pulses are coupled with an optical delay-line. We assume that the two beam splitters (BSs) have reflectance $R_{B}$, the two half-beam splitters (HBSs) have $50\%$ reflectance, and the mirrors (Ms) have $100\%$ reflectance. The setting shown in Fig.6 has HBSs in the delay-line and is, as we show below, equivalent to the Liouvillian coupling of Eq.(38). We note that some models Hamerly16 ; Inagaki16 discussed in earlier papers are not equivalent to Eq.(38) because they have no HBSs in the delay-line. The traveling-pulse model is easily implemented with the Heisenberg Langevin model or SDE for the Wigner representation.

In the Wigner representation, the $r$-th signal pulse $\alpha_{r}$ is mixed with the vacuum noise $f_{1r}$ at a beam splitter (BS) and transformed into the transmitted $\alpha_{T,r}$ and reflected $\alpha_{R,r}$ modes: $\alpha_{T,r}=\sqrt{1-R_{B}}\alpha_{r}+\sqrt{R_{B}}f_{1r}$ and $\alpha_{R,r}=\sqrt{R_{B}}\alpha_{r}-\sqrt{1-R_{B}}f_{1r}$. The reflected light is then separated by an HBS into $\alpha_{D,r}=\frac{1}{\sqrt{2}}f_{2r}+\frac{1}{\sqrt{2}}\alpha_{R,r}$ and $\alpha_{A,r}=\frac{1}{\sqrt{2}}f_{2r}-\frac{1}{\sqrt{2}}\alpha_{R,r}$, where $f_{2r}$ is the independent vacuum noise. Here, $\alpha_{D}$ has a $1$-bit delay and $\alpha_{A}$ has a $-1$-bit delay. The optical delayline injects $\alpha_{A,r+1}$ and $\alpha_{D,r-1}$ into the $r$-th pulse $\alpha_{T,r}$. $\alpha_{F,r}=\frac{1}{\sqrt{2}}\alpha_{D,r-1}-\frac{1}{\sqrt{2}}\alpha_{A,r+1}$ is the mixed injection field from the delayline. The combined field at the injection beam splitter is represented by $\alpha^{\prime}_{r}=\sqrt{1-R_{B}}\alpha_{T,r}+\sqrt{R_{B}}\alpha_{F,r}$. Here, $f_{1r}$ and $f_{2r}$ are vacuum fields in the Wigner representation, satisfying $\langle f_{ar}^{*}f_{br^{\prime}}\rangle=\frac{1}{2}\delta_{ab}\delta_{r,r^{\prime}}$. Next, we assume that the DOPO cavity round-trip time is $\Delta t$, and that the delay-line’s reflectance is effectively simulated by a distributed coupling constant $J$, via $R_{B}=J\Delta t$. The traveling pulse with optical delay-line coupling model of Fig.6 is then represented as

[TABLE]

Here, $\xi_{C1,r}$ and $\xi_{C2,r}$ come from the input noise of the BS and HBS respectively, and satisfy $\langle\xi_{Ca,r}(t)\xi_{Cb,r^{\prime}}^{*}(t^{\prime})\rangle=2\delta_{ab}\delta_{r,r^{\prime}}\delta(t-t^{\prime})$. We can see that this traveling-pulse model produces the same diffusion terms in the Fokker-Planck equation, as Eq.(40) does. We assume here that $R_{B}=J\Delta t\ll 1$.