## Abstract

We propose a mapping protocol to implement Ising models in injection-locked laser systems. The proposed scheme is based on optical coherent feedback and can be potentially applied for large-scale Ising problems.

©2011 Optical Society of America

## 1. Introduction

The emergence of optimization problems is ubiquitous in our modern life, for which we have to find the best solution by choosing the particular combination of *M* variables *σ*
_{1}, *σ*
_{2}, ... *σ _{M}* under some constraints. In most cases, such a mission can be formulated as a computational problem of minimizing a cost function

*E*(

*σ*

_{1},

*σ*

_{2}, ...

*σ*). The mission is to find the specific values for those

_{M}*M*variables for which the function

*E*(

*σ*

_{1},

*σ*

_{2}, ...

*σ*) has the minimum value. If such an optimization problem can be solved in polynomial time, i.e. the computational time varies polynomially as a problem size

_{M}*M*, by only nondeterministic machines, such a problem is said to belong to the class NP (NP for nondeterministic polynomial). Among NP problems, there are certain subsets of problems known as NP-complete problems. Any NP problems can be mapped to NP-complete problems using a polynomial step. The Graph Partition Problem (GPP) or MAX-CUT problem is a representative of the NP-complete problem. If one has an efficient machine to solve a certain NP-complete problem of size

*M*, then one can solve any NP problems with an extra cost in time that scales only polynomially with

*M*. However, problems in NP-complete class are considered to be computationally hard, since such a nondeterministic machine cannot be simulated by a deterministic Turing machine without an exponential growth of computational time [1].

The above mentioned GPP or MAX-CUT problem can be formulated as a ground state search problem of an Ising spin model [2]:

where*σ*denotes the Ising spin (projection along

_{iz}*z*-axis) at a site

*i*, whose values can be +1 (up) or −1 (down), and

*J*the magnetic interactions between them. An Ising model also provides a prototype framework for studying various magnetic orders in frustrated spin lattice and random spin glasses [3, 4]. Hence, an Ising machine that can find a ground state of Eq. (1) efficiently has been extensively searched in both classical and quantum domains [5].

_{ij}Quantum annealing (or quantum adiabatic computation) is proposed to solve Ising models by utilizing quantum uncertainty, more specifically quantum mechanical tunneling across a potential energy landscape (PEL) [6–11]. Experimental realization of quantum annealing employed either a sample of real magnetic crystal [12, 13] or molecular NMR technique [14]. However, in order to map a given mathematical problem onto such a quantum annealing Ising machine, non-local interaction *J _{ij}* must be implemented artificially irrespective of actual distance between two sites. This is extremely hard to achieve in real crystals or molecules.

An Ising machine that harnesses the effect of Bose-Einstein condensation (BEC) and measurement-feedback control has been recently proposed [15]. The proposed machine utilizes the bosonic final state stimulation property of BEC to attain the speedup in cooling down to the ground state by a factor of bosonic particle number per site [16]. However, a measurement-feedback circuit which is required to implement nonlocal interactions *J _{ij}* is an inherently incoherent device and cannot create the quantum coherence between different sites. The cost resulting from this fact is that the computational time scales exponentially as a problem size of

*M*[16]. After all, the Ising model, itself, is NP-complete [17].

In this paper, we propose a new Ising machine based on one master laser and *M* mutually injection-locked slave lasers. An Ising model is implemented by coherent feedback network using optical interference circuits instead of incoherent electrical measurement-feedback circuits. A spin degree of freedom *σ _{iz}* at each site is represented by right or left circular polarization states of each slave laser. The ground state of an implemented Ising model which is represented with polarization configuration of

*M*slave lasers emerges spontaneously through the natural mode competition induced by cross-gain saturation among all candidate polarization configurations. The cost function

*E*(

*σ*

_{1},

*σ*

_{2}, ...,

*σ*) mentioned above corresponds to the over-all photon decay rate in our system. The injection-locked laser system oscillates with a specific polarization configuration which minimizes the cost function (photon decay rate).

_{M}A semiconductor laser is particularly attractive for this application because of its rapid intra-band spin relaxation, small saturation power, short photon lifetime and compact size. A rapid spin relaxation process realizes equally populated conduction electron spins and valence hole spins, so that the particular polarization configuration with a minimum photon decay rate can oscillate alone by suppressing the oscillation of all the other 2* ^{M}* – 1 polarization configurations. A saturation photon number

*n*(or inverse fractional spontaneous emission coupling efficiency $\frac{1}{\beta}={n}_{s}$) of a semiconductor laser is many orders of magnitude smaller than that of the other types of laser systems. Therefore, the above cross-gain saturation can be switched on by an extremely small injection power. A photon lifetime of typical semiconductor lasers is ∼ 1

_{s}*psec*, which allows a very large injection-locking bandwidth even with a very small injection power. A stable and fast operation of this computing machine originates from this fact that an injection-locking bandwidth is extremely broad. Finally a large number of arrayed vertical cavity surface emitting semiconductor lasers (VCSEL) can be integrated into one chip. It is even possible to integrate hybridly those VCSELs with a master laser, other linear optical circuits and detectors, which makes a whole computing system very compact even for a large-scale Ising problem

*M*≫ 1.

## 2. Proposed system

Our goal is to map the following Ising Hamiltonian into an externally controllable physical system:

*σ*describes an Ising spin projection on z-axis.

_{iz}*J*is a mutual interaction term between spin

_{ij}*i*and spin

*j*, and

*λ*is a supplemental Zeeman (external magnetic field) term. The injection-locked laser system we propose in this paper can find the ground state of the Hamiltonian Eq. (2) by the spontaneous formation of a macroscopic population of photons in a specific polarization configuration through a laser phase transition. Every photon in the lasing mode is not localized in any specific slave laser but coherently spreads its wavefunction into all slave lasers as partial waves. The polarization state of each slave laser represents the ground state {

_{i}*σ*,

_{iz}*i*= 1 ∼

*M*}.

We use circular polarization degrees of freedom in VCSELs to represent the Ising spin information. An effective spin state *σ _{iz}* at site

*i*is experimentally determined by a majority vote using photodetection signals at the final step of this computational scheme, i.e.,

*n*and

_{Ri}*n*are the number of photons with right and left circular polarizations detected at site

_{Li}*i*(from

*i*-th slave laser.) We assume that all slave lasers are driven by the same pump power.

A proposed injection-locked laser system is shown in Fig. 1. All slave lasers, which play a role of respective Ising sites, are injection-locked by a single master laser with vertical linear polarization. VCSELs can be used as slave lasers. An external cavity controlled diode laser (ECDL) can be used as a master laser. The injection-locking bandwidth is given by [18, 19]

where $\frac{\omega}{Q}$ is the cavity photon decay rate of the slave laser and on the order of ∼10^{12}( $\frac{1}{s}$) for the case of a VCSEL,

*P*is the injection signal power from the master laser into the slave laser and

_{in}*P*is the self-oscillation power of the slave laser. We assume that the linewidth enhancement factor (or anomalous dispersion parameter) is zero for simplicity [20]. The locking bandwidth is on the order of 0.1 ∼ 1

_{out}*GHz*in the later discussions on numerical examples.

A master laser output is split into *M* paths and injected into *M* slave lasers. Linear optical circuits can implement the Zeeman term and Ising interaction term in Eq. (2). To implement the Zeeman term *λ _{i}*, the master laser output is injected into the

*i*-th slave laser with a weak horizontal linear polarization component as well as a strong vertical polarization component. The amplitude and phase of the horizontally polarized injection signal is controlled with half-wave plates (HWP) and quarter-wave plates (QWP) as shown in Fig. 1. The Ising interaction term between site

*i*and site

*j*is implemented with the mutual injection of the slave laser outputs via an attenuator, phase shifter and horizontal linear polarizer. In total, $\frac{1}{2}M(M-1)$ paths must be optically connected between all slave lasers in the worst case. The details of the above implementation scheme will be presented in Sec. 7.

Two circular polarization modes in all slave lasers have an identical phase due to injection-locking by the master laser with vertical linear polarization. Those 2 × *M* modes are coherently connected by mutual injection-locking among all slave lasers. Before switching the Hamiltonian Eq. (2) on at *t* = 0, the entire slave laser systems are prepared in 2* ^{M}* linear superposition states,

*R*〉

_{1}and |

*L*〉

_{1}represent the orthogonal polarization states (the right and left circular polarizations) which exist in the slave laser 1. Those 2

*different polarization configurations compete with each other through the cross gain saturation. When the whole slave laser systems reach a steady state condition after the Hamiltonian Eq. (2) is turned on, it is expected that each photon in the*

^{M}*M*slave laser systems is in a specific polarization configuration with a high probability, such as |

*R*〉

_{1}|

*L*〉

_{2}... |

*R*〉

*as shown in Fig. 2, which has the minimum overall loss among 2*

_{M}*polarization configurations. In order to readout the computational result, we can access each partial wave existing in a specific slave laser, as indicated in Fig. 1. Then the computational result can be determined by a majority vote for each slave laser as in Eq. (3).*

^{M}## 3. Theoretical model

We will introduce the coupled rate equations for analyzing our Ising machine, starting with the quantum mechanical Langevin equation for an injection-locked laser [20, 21]. The Heisenberg-Langevin equation for a (non-Hermitian) cavity field operator is described as below,

*μ*is a non-resonant refractive index.

*ω*is the empty cavity resonance frequency of the slave laser, while

_{r}*ω*is the frequency of the injection signal. The operator

*$\tilde{\chi}$*represents the net stimulated emission gain and

_{i}*$\tilde{\chi}$*represents the nonlinear dispersion, both of which include the electron population operator.

_{r}*f̃*is the Langevin noise operator for the electric dipole moment which is associated with random photon emission and absorption by the gain medium.

_{G}*f̂*is the Langevin noise operator for the cavity field, which is associated with the injection signal noise including a vacuum fluctuation.

_{L}*F*

_{0}is a

*c*-number injection signal amplitude.

The net gain operator
$\frac{\omega}{{\mu}^{2}}{\tilde{\chi}}_{i}={\tilde{E}}_{CV}-{\tilde{E}}_{VC}$ is composed of the photon emission rate operator *Ẽ _{CV}* into a lasing mode and absorption rate operator

*Ẽ*. We assume an absorption loss is negligible,

_{VC}*Ẽ*–

_{CV}*Ẽ*≃

_{VC}*Ẽ*, here after. Using the relation between a photon number operator and field operator as

_{CV}*n̂*(

*t*) =

*Â*

^{†}(

*t*)

*Â*(

*t*) in Eq. (6), a quantum mechanical rate equation for photon number operator

*n̂*(

*t*) is derived as follows [20, 23]:

*F̂*(

_{n}*t*) for the photon number operator satisfies the following two-time correlation function [23]:

*κ*is a squeezing parameter and

*κ*< 1 (or

*κ*> 1) represents the amplitude (or phase) squeezed injection signal [20].

Quantum mechanical rate equation for the total electron number operator *Ñ* in a slave laser is [20, 23],

*P*is an average pump rate, $\frac{\tilde{N}}{{\tau}_{sp}}$ is a total spontaneous emission term,

*τ*is a spontaneous emission lifetime. The noise term

_{sp}*F̃*(

_{c}*t*) for the electron number operator is also characterized by [20, 23]

*F̂*and

_{n}*F̃*are negatively correlated, i.e. [20, 23]

_{c}*β*as a fractional coupling efficiency of spontaneous emission into a lasing mode, then a photon emission rate operator into a lasing mode is redefined as ${\tilde{E}}_{CV}=\beta \frac{\tilde{N}}{{\tau}_{sp}}$. A typical value of

*β*is in the range of

*β*= 10

^{−5}∼ 10

^{−4}for VCSELs [24].

When the noise operator *f̂ _{L}* in Eq. (6) represents a vacuum fluctuation, this means the quantum state of the injection signal is a coherent state. In a such a case, the amplitude (or photon number) noise of the injection-locked slave laser output is squeezed to below the standard quantum limit (SQL) under quiet electrical pumping [25], while the phase noise is above the SQL due to random walk phase diffusion [20]. The product of the photon number noise and the phase noise is close to the Heisenberg minimum uncertainty product [20]. That is, the output field of the injection-locked semiconductor laser is close to the minimum uncertainty wavepacket. In such a case, the amplitude and phase noise of the output-field as well as the electron number noise are small compared to their average values and we can safely neglect the quantum noise in the first-order approximation and replace the operators

*Â*(

*t*),

*n̂*(

*t*) and

*N̂*(

*t*) by the corresponding c-numbers. Due to this approximation, we lose the information about the small quantum noise but can still describe the ensemble averaged measurement results for our Ising machine.

We now expand the c-number photon field amplitude as

where*A*

_{0}(

*t*) is the slowly varying amplitude (positive real) and

*ϕ*

_{0}(

*t*) is a slowly varying phase. Substituting Eq. (12) into the classical analog of Eqs. (6) and (7) in which the operators

*Â*and

*n̂*are replaced by the corresponding c-numbers and the noise terms are neglected, we obtain the following three equations:

*E*= 〈

_{CV}*Ẽ*〉 and

_{CV}*χ*= 〈

_{r}*$\tilde{\chi}$*〉. Notice that the third term

_{r}*E*of R. H. S. in Eq. (14) represents the spontaneous emission rate into the lasing mode induced by the quantum noise [22, 23], which we keep in our analysis. The injection-locking bandwidth Eq. (4) is derived from Eq. (15) by requiring |sin[

_{CV}*ϕ*

_{0}(

*t*)]| ≤ 1 and noting

*P*=

_{in}*h̄ω*|

*F*

_{0}|

^{2}and ${P}_{out}=\frac{\omega}{Q}\overline{h}\omega |{A}_{0}(t){|}^{2}$. When

*ω*

_{0}≠

*ω*, the slave laser phase

*ϕ*

_{0}(

*t*) is shifted away ( $-\frac{\pi}{2}<{\varphi}_{0}(t)<\frac{\pi}{2}$) from the master laser phase

*ϕ*= 0. In the following discussions, we assume

_{M}*ω*=

*ω*

_{0}so that the slave laser phase is identically equal to the master laser phase

*ϕ*

_{0}(

*t*) =

*ϕ*= 0.

_{M}Since the phase of a slave laser is locked to that of the dominant vertically polarized master laser signal during an entire computational time, the phases of right and left circular polarization modes in all slave lasers are identical to that of the master laser. Then, the last term of R. H. S. of Eq. (14) can be rewritten as
$2\sqrt{\frac{\omega}{Q}n(t)}{F}_{0}$. If an injection signal is only from the master laser with vertical linear polarization
${F}_{0}=\zeta \sqrt{\frac{\omega}{{Q}_{M}}}\sqrt{{n}_{M}}$, where
$\frac{\omega}{{Q}_{M}}$ is the master laser photon decay rate, *n _{M}* is the average photon number of the master laser and

*ζ*is an amplitude attenuation coefficient including a factor of $\frac{1}{\sqrt{2}}$ for the projection factor form the vertical linear polarization to the right (or left) circular polarization. For simplicity, we assume that the cavity quality factor of the master laser

*Q*is equal to that of the slave lasers,

_{M}*Q*=

_{M}*Q*. Then Eq. (14) can be simplified as,

We now consider two sets of coupled rate equations for the right and left circularly polarized modes of each slave laser. In a proposed injection-locked laser system before the Hamiltonian Eq. (2) is turned on (*t* < 0), the rate equations are

*i*designates a specific slave laser. The vertically polarized injection signal from the master laser into the slave laser

*i*achieves the proper initialization, i.e.,

*n*=

_{Ri}*n*and

_{Li}*ϕ*=

_{Ri}*ϕ*= 0.

_{Li}The right and left circularly polarized modes couple to the conduction electrons with opposite spin angular momenta (
${J}_{e}=\pm \frac{1}{2}$) and the valence holes with opposite orbital and spin angular momenta (
${J}_{h}=\pm \frac{3}{2}$) due to the selection rule [26]. We assume the active quantum well layer is doped with acceptor impurities heavily, so that the gain is governed solely by the minority carrier, conduction electrons. If we denote the electron number of two spin components by *N _{Ri}* and

*N*, we have

_{Li}*τ*is the spin relaxation time in the conduction band. Note that a short

_{spin}*τ*contributes to equalize the populations of two spin components. Our proposed Ising machine relies on this very fast spin relaxation process in VCSELs at room temperature. When

_{spin}*τ*is much shorter than the spontaneous emission lifetime

_{spin}*τ*and the stimulated emission lifetime

_{sp}*τ*[20], Eqs. (19) and (20) are reduced to the single rate equation for the total electron number

_{st}*N*=

_{i}*N*+

_{Ri}*N*:

_{Li}*N*=

_{Ri}*N*, it follows that

_{Li}*E*=

_{CVRi}*E*=

_{CVLi}*E*.

_{CVi}At *t* = 0, we inject a small amount of horizontally polarized signal to the slave laser *i* from the master laser to implement the Zeeman term *λ _{i}* in Eq. (2). At

*t*= 0, we also inject a small amount of horizontally polarized signal from a slave laser

*i*to a slave laser

*j*and vice versa to implement the Ising interaction term

*J*in Eq. (2). Then the rate equations for

_{ij}*n*and

_{Ri}*n*after

_{Li}*t*= 0 are given as follows,

*η*for the horizontally polarized injection signal from the master laser. The mutual injection term is optically implemented by inserting a horizontal linear polarizer, attenuator and phase shifter in the optical path between two slave lasers

_{i}*i*and

*j*. The amplitude attenuation coefficient is $\frac{1}{2}{\xi}_{ij}$ for the horizontally polarized signal.

*η*is positive and real when the phase of the horizontal polarization component is advanced by $\frac{\pi}{2}$ with respect to the vertical polarization component. In this case, the overall injection signal from the master laser is in a left circular elliptic polarization.

_{i}*ξ*is positive and real when the two slave lasers are connected by

_{ij}*π*phase difference. The detailed explanation for implementation of the Ising model will be presented in Sec. 7.

## 4. Mapping protocol

When the system reaches a steady state condition after *η _{i}* and

*ξ*are switched on at

_{ij}*t*= 0, $\frac{d}{dt}{n}_{Ri}=\frac{d}{dt}{n}_{Li}=0$ and $\frac{d}{dt}{N}_{i}=0$ should hold in Eqs. (21)–(23). Then, we can solve for

*E*for the sum of Eqs. (22) and (23) as follows,

_{CVi}*E*, the third terms of R.H.S. in Eqs. (22) and (23), which is much smaller than the stimulated emission term

_{CVi}*E*(or

_{CVi}n_{Ri}*E*), the second terms of R. H. S. in Eqs. (22) and (23) where

_{CVi}n_{Li}*n*and

_{Ri}*n*are of the order of 10

_{Li}^{4}∼ 10

^{5}.

We expect that the coherently coupled slave laser systems oscillate with a polarization configuration which minimizes the overall loss by optimizing the polarization configuration (*σ*
_{1}, *σ*
_{2}, ... *σ _{M}*). The minimized overall loss is identical to the threshold gain,
${\mathrm{\Sigma}}_{i=1}^{M}{E}_{CVi}$. In a standard laser oscillator, there are numerous cavity modes with different eigen-frequencies within a gain bandwidth, so that those multi-modes can oscillate simultaneously if a gain medium is inhomogeneously broadened. However, if we differentiate the loss rates among those cavity modes, for instance by placing a mode (or frequency) selective element inside a cavity and a gain medium is homogeneously broadened, a particular single cavity mode with a minimum loss oscillates alone via cross-gain saturation. We expect the same single mode oscillation is realized in the mutually injection-locked laser systems, which is indeed confirmed by the numerical simulations discussed later.

The first and second terms of R.H.S. of Eq. (24) are almost independent of the polarization configurations, so that the spontaneously selected polarization configuration by the single lasing mode in the entire injection-locked laser systems is expected to minimize

*n*=

_{Ti}*n*+

_{Ri}*n*.

_{Li}If we define an effective Ising spin by ${\sigma}_{iz}=\frac{\sqrt{{n}_{Ri}}-\sqrt{{n}_{Li}}}{\sqrt{{n}_{Ti}}}$ in Eq. (25), the proposed systems oscillate with the polarization configuration that minimizes the quantity

*σ*≤ 1. If we interpret ${\eta}_{i}\frac{\sqrt{{n}_{M}}}{\sqrt{{n}_{Ti}}}$ as a Zeeman term

_{iz}*λ*and

_{i}*ξ*as an Ising interaction term

_{ij}*J*, it is concluded the proposed injection-locked laser systems can find the ground state of an Ising model Eq. (2). When

_{ij}*n*=

_{M}*n*, the two parameters

_{Ti}*η*and

_{i}*ξ*in the rate equations Eqs. (22) and (23) are determined by the relations: where

_{ij}*α*and

*α*′ are extra attenuation parameters that are chosen as small quantities (0 <

*α*≪ 1, 0 <

*α*′ ≪ 1) to ensure the stable operation of a whole system.

In the numerical simulation below (Figs. 3–5), we have used these two parameters *α* and *α*′ which determine the attenuation coefficients for the two injection signals into the slave laser, one from the master laser and the other from other slave lasers. We numerically found that *α*′ must be equal to 2*α* in order to obtain the correct ground state.

## 5. Numerical simulation results

In order to confirm the validity of the proposed Ising machine based on mutually injection-locked laser systems, we perform a series of numerical simulations directly using Eqs. (21), (22) and (23). For the numerical integration of Eqs. (21), (22) and (23), we used the fourth-order Runge-Kutta (RK4) method. The time step is fixed at Δ*t* = 10^{−14}
*sec*, which is sufficiently short compared to any time constants appearing in Eqs. (21), (22) and (23). For a certain parameter range {*J _{ij}*,

*λ*,

_{i}*M*}, we obtained the correct ground states. Since we neglect the quantum noise terms in the quantum mechanical rate equations Eqs. (7) and (9), the numerical results shown in this section are for the ensemble averaged quantities. These numerical simulation results are shown in Figs. 3(a)–3(d) and Figs. 4(a)–4(d). The site number

*M*is varied from

*M*= 2 to

*M*= 10, where the Ising model parameters

*λ*and

_{i}*J*are chosen so that the given Ising model has local minima separated from the ground state by a macroscopic Hamming distance (number of different spins between two spin configurations). We got the right answers for all cases from

_{ij}*M*= 2 to

*M*= 10.

In Fig. 3, the average photon numbers in two polarization modes are equal, *n _{Ri}* =

*n*, at

_{Li}*t*= 0, due to the initialization induced by the vertically polarized injection-locking signal from the master laser. After a short time delay less than 1

*nsec*, the average photon numbers in two polarization modes depart with each other due to the effect of mutual injection locking among slave lasers. The system reaches a steady state with a delay time of a few

*nsec*. The physical origin for this delay time is the inverse of locking-bandwidth and will be discussed elsewhere [27]. We also notice that the average electron numbers

*N*shown in Fig. 4 decrease and reach the steady state values in the same time delay for all slave lasers. This fact confirms our claim that the proposed Ising machine spontaneously selects the specific polarization configuration which minimizes the overall loss and thus achieves the minimum threshold gain Σ

_{i}

_{i}*E*.

_{CVi}We notice that the transient time to reach a steady state does not show a strong dependence on the problem size *M*. It is nearly constant and on the order of a few *nsec* for *M* = 2 to *M* = 10. However, we do not confirm yet if this encouraging result is valid or not for a larger problem size.

The proposed Ising machine also reports a wrong answer for a specific parameter case. Such an example is shown in Fig. 5(a), where the state (*σ*
_{1z}, *σ*
_{2z}, *σ*
_{3z}) = (−1, 1, 1) (= (*L*
_{1}, *R*
_{2}, *R*
_{3})) is chosen by our Ising machine, even though the correct answer is (1, −1, −1). The Zeeman term *η _{i}* in Eqs. (22) and (23) drives the photon numbers

*n*and

_{Ri}*n*to certain directions immediately after the term is switched on. On the other hand, the Ising interaction term

_{Li}*ξ*in Eqs. (22) and (23) does not drive

_{ij}*n*and

_{Ri}*n*initially because of the initial condition

_{Li}*n*=

_{Ri}*n*at

_{Li}*t*= 0. The Zeeman term is more powerful than the Ising term in the beginning of computation. If the system without quantum noise sources is trapped at a local minimum rather than the global minimum (ground state), the system cannot escape from the wrong answer. However the real injection-locked laser systems have inherent quantum noise sources, by which the system has a possibility to escape from the local minimum or the system is dragged toward the global minimum from the beginning. Figure 5(c) shows the numerical results for the same parameter as Fig. 5(a) when we introduce the Poissonian photon number fluctuation to

*n*and

_{Ri}*n*by adding the Gaussian random numbers in the numerical integration of Eqs. (21), (22) and (23). We recover the correct answer with a finite probability as shown in Fig. 5(c). As a comparison we show other examples that the system can find the ground state (−1, 1, 1) for both cases of with and without Poissonian noise in Figs. 5(b) and 5(d), respectively.

_{Li}## 6. Alternative picture of the proposed Ising machine

So far we have described the operational principle of the proposed Ising machine in terms of the probability amplitude modulation and cross gain saturation between right and left circular polarization modes in slave lasers. That is, the probability amplitudes of 2* ^{M}* polarization configurations are initially
$C({\sigma}_{1},\hspace{0.17em}{\sigma}_{2},\dots {\sigma}_{M})=\frac{1}{\sqrt{{2}^{M}}}$ but at the end of computation,

We can alternatively expand the partial wave in a slave laser *i* by the diagonal linear polarization basis, |*D*〉 and |*D̄*〉, as shown in Fig. 2:

*A*

_{D̄i0}(

*t*) and

*ϕ*(

_{D̄i}*t*). Here we assume again

*Q*=

*Q*and $\delta ={\text{tan}}^{-1}\left(\frac{{\eta}_{i}}{\zeta}\right)$. We can numerically integrate those four equations for

_{M}*A*

_{Di}_{0}(

*t*),

*A*

_{D̄i}_{0}(

*t*),

*ϕ*(

_{Di}*t*) and

*ϕ*(

_{D̄i}*t*) together with the corresponding equation of motion for

*N*(

_{i}*t*),

*n*(

_{Di}*t*) =

*A*

_{Di0}(

*t*)

^{2}and

*n*(

_{D̄i}*t*) =

*A*

_{D̄i}_{0}(

*t*)

^{2}.

Figures 6(a) and 6(b) show the average amplitudes *A _{Di}*

_{0}and

*A*

_{D̄i}_{0}, and phase

*ϕ*and

_{Di}*ϕ*in the two diagonal polarization modes for the three sites (

_{D̄i}*M*= 3) problems, calculated by Eqs. (32), (33) and (34). As expected from the picture shown in Fig. 2,

*A*

_{Di0}=

*A*

_{D̄i0}is satisfied at all time (no amplitude modulation between

*A*

_{Di0}and

*A*

_{D̄i0},) but

*ϕ*and

_{Di}*ϕ*depart with each other in a time scale shorter than 1

_{D̄i}*nsec*and this phase modulation reaches the steady state within a few

*nsec*. The slight and simultaneous increase in the amplitudes

*A*

_{Di}_{0}and

*A*

_{D̄i}_{0}results from the overall reduction in the photon loss rate. If we project the complex amplitudes

*A*(

_{Di}*t*) and

*A*(

_{D̄i}*t*) given by Eqs. (30) and (31) onto the photon numbers

*n*and

_{Ri}*n*in circular polarization basis, the two complex amplitudes interfere with each other and we can recover the amplitude modulation and mode competition behaviour between two circular polarization modes:

_{Li}*R*〉 and |

*L*〉 basis. We confirm that the two curves shown in Figs. 6(c) and 6(d) are completely identical. The two complementary pictures, the mode competition between |

*R*〉 and |

*L*〉 states and the mutual interference between |

*D*〉 and |

*D̄*〉 states, can explain the same operational principle of the proposed Ising machine equally well.

## 7. Implementation

In this section we describe how to implement the Ising model Eq. (2) in injection-locked laser systems. As shown in Fig. 1, the output from the master laser is equally split and injected into all slave lasers via an optical isolator. The master laser output is used for two different purposes. One is to initialize all slave lasers in a linear superposition of right and left circular polarization states or the two diagonal linear polarization states, i.e., in a vertical linear polarization, and another is to implement the Zeeman term in the generalized Ising model Eq. (2). The output from each slave laser is also split and injected into the other slave lasers (without an isolator in an optical path). The final polarization state at each site (slave laser) is read out by using a polarization beam splitter (PBS) and dual-photodetectors. The implementation at each computational step is described below:

- Initialization (t<0)
Before introducing the Zeeman term and Ising interaction term, we prepare the initial states in all slave lasers. The frequency, phase and polarization of all slave lasers are identical to those of the master laser. In Eqs. (22) and (23) or Eqs. (32) and (33),

*ζ*determines the amount of vertical linearly polarized injection signal. This initialization step is implemented by an isolator and simple attenuators. - Zeeman term
*λ*(t≥0)_{i}At a time

*t*= 0, we implement the amplitude and phase of the Zeeman term*λ*for each site_{i}*i*. In the rate equations Eqs. (22) and (23) for the photon numbers, the Zeeman term*λ*is represented by the horizontal polarization component of the injection signal from the master laser, which is expressed by the term proportional to_{i}*η*. As shown in Fig. 7(a), the magnitude of_{i}*λ*(|_{i}*λ*|) is implemented by choosing the proper rotation angle_{i}*θ*for the initial vertical linear polarization signal to include the horizontal polarization component. This is easily achieved by inserting a HWP on a path from the master laser to each slave laser. The horizontal polarization component $|H\u3009=\frac{\u0131}{\sqrt{2}}\left(|R\u3009-|L\u3009\right)$ of the master laser signal couples to right and left circular polarization modes of the slave laser with $\frac{\pi}{2}$ and $-\frac{\pi}{2}$ phases, respectively. Notice that the sign in front of the constant*η*is opposite for_{i}*n*and_{Ri}*n*in Eqs. (22) and (23). Here the rotation angle_{Li}*θ*is determined by the amount of required horizontal polarization component |*H*〉. The sign of*η*is minus for negative_{i}*λ*and plus for positive_{i}*λ*in Eq. (27). To implement the polarity of_{i}*λ*, we put a QWP after HWP to introduce $\pm \frac{\pi}{2}$ phase shift. Then the polarization state changes into the right or left ellipsoidal polarization cos_{i}*θ*|*V*〉 ± ýsin*θ*|*H*〉. The QWP rotation angle determines the dominant polarization in a slave laser. When the rotation angle is $\frac{\pi}{2}$,*η*is positive so that the left circular polarization will be dominant in the slave laser_{i}*i*. - Ising interaction term
*J*(t≥0)_{ij}The emitted light from each slave laser is injected to other slave lasers via a horizontal linear polarizer, phase shifter and attenuator as shown in Fig. 7(b). The magnitude of

*J*is implemented by the transmission coefficient of the attenuator. The horizontal polarizer output of the slave laser_{ij}*j*is proportional to $\sqrt{{n}_{Rj}}-\sqrt{{n}_{Lj}}$. This mutual injection signal with horizontal linear polarization couples to the right and left circular polarization modes of the slave laser*i*with equal amplitudes and opposite phases. Notice that the sign in front of the transmission constant*ξ*is opposite for_{ij}*n*and_{Ri}*n*in Eqs. (22) and (23). The sign of_{Li}*J*is implemented by the optical phase shift_{ij}*ϕ*between the two slave lasers._{ij}*ϕ*= 2_{ij}*πl*is chosen for a negative*J*and_{ij}*ϕ*= 2_{ij}*πl*+*π*for a positive*J*, where_{ij}*l*is a positive integer. - Readout
After a steady state condition is reached, the outputs from all slave lasers are input into polarization beam splitters (PBS) to separate the right and left circular polarization components. The photodetectors at each output port of the PBS provides the computational result through a majority vote.

- Signal to noise ratio
Let

where $Q=\u3008{m}_{Ri}-{m}_{Li}\u3009/2\sqrt{\u3008\mathrm{\Delta}{m}_{Ri}^{2}\u3009+\u3008\mathrm{\Delta}{m}_{Li}^{2}\u3009}$ is the Personick Q parameter. The required (peak-to-peak) signal to noise ratio is approximately 4*m*and_{Ri}*m*denote the right and left circularly polarized and time integrated photon numbers. The electrical signal power is proportional to the squared difference of the average photon numbers 〈_{Li}*m*–_{Ri}*m*〉_{Li}^{2}. If the two polarization modes are in independent coherent states, the noise power is proportional to the sum of the variances $\u3008\mathrm{\Delta}{m}_{Ri}^{2}\u3009+\u3008\mathrm{\Delta}{m}_{Li}^{2}\u3009$ where ${m}_{Ri}={\eta}_{D}\left(\frac{\omega}{Q}\right){n}_{Ri}T$, ${m}_{Li}={\eta}_{D}\left(\frac{\omega}{Q}\right){n}_{Li}T$,*η*is an overall detection quantum efficiency and_{D}*T*is an integration time. Due to the above assumption, the slave laser output photon statistics obey independent Poisson distributions, i.e., $\u3008\mathrm{\Delta}{m}_{Ri}^{2}\u3009=\u3008{m}_{Ri}\u3009$ and $\u3008\mathrm{\Delta}{m}_{Li}^{2}\u3009=\u3008{m}_{Li}\u3009$. Since the Poisson distribution with a large average photon number is well approximated by the Gaussian distribution with the same variance, the bit error rate in this detection scheme is given by the complementary error function of the signal-to-noise ratio [28]:*Q*^{2}∼ 144 (21.6*dB*) to achieve a prescribed error rate, for instance,*P*= 10_{e}^{−9}. Then, we have$$4{Q}^{2}=\frac{{\u3008{m}_{Ri}-{m}_{Li}\u3009}^{2}}{\u3008\mathrm{\Delta}{m}_{Ri}^{2}\u3009+\u3008\mathrm{\Delta}{m}_{Li}^{2}\u3009}=\frac{{\eta}_{D}(\frac{\omega}{Q})\u3008{n}_{Ri}\u3009T{(1-R)}^{2}}{1+R}\ge 144,$$where $R=\frac{\u3008{n}_{Li}\u3009}{\u3008{n}_{Ri}\u3009}<1$ is the ratio of the average photon numbers of the minor to major polarization modes. If we use the numerical examples of*η*= 0.01, $\frac{\omega}{Q}={10}^{12}\left(\frac{1}{s}\right)$, 〈_{D}*n*〉 = 10_{Ri}^{4}and*R*= 0.5, the required integration time to satisfy Eq. (38) is*T*= 8.4*psec*, which is much shorter than the computational time. If the photon number noise is above the Poisson limit due to the mode partition noise [29], we can always satisfy Eq. (38) by employing a longer integration time. Note that the overall error rate is only additive with the site number*M*, so that the over all error rate is still bounded by*P*= 10_{e}^{−6}even for a large-scale problem*M*= 1000.

## 8. Conclusion

We proposed the mapping protocol of an Ising model onto injection-locked laser systems. Each site of an Ising model is represented by an injection-locked slave laser, in which the Ising spin (up or down) corresponds to the right or left circular polarization of the partial wave of each individual photon which coherently spreads in the entire slave laser systems. Before a computation will start, all slave lasers are injection-locked by a vertical linearly polarized master laser output, so that the two polarization modes of all slave lasers have identical average photon numbers and phases. The Ising interaction term *J _{ij}* can be implemented by mutual injection between the two slave lasers

*i*and

*j*via a horizontal linear polarizer, phase shifter (0 or

*π*phase) and attenuator. The Zeeman term

*λ*can be implemented by injecting a horizontal linearly polarized master laser signal in addition to the vertical linearly polarized master laser output into the slave laser

_{i}*i*.

We have confirmed numerically that the proposed Ising machine outputs correct answers for most cases even if the quantum noise source is neglected. For some other cases, the proposed system outputs wrong answers due to the trapping problem in local minimum states. It is found that the Poissonian photon number noise can make the system to escape from the local minimum state and to find the global minimum state. We have analyzed the proposed Ising machine in terms of complementary diagonal linear polarization basis, in which the amplitudes of the two polarization modes |*D*〉 and |*D̄*〉 are identical but the phases depart with each other. When such phase modulated diagonal linear polarization modes are projected onto the two circular polarization modes |*R*〉 and |*L*〉, we recover the identical results and confirm the two alternative pictures explain the same operational principle of the proposed Ising machine equally well.

A time delay to reach a steady state condition after the Ising and Zeeman energy terms are switched on can be considered as a lower bound of computational time. The transient time does not show a strong dependence on the problem size *M* according to numerical simulations. It is on the order of a few *nsec* from *M* = 2 to *M* = 10. However, even an efficient physical simulator is often limited by getting trapped in one of many metastable states (local minima) in the case of NP-complete problems when the problem size becomes large [30]. The relation of the computational power of the proposed Ising machine to the NP completeness in the limit of large problem size will be a subject of future study.

## 9. Summary

The proposed Ising machine can be potentially implemented by semiconductor VCSELs, ECDL and linear optical components even for a very large-scale problem *M* ∼ 1000. Even though a preliminary numerical simulation result is encouraging in terms of the computational time vs. problem size *M*, its asymptotic behaviour for NP-complete problems with large *M* should be carefully studied in a future publication. After completing this work, we were informed of the quantum optical simulation of Ising models using atom-cavity QED systems proposed in [31].

## Acknowledgments

The authors wish to thank a useful discussion with Kai Wen. This work is supported by the JSPS through its FIRST program, Navy/SPAWAR Grant and the Special Coordination Funds for Promoting Science and Technology.

## References and links

**1. **M. R. Garey and D. S. Johnson, *Computers and Intractability: A Guide to the Theory of NP-Completeness* (W. H. Freeman and Company, 1979).

**2. **K. Binder and A. A. Young, “Spin glasses: experimental facts, theoretical concepts, and open questions,” Rev. Mod. Phys. **58**, 801–976 (1986). [CrossRef]

**3. **V. Dotsenko, *Introduction to the Replica Theory of Disordered Statistical Systems* (Cambridge University Press, 2000). [CrossRef]

**4. **H. Nishimori, *Statistical Physics of Spin Glasses and Information Processing: An Introduction* (Oxford University Press, 2001). [CrossRef]

**5. **A. Das and B. K. Chakrabarti, “Colloquium: quantum annealing and analog quantum computation,” Rev. Mod. Phys. **80**, 1061–1081 (2008). [CrossRef]

**6. **P. Ray, B. K. Chakrabarti, and A. Chakrabarti, “Sherrington–Kirkpatrick model in a transverse field: absence of replica symmetry breaking due to quantum fluctuations,” Phys. Rev. B **39**, 11828–11832 (1989). [CrossRef]

**7. **B. Appoloni, C. Carvalho, and D. de Falco, “Quantum stochastic optimization,” Stochastic Proc. Appl. **33**, 233–244 (1989). [CrossRef]

**8. **R. Martonak, G. E. Santoro, and E. Tosatti, “Quantum annealing by the path-integral Monte Carlo method: the two-dimensional random Ising model,” Phys. Rev. B **66**, 094203 (2002). [CrossRef]

**9. **G. E. Santoro, R. Martonak, E. Tosatti, and R. Car, “Theory of quantum annealing of an Ising spin glass,” Science **295**(5564), 2427–2430 (2002). [CrossRef] [PubMed]

**10. **G. E. Santoro and E. Tosatti, “Quantum to classical and back,” Nat. Phys. **3**, 593–594 (2007). [CrossRef]

**11. **R. D. Somma and C. D. Batista, “Quantum approach to classical statistical mechanics,” Phys. Rev. Lett. **99**, 030603 (2007). [CrossRef] [PubMed]

**12. **J. Brooke, D. Bitko, T. F. Rosenbau, and G. Aeppli, “Quantum annealing of a disordered magnet,” Science **284**(5415), 779–781 (1999). [CrossRef] [PubMed]

**13. **G. Aeppli and T. F. Rosenbaum, in *Quantum Annealing and Related Optimization Methods*, A. Das and B. K. Chakrabarti, eds. (Springer Verlag, 2005).

**14. **M. Steffen, W. van Dam, T. Hogg, G. Breyta, and I. Chuang, “Experimental implementation of an adiabatic quantum optimization algorithm,” Phys. Rev. Lett. **90**, 067903 (2003). [CrossRef] [PubMed]

**15. **T. Byrnes, K. Yan, and Y. Yamamoto, “Optimization using Bose-Einstein condensation and measurement-feedback circuits,” arXiv:0909.2530v2.

**16. **K. Yan, T. Byrnes, and Y. Yamamoto, “Kinetic Monte Carlo study of accelerated optimization problem search using Bose-Einstein condensates,” Prog. Inform. **8**, 1–9 (2011).

**17. **F. Barahona, “On the computational complexity of Ising spin glass models,” J. Phys. A: Math. Gen. **15**, 3241–3253 (1982). [CrossRef]

**18. **S. Kobayashi and T. Kimura, “Injection locking in AIGaAs semiconductor laser,” IEEE J. Quantum Electron. **17**(5), 681–689 (1981). [CrossRef]

**19. **S. Kobayashi, Y. Yamamoto, and T. Kimura, “Optical FM signal amplification and FM noise reduction in an injection locked AlGaAs semiconductor laser,” Electron. Lett. **17**(22), 849–851 (1981). [CrossRef]

**20. **L. Gillner, G. Bjork, and Y. Yamamoto, “Quantum noise properties of an injection-locked laser oscillator with pump-noise suppression and squeezed injection,” Phys. Rev A **41**(9), 5053–5065 (1990). [CrossRef] [PubMed]

**21. **H. A. Haus and Y. Yamamoto, “Quantum noise of an injection-locked laser oscillator,” Phys. Rev. A **29**, 1261–1274 (1984). [CrossRef]

**22. **M. Surgent, M. O. Scully, and W. E. Lamb, *Laser Physics* (Westview Press, 1978), Chap. 20, pp. 331–335.

**23. **H. Haug, “Quantum-mechanical rate equations for semiconductor lasers,” Phys. Rev.184, 338–348 (1969). [CrossRef]

**24. **S. F. Yu, *Analysis and Design of Vertical Cavity Surface Emitting Laser* (Wiley-Interscience, 2003), Chap. 8. [CrossRef]

**25. **Y. Yamamoto, S. Machida, and O. Neilsson, “Amplitude squeezing in a pump-noise-suppressed laser oscillator,” Phys. Rev. A **34**, 4025–4042 (1986). [CrossRef] [PubMed]

**26. **C. Weisbuch and B. Vinter, *Quantum Semiconductor Structures: Fundamentals and Applications* (Academic Press, 1991), Chap. 3, pp. 65–69.

**27. **K. Takata, S. Utsunomiya, and Y. Yamamoto, “Transient time of an Ising machine based on injection-locked lasers: contribution of locking bandwidth and Zeeman component,” (to be submitted).

**28. **S. D. Personick, “Receiver design for digital fiber optic communication systems, I,” Bell Syst. Tech. J. **52**(6), 843–874 (1973).

**29. **S. Inoue, H. Ohzu, S. Machida, and Y. Yamamoto, “Quantum correlation between longitudinal mode intensities in a multimode squeezed semiconductor laser,” Phys. Rev. A **46**, 2757–2765 (1992). [CrossRef] [PubMed]

**30. **A. P. Young, S. Knysh, and V. N. Smelyanskiy, “First-order phase transition in the quantum adiabatic algorithm,” Phys. Rev. Lett. **104**, 020502 (2010). [CrossRef] [PubMed]

**31. **Y. Zhang, Y. Xia, Z. Man, and G. Guo, “Simulation of the Ising model, memory for Bell states and generation of four-atom entangled states in cavity QED,” Sci. China, Ser. G **52**, 700–707 (2009). [CrossRef]