A Quantum Model for Coherent Ising Machine: Stochastic Differential Equations with Replicator Dynamics

TL;DR

This paper develops a quantum theoretical framework for coherent Ising machines using stochastic differential equations and replicator dynamics, enabling simulation of quantum effects in Ising spin models.

Contribution

It introduces a novel quantum model combining stochastic differential equations with replicator dynamics for coherent Ising machines.

Findings

Simulated simple Ising spin models demonstrating quantum effects.

Elucidated unique features of quantum coherent Ising machine.

Provided a theoretical basis for future quantum optimization algorithms.

Abstract

The quantum theory of coherent Ising machines, based on degenerate optical parametric oscillators and measurement-feedback circuits, is developed using the positive $P({\alpha},{\beta})$ representation of the density operator and the master equation. The theory is composed of the c-number stochastic differential equations for describing open dissipative quantum dynamics and the replicator dynamics equations for handling measurement-induced collapse of the density operator. We apply the present theory to simulate two simple Ising spin models and elucidate the unique features of this computing machine.

Equations

$$\mathcal{H}=\sum _{i,j\in E}J(i,j)\sigma (i)\sigma (j),$$

$$\rho ={\rho }_1\otimes {\rho }_2,$$

$$\hat{H}$$

$${\hat{H}}_{free}$$

$${\hat{H}}_{int}$$

$${\hat{H}}_{pump}$$

$${\hat{H}}_{FB}$$

$${\hat{H}}_{loss}$$

$$\frac{d{\hat{\rho }}_{DOPO}}{dt}$$

$$\frac{d{\hat{\rho }}_{meas}}{dt}$$

$$d{W}_i(t)$$

$$⟨d{W}_i(t)d{W}_j(t^{\prime })⟩$$

$$X_{i}dt=⟨{\hat{a}}_{si}^{†}+{\hat{a}}_{si}⟩dt+\frac{d{W}_i}{\sqrt{\xi }}.$$

$${ϵ}_{si}=\zeta \sum _j{J}_{ij}{X}_j,$$

$$\frac{d\hat{\rho }}{dt}$$

$$\hat{\rho }=\int \int P(\alpha ,\beta )\hat{\Lambda }(\alpha ,\beta )d^2\alpha d^2\beta ,$$

$$\hat{\Lambda }(\alpha ,\beta )=\frac{|\alpha ⟩⟨\beta |}{⟨\alpha |\beta ⟩}$$

$$\int d^2\alpha d^2\beta P(\alpha ,\beta )=1.$$

$$⟨{\hat{a}}^{†n}{\hat{a}}^m⟩=\int d^2\alpha d^2\beta {\alpha }^m{\beta }^{n}P(\alpha ,\beta ).$$

$$Prob(X)$$

$$P(\alpha ,\beta )$$

$$\frac{\partial P(\alpha ,\beta )}{\partial t}$$

$$⟨{\alpha }_{si}⟩=\int d^2\alpha d^2\beta {\alpha }_{si}P(\alpha ,\beta ),$$

$$⟨{\beta }_{si}⟩=\int d^2\alpha d^2\beta {\beta }_{si}P(\alpha ,\beta ).$$

$$\displaystyle\left[\begin{array}[]{c}d\alpha_{si}\\ d\beta_{si}\end{array}\right]$$

$$+$$

$$\displaystyle\left[\begin{array}[]{c}d\alpha_{pi}\\ d\beta_{pi}\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}d\alpha_{si}\\ d\beta_{si}\end{array}\right]$$

$$+$$

$$\displaystyle\left[\begin{array}[]{c}d\eta_{i}\\ d\mu_{i}\end{array}\right]=\left[\begin{array}[]{c}-(1+\xi^{\prime})\eta_{i}+\mu_{i}(p_{i}-\eta_{i}^{2})+f_{i}\\ -(1+\xi^{\prime})\mu_{i}+\eta_{i}(p_{i}-\mu_{i}^{2})+f_{i}\end{array}\right]d\tau$$

$$\displaystyle+\left[\begin{array}[]{c}g\sqrt{p_{i}-\eta_{i}^{2}}d\omega_{\eta_{i}}\\ g\sqrt{p_{i}-\mu_{i}^{2}}d\omega_{\mu_{i}}\end{array}\right].$$

$$\frac{\partial P(\alpha ,\beta )}{\partial t}=\lambda (\alpha ,\beta )P(\alpha ,\beta ).$$

$$\displaystyle\left\{\begin{array}[]{ll}\mathrm{copied\ with\ probability\ }\lambda(\alpha,\beta),&(\lambda(\alpha,\beta)>0),\\ \mathrm{vanished\ with\ probability\ }-\lambda(\alpha,\beta),&(\lambda(\alpha,\beta)<0),\end{array}\right.$$

$$\frac{\partial P(\alpha ,\beta )}{\partial t}$$

$$d⟨{\alpha }_{si}⟩$$

$$d⟨{\beta }_{si}⟩$$

Full text

I Introduction

Various combinatorial optimization problems belong to the NP-hard or NP-complete class and are difficult to solve in a polynomial time with a deterministic Turing machineNPhard . To find exact or approximate but satisfying solutions for such problems, many heuristic algorithms, such as classical neural networksHopfield HopfieldTank , simulated annealingSimulatedAnnealing , and quantum annealingFarhi QuantumAnnealing , are proposed.

The Ising problem is the simplest model for finding the minimum energy ground state of spin glassesIsingOriginal . The Ising Hamiltonian is given by

[TABLE]

where $G=(V,E)$ is a given graph, $J:E\to\mathbb{R}$ is the weight of an edge, and $\sigma:V\to\{-1,1\}$ is the vertex value, called an “Ising spin.” It is known that a three-dimensional Ising model and two-dimensional Ising model with a Zeeman field are NP-hard problems.Ising

There have been several attempts to solve Ising problems with actual physical devices, called Ising machines, rather than algorithms inspired by physical phenomena. A coherent Ising machine (CIM) is one such physical deviceInjectionCIM1 InjectionCIM2 OPOCIM1 . The first generation of CIM uses injection-locked lasers to represent Ising spins,InjectionCIM1 InjectionCIM2 while the second generation employs degenerate optical parametric oscillators (DOPOs)OPOCIM1 . The coupling $J(i,j)$ between Ising spins is implemented by the optical delay line coupling between the oscillatorsNaturePhotonics TakataSciRep InagakiNature .

Recently, a measurement-feedback circuit has been used to implement $J(i,j)$ efficiently. The exprimental system with N DOPO pulses shown in Fig.1(a) Science1 Science2 is modeled and sinplified as shown in Fig 1(b). In this study, we develop the quantum theory of such measurement-feedback based CIMs based on the positive $P(\alpha,\beta)$ representation of the density operator and the c-number stochastic differential equations (CSDE) with replicator dynamics. The accompanying paper addresses the same problem with a different theoretical model based on a completely positive trace preserving mapYamamura .

II Master equation for CIM

In this section, we will derive the master equation of a measurement-feedback-based CIM. For the sake of simplicity, we describe the derivation of the master equation for only a two-spin system, but the theory can be easily extended to many-spin systems.

II.1 A simple model of CIM

Figure 1(b) shows a simplified model of the CIM. There are two DOPO cavities with the identical signal frequency $\omega_{s}$ and pump frequency $\omega_{p}=2\omega_{s}$. The photon annihilation operators of the signal and pump fields are denoted by $\hat{a}_{s1},\hat{a}_{p1},\hat{a}_{s2}$, and $\hat{a}_{p2}$. There are also two external fields injected into the cavities. One is an excitation pump field at $\omega_{p}$. Another is a feedback signal field at $\omega_{s}$, which is prepared by the measurement-feedback circuit.

The intra-cavity pump field and signal field have the loss rates denoted by $\gamma_{p}$ and $\gamma_{s}$, respectively. To measure the in-phase component $\hat{X}_{i}=\hat{a}_{si}+\hat{a}_{si}^{\dagger}(i=1,2)$ of the signal field, a part of the intra-cavity field is picked off and measured by homodyne detectors. The feedback signal is prepared based on the measurement results.

Since the measurement-feedback process is local operation and classical communication (LOCC), the density matrix of the total system stays in a product state during the whole computation process and is given as follows:

[TABLE]

where $\rho_{i}$ is a density matrix of the DOPO $i$.

Our theretical model is a continuous time evolution model, in which all the quantum operations proceed simultaneously, while a real measurement-feedback-based CIM Science1 Science2 is based on discrete quantum operations, as shown in Fig. 1(a). Extension of the present work to describe a discrete model is straightforward and will be discussed elsewhere.

In a real CIM, each DOPO takes the form of a pulse circulating in a fiber ring resonator and is operated sequentially in time. Let $\tau$ be time interval of DOPO pulses and consider a set of DOPO pulse amplitudes $\{\hat{a}_{s1}(t),\hat{a}_{s2}(t+\tau),...,\hat{a}_{sN}(t+(N-1)\tau)\}$, which appear at a same spatial point (injection coupler) in a real experimental system. The internal target pulse and externally prepared feedback pulse collide at the injection coupler with the exactly same delay time, so that a delay of feedback signal in the theoretical morel can be assumed to be zero.

II.2 Derivation of the master equation

To acquire the stochastic differential equations for the CIM, we treat a DOPO and a measurement-feedback circuit separately.

The Hamiltonian of DOPOs is given by

[TABLE]

where $\kappa$ is a parametric coupling constant between the signal field and the pump field in a nonlinear crystal, and $\Gamma_{pi}$ and $\Gamma_{si}$ are the external reservoir field operators, which account for the fluctuation forces injected from the external environment. By tracing out these external fields by the standard Born–Markov approximationGerdiner , we can obtain the following master equation of the DOPOs:

[TABLE]

By taking a rotating reference frame properly, we can eliminate the two terms in the first line of Eq.(9).

To describe the nonunitary reduction of a wave function by the homodyne measurement of $\hat{x}$, Wiseman and Milburn proposed the following equationWisemanFeedback :

[TABLE]

where $\left<\hat{X}_{i}\right>=\left<\hat{a}_{si}^{\dagger}+\hat{a}_{si}\right>$ is the expectation value of the in-phase amplitude of the signal field and $dW$ is the Wiener increment, which satisfies

[TABLE]

In this model, the actually measured value $X_{i}$ is given by

[TABLE]

A feedback signal $\epsilon_{si}$ is now prepared according to the formula,

[TABLE]

where $\zeta$ is the strength of the feedback coupling between DOPOs.

Finally, we obtain the overall master equation of the measurement-feedback-based CIM by combining Eq.(9) and (10) as follows:

[TABLE]

III Stochastic differential equations

III.1 Positive $P(\alpha,\beta)$ representation

The $P(\alpha,\beta)$ representation of the density operator is defined byOPOhamiltonian

[TABLE]

where $\alpha,\beta\in\mathbb{C}$ and

[TABLE]

is the off-diagonal projector in terms of coherent states. Here, $|{\bf\alpha}\left.\right>$ and $|{\bf\beta}\left.\right>$ are the tensor product coherent states: $|{\bf\alpha}\left.\right>=|\alpha_{s1}\left.\right>\otimes|\alpha_{p1}\left.\right>\otimes|\alpha_{s2}\left.\right>\otimes|\alpha_{p2}\left.\right>$, and $|{\bf\beta}\left.\right>=|\beta_{s1}\left.\right>\otimes|\beta_{p1}\left.\right>\otimes|\beta_{s2}\left.\right>\otimes|\beta_{p2}\left.\right>$.

An important property of the positive $P(\alpha,\beta)$ representation is that we can always define $P({\bf\alpha},{\bf\beta})$ as a positive and real function for arbitrary quantum states of fields, which satisfies

[TABLE]

Therefore, we can regard $P(\alpha,\beta)$ as a probability distribution function for finding the projector $|\alpha\left.\right>\left<\right.\beta|$ in the density matrix.

Since a coherent state is an eigenstate of an annihilation operator, the moment of the density matrix is easily evaluated by

[TABLE]

Using this relation, we can obtain any statistics of an observable, which is composed by the creation and annihilation operators such as $\left<\hat{X}\right>=\left<\hat{a}+\hat{a}^{\dagger}\right>$.

The probability density function of finding the in-phase component $X$ can be also computed from the positive P representation as follows:

[TABLE]

Because the density matrix of the measurement-feedback-based CIM is given as the product state of each DOPO density matrix, we can express $P(\alpha,\beta)$ as

[TABLE]

Therefore, we can describe the total system with the partial differential equation (PDE) for each DOPO.

Using the properties of the coherent states, we obtain the PDE of $P(\alpha_{si},\beta_{si},\alpha_{pi},\beta_{pi})$ from the master equation as follows:

[TABLE]

where

[TABLE]

III.2 Stochastic differential equations and replicator dynamics

Except for the first line in Eq. (22), the PDE has an identical form as the Fokker–Planck equation. It is well established that the Fokker–Planck equation can be transformed to the following stochastic differential equations (SDE)Gerdiner :

[TABLE]

When $\gamma_{p}>>\gamma_{s}$, the pump field decays more rapidly than the signal field, so that the pump field follows the dynamics of the signal field (the slaving principle). We can eliminate the pump field by assuming $d\alpha_{pi}=d\beta_{pi}=0$,

[TABLE]

By introducing $\eta_{i}=g\alpha_{si},\mu_{i}=g\beta_{si},d\tau=\gamma_{s}dt,d\omega_{\eta_{i}}=\sqrt{\gamma}_{s}dW_{\alpha_{i}},d\omega_{\mu_{i}}=\sqrt{\gamma}_{s}dW_{\beta_{i}},p_{i}=\kappa\epsilon_{pi}/\gamma_{s}\gamma_{p},f_{i}=g\epsilon_{si}/\gamma_{s}$, and $g=\kappa/\sqrt{2\gamma_{p}\gamma_{s}},\xi^{\prime}=\xi/\gamma_{s}$, we obtain the normalized SDE

[TABLE]

Note that above the oscillation threshold, fixed points of in-phase component $\left<\hat{X}_{i}\right>=\left<\alpha_{i}+\beta_{i}\right>$ are $\pm\sqrt{p-1}/g$, and that of mean photon number $\left<\hat{n}\right>$ is $(p-1)/g^{2}$. Therefore, to tune the present theoritical model to a real experiment, we can set $g$ from the measured photon number at a pump rate $p=2$.

On the other hand, the first line of Eq. (22), which describes the reduction of the signal density operator induced by the measurement, cannot be simulated by the standard method using the SDE. In previous worksWisemanFeedback , by assuming that the measurement result is incidentally identical to the expectation value, they ignored this term. In this study, we need to know the measurement effect on the evolution of the DOPO state, so that we keep the random measurement results on $\eta_{i}(t)$ and $\mu_{i}(t)$ by taking pseudo-random numbers $dW$.

Because the first line of Eq. (22) is a replicator equation, we extend the branching Brownian motion model BranchingBrownianMotion , which is called replicator dynamics in our case. In replicator dynamics, the change of $P(\alpha,\beta)$ is governed by

[TABLE]

Here a Brownian particle at $(\alpha,\beta)$ is

[TABLE]

where $\lambda(\alpha,\beta)=X_{i}-\left<\hat{X}_{i}\right>$.

Because the expectation values $\left<\alpha_{si}+\beta_{si}\right>$ are needed to compute $\lambda(\alpha,\beta)$, we run many Brownian particles according to the identical SDEs and the same measurement results $X_{i}$ at the same time.

III.3 Gaussian approximation

In this section, we derive an approximation method to describe the measurement-feedback-based DOPO system.

We start from the following PDE of the signal fields after the adiabatic elimination of the pump field:

[TABLE]

By partial integration of Eq. (31), the equations of motion for the expectation values are obtained as follows:

[TABLE]

Similarly, we can derive the equations of motion for the higher order statistics such as $\left<\alpha_{si}^{2}\right>,\left<\beta_{si}^{2}\right>$, and $\left<\alpha_{si}\beta_{si}\right>$. Even though $\epsilon_{si}$ contains the statistics of other DOPOs such as $\left<\alpha_{sj}\right>(i\neq j)$, $\left<\alpha_{si}\alpha_{sj}\right>$ is expressed as $\left<\alpha_{si}\right>\left<\alpha_{sj}\right>$ since the total density matrix is separable. By changing the basis from $\left<\alpha_{si}\right>and\left<\beta_{si}\right>$ to $\left<X_{i}\right>=\left<\alpha_{si}\right>+\left<\beta_{si}\right>and\left<iP_{i}\right>=\left<\alpha_{si}\right>-\left<\beta_{si}\right>$, the equations of motion of $\left<X_{i}\right>$ and $\left<P_{i}\right>$ are acquired.

Though the dynamical equations are acquired, we cannot simulate Eq. (32) and (33) immediately because of higher order terms like $\left<\alpha_{si}^{2}\beta_{si}\right>$. To avoid this difficulty, we consider the approximate wave function given by a displaced squeezed vacuum state

[TABLE]

where $\hat{D}(\mu_{i})=\exp(\mu_{i}\hat{a}^{\dagger}-\mu_{i}^{*}\hat{a})$ is the displacement operator, and $\hat{S}(\sigma_{i}^{2})=\exp(1/2(\sigma_{i}^{2*}\hat{a}^{2}-\sigma_{i}^{2}\hat{a}^{{\dagger}2}))$ is the squeezing operator. For simplicity, both $\mu$ and $\sigma_{i}^{2}$ are real. This approximation means that the DOPO state is always described by a pure squeezed state and higher order statistics has no effect on the dynamics of the system. By this approximation, we can finally get the dynamical equations of motion of the DOPOs as follows:

[TABLE]

where

[TABLE]

The first term in Eq. (35) describes the shift of the center position of the wave function by the measurement. The last term in Eq.(36) describes the reduction of the variance by the measurement.

IV Numerical simulation results

In this section, we will show the numerical simulation results based on exact CSDE (28) and replicator dynamics of Eq. (29) and compare them with the simulation results by the Gaussian approximation based on Eqs. (35) and (36).

IV.1 N=2 DOPO model

First, we study the system of two DOPOs with antiferromagnetic coupling. Figure 2 shows the time evolution of the average in-phase amplitudes $\left<X_{i}\right>$ vs. normalized pump rate $\epsilon_{p}/\epsilon_{th}$, where $\gamma_{s}=1,\gamma_{p}=10,\kappa=0.1,\xi=0.1$, and $\zeta=0.3$. The corresponding saturation parameter is $g=0.02$. Figure 2(a)-(f) shows the result of exact replicator dynamics discussed in sec.II-B, and Figure 2(g)-(h) shows that of Gaussian approximation discussed in sec.II-C. The external pump rate $\epsilon_{p}$ is linearly increased from 0 to 1.5 times the threshold value $\epsilon_{th}$. The two DOPOs are coupled by the antiferromagnetic interaction ($J_{ij}=-1$). Figure 2(b) expands the average in-phase amplitudes $\left<X_{1}\right>$ and $\left<X_{2}\right>$ near the bifurcation point (decision making point). The two DOPOs point toward one ground state $\left.|\uparrow\downarrow\right>$ at one time but switch back to the other ground state $\left.|\downarrow\uparrow\right>$ at another time. This random search process continues until the final decision is made at $\epsilon_{p}/\epsilon_{th}\simeq 0.9$.

Figure 2(c) shows the average photon number $\left<n\right>$ versus the normalized pump rate $\epsilon_{p}/\epsilon_{th}$. The lower bound of the average photon number below the threshold in Fig. 2(f) is associated with the squeezed vacuum state, while the noisy spikes correspond to the finite in-phase amplitude induced by the measurement-feedback process. Note that the average photon number per DOPO is on the order of one at this decision making instance, as shown in Fig.2(c).

Figure 2(d) shows the measurement results actually reported by the homodyne detectors. We conclude that the negative correlation is formed between $\left<X_{1}\right>$ and $\left<X_{2}\right>$ at very early stages by the measurement-feedback process as shown in Fig.2(b), but the actual measurement results are too noisy to disclose those quantum search processes. A final solution, which the CIM will report should be determined at the effective threshold pump rate $\epsilon_{p}/\epsilon_{th}=1-\zeta=0.7$. We will discuss in the next section that the late bifurcation at $\epsilon_{p}/\epsilon_{th}=0.9$ in Fig.2(b) rather than $0.7$ stems from the quantum tunneling.

The variance $\left<\Delta X_{i}^{2}\right>$ and skewness$\left<\Delta X_{i}^{3}\right>$ in the anti-squeezed in-phase amplitudes are shown in Fig.2(e) and (f). The DOPO wavepackets near and above the threshold are clearly deviated from the Gaussian wavepackets, for which $\left<\Delta X_{i}^{3}\right>=0$ holds. This is because the bottom of the potential function, $V_{b}(X_{i})=-\frac{1}{2}(\epsilon_{p}/\epsilon_{th}-1)X_{i}^{2}+\frac{1}{4}X_{i}^{4}$, above the threshold does not have a symmetric barrier, i.e. steep barrier toward a large amplitude $|X_{i}|\neq 0$ and gradual barrier at zero amplitude $|X_{i}|=0$. However, the wavepackets below the threshold differ only slightly from the Gaussian wavepackets. Figures 3(a) and 3(b) show the probabilities $P(x>0)$ and $P(x<0)$ of the DOPO wave functions, respectively, and compare them with those of the Gaussian wavepackets. These results suggest that the tails of the DOPO wavepackets toward $|X_{i}|=0$ are broader than those of the Gaussian wavepackets even at a pump rate below the threshold.

Figures 2(g) and (h) show $\left<X_{i}\right>$ based on the Gaussian approximation described above, and Figure 4(i) shows the variances $\sigma_{i}^{2}$s. We can find that the final correlation between $\left<X_{1}\right>$ and $\left<X_{2}\right>$ is formed already at the effective threshold $\epsilon_{p}/\epsilon_{th}=0.7$. Compared with the numerical results of the exact replicator dynamics, the bifurcation of $\left<x_{1}\right>$ and$\left<x_{2}\right>$, namely a final decision making occurs at the threshold $\epsilon_{p}/\epsilon_{th}=1-\zeta=0.7$. Please note a remarkable difference between Fig.2(b) for exact dynamics and Fig.2(h) for Gaussian approximation.

Note that the deviation from Gaussian wavepackets and the occurrence of skewness are also observed in non degenerate optical parametric oscillator (NDOPO) systems and contribute the entanglement generation in NDOPO systemsReferees1 . Existence of nonlinear effect in optical systems yields non Gaussian distribution which leads to deference from classical systemsReferees2 Referees3 .

IV.2 N=16 DOPO system

To reveal the unique capability of the CIM as an optimizer, we simulated N=16 DOPOs coupled by the nearest neighbor antiferromagnetic interaction in a one-dimensional ring configuration:

[TABLE]

The two degenerate ground states of this model are $(-1,1,-1,...,-1,1)$ and $(1,-1,1,...,1,-1)$.

We compared three models. The first and second models are based on the exact replicator dynamics and the Gaussian approximation. The third model is based on the quantum theory of an optical delay line coupling CIM analyzed by TakataPRA Maruo .

Figure 4(a) shows the success rates of finding either one of the two degenerate ground states in 1000 trials, where $\gamma_{s}=1,\gamma_{p}=10,\kappa=0.1,\xi=0.1$ and the injection rate $\zeta$ changes from $0.01$ to $1$. $\epsilon_{p}$ is linearly increased from 0 to 1.2 times $\epsilon_{th}$.

When the mutual coupling parameter $\zeta$ is small, the potential landscape for the DOPO field is almost symmetric with respect to $X=0$ as shown in Fig. 6. In such a case, the measurement-induced wavepacket reduction and the feedback-induced wavepacket displacement play major roles in the solution search process. In this case, the tightly confined Gaussian wavepacket is more advantageous than the broadly spread exact wavepacket, because the latter introduces more random measurement results and takes a longer time to reach a final decision. However, when the mutual coupling parameter $\zeta$ is large, the potential landscape of the DOPO field is highly asymmetric. In such a case, the decision making process in the Gaussian approximation happens too quickly so that the system is easily trapped in a wrong state as shown in Fig.5(b). The non-Gaussian wavepacket induced quantum tunneling in the exact replicator dynamics plays an important role for the system to escape from the wrong solution. In this case, the broadly spread exact wavepacket is more advantageous than the tightly confined Gaussian wavepacket, as shown in Fig. 5(b). The numerical simulation results in Fig.4(a) confirm this trade-off relation and also suggest the importance of quantum tunnelingKimsler in the solution search process of the CIM near the threshold.

The numerical results in Fig.4(b) shows that the optical coupling CIM is more efficient than the measurement-feedback-based CIM when $\zeta$ is small. However, in the case of $\zeta>0.6$, the measurement-feedback-based CIM has a higher success rate than the optical delay line coupling CIM. When the connection strength $\zeta$ is close to one, it means that the extracted signal field is boosted by a high-gain phase sensitive amplifier before it is injected back to the DOPO cavity (Fig.1 of ref.Maruo ). This is necessary since the injection coupler has a very small coupling constant. During this external amplification process, the vacuum fluctuation added to the extracted signal field is also amplified and contributes to the degradation of the degree of negative correlation among neighboring DOPOs. Because of this reason, there is an optimum coupling strength to maximize the degree of correlation in the optical delay line coupling CIM (see Fig.3(b) of ref.TakataPRA ). The maximum success rate at $\zeta\simeq 0.5$ corresponds to this optimum coupling strength. In the case of the measurement-feedback-based CIM, the search mechanism is not the formation of correlation between DOPOs but the feedback signal-induced quantum tunneling so that a higher coupling strength $\zeta$ always improves the success rate.

Fig.4 (c) and (d) shows the same comparison as Fig.4 (a) and (b) where $\kappa=0.01$ and $g=0.002$, which means that non-linear coupling constant $\kappa$ and the saturation parameter g are ten times smaller. The tendency of the results are not so different from the case of $\kappa=0.1$ and $g=0.02$. The success rate of $\kappa=0.01$ and $g=0.002$ case is higher than that of $\kappa=0.1$ and $g=0.02$ case, Therefore, we can conclude that the saturation parameter $g$ plays a crucial role in the computation.

V Conclusion

We presented the two theoretical methods for numerically simulating the measurement-feedback-based CIM with exact replicator dynamics and Gaussian approximation.

The Gaussian approximation method is computationally inexpensive, that is, we can efficiently implement this model in modern digital computers, while the exact replicator dynamics is computationally expensive because we have to run many Brownian particles simultaneously to properly account for the non-unitary state reduction of non-Gaussian wavepackets. As we have shown, the exact replicator dynamics predicts a better performance than the Gaussian approximation when the mutual coupling strength is large. Therefore, the present result supports the current experimental effort to develop an actual physical device solving replicator dynamicsScience1 Science2 rather than using the CSDE under Gaussian approximation as a new algorithm.

Acknowledgement

This research was funded by ImPACT Program of Council for Science, Technology and Innovation (Cabinet Office, Government of Japan). The authors wish to thank P. Drummond, H. Mabuchi, R. Hamerly and A. Yamamura for their critical discussions.