Transient time of an Ising machine based on injection-locked laser network

We numerically study the dynamics and frequency response of the recently proposed Ising machine based on the polarization degrees of freedom of an injection-locked laser network (Utsunomiya et al 2011 Opt. Express 19 18091). We simulate various anti-ferromagnetic Ising problems, including the ones with symmetric Ising and Zeeman coefficients, which enable us to study the problem size up to M = 1000. Transient time, to reach a steady-state polarization configuration after a given Ising problem is mapped onto the system, is inversely proportional to the locking bandwidth and does not scale exponentially with the problem size. In the Fourier analysis with first-order linearization approximation, we find that the cut-off frequency of a system's response is almost identical to the locking bandwidth, which supports the time-domain analysis. It is also shown that the Zeeman term, which is created by the horizontally polarized injection signal from the master laser, serves as an initial driving force on the system and contributes to the transient time in addition to the inverse locking bandwidth.


Introduction
Since Feynman pointed out the exponential growth of physical resources to simulate a quantum system with classical computing machines [1], many studies have focused on realizing quantum simulators or analogue quantum computers that efficiently simulate a quantum system [2,3]. Distinct from a digital quantum simulator, which essentially operates a phase estimation algorithm in a digital quantum computer [4], an analogue quantum simulator can directly map a quantum system of interest onto another quantum system with better controllability [5][6][7][8][9]. A quantum simulator has a much simpler structure with reduced resources compared to a faulttolerant digital quantum computer with a quantum error correcting code. Nevertheless, such machines are expected to solve hard mathematical, chemical and physical problems which are intractable with current computational power.
An Ising problem is a classic example of an NP-complete problem [10]. It is the searching problem for the ground state of an Ising Hamiltonian which is a relatively simple theoretical model but elucidates the orders in various magnetic materials, frustrated spin lattices and spin glasses [11]. Even though the two-dimensional (2D) Ising model can be solved in a polynomial time with a deterministic Turing machine, solving the 2D Ising model with external magnetic field and the 3D Ising model is known to be NPcomplete [10]. Other NP-complete problems, such as the graph partitioning problem and MAX-CUT problem, can be mapped on the Ising model (1) with a polynomial time [11]. Therefore, developing an efficient Ising machine is a relevant challenge in quantum information science.
3 As a candidate for such an Ising machine, a system of trapped bosonic particles under Bose-Einstein condensation with a measurement-feedback circuit was proposed [12]. A measurement-feedback circuit introduces non-local interaction coefficients J i j between any two sites. However, it was numerically found that this simulator cannot prevent the exponential growth of the computational time with the problem size [13], possibly because a measurementfeedback circuit cannot create quantum coherence between the different Ising sites.
In order to overcome these difficulties, we have recently proposed a new Ising machine based on an optical coherent feedback circuit [14]. It consists of a master laser, M slave lasers and linear optical components to control the amplitude and phase of optical signals in the laser system. The master laser locks the optical frequency and phase of all slave lasers [15][16][17][18][19] to keep quantum coherence in the whole system and also provides the Hamiltonian with the effective magnetic field (Zeeman term). Each slave laser represents an Ising spin in terms of its circular polarization state. The Ising interaction term can be implemented by mutual injection between two slave lasers [20]. We have numerically found that a transient time, which is one fundamental time scale to determine a computational time of the Ising machine, does not scale exponentially with the problem size for particular Ising problems. However, we have also found that there are certain problems for which this Ising machine can find a correct answer only statistically with the help of quantum noise [14].
In this paper, we focus on the scaling law of a transient time for the problems for which a correct answer is obtained deterministically. We numerically show that the order of transient time is dependent not so much on the problem size M, but determined by the locking bandwidth of the slave laser. We take two complementary approaches to reach the conclusion. In the first approach, we simulate the time evolution of the system with coupled rate equations. In the second approach, we calculate the frequency response to the small optical modulation signal using the linearization method [21][22][23]. It shows that the cut-off frequency of the response function is almost identical to the locking bandwidth. The fact has already been demonstrated experimentally for the single injection-locked laser [24], and we find that the response of the M mutually coupled slave lasers is also governed by this single parameter. We also note that the small initial kick provided by the Zeeman component introduces an additional delay on the transient time.

Operator Langevin equations for an injection-locked semiconductor laser
The model of an injection-locked slave laser is shown in figure 1. The operator Langevin equation for the cavity internal fieldÂ(t) is [22] where F 0 andf (t) are the slowly varying c-number average amplitude and noise operator of an injection signal. The optical angular frequency of the injection signal is denoted by ω, while the (empty) cavity resonance angular frequency is ω r . The total Q of the slave laser cavity is split r (output field) into the external (mirror) loss Q e and the internal loss Q 0 according to µ is the (non-resonant) refractive index andχ =χ r + iχ i is the electric susceptibility operator. The net stimulated emission gain is expressed in terms of the imaginary part ofχ: whereẼ C V andẼ V C are the stimulated emission and absorption rate operators. The (active) cavity resonant frequency or self-oscillation frequency of the slave laser is expressed in terms of the real part ofχ: where the second term represents the frequency pulling and pushing effect induced by the anomalous dispersion of the gain medium [25]. The noise operatorsf G (t),f L (t) andf (t)e −iωt are associated with the stimulated emission/absorption, internal loss and injection signal fluctuations, respectively. We use the tilde to denote the electronic operators and the hat to express the photonic operators. We assume that the active layer of a slave laser is heavily doped with acceptor impurities so that the gain is determined solely by an injected electron number. The operator Langevin equation for the electron numberÑ C (t) is [22] d dtÑ where p = I /e is the pumping rate, I is an injection current, e is an electronic charge, τ sp is the spontaneous emission lifetime andn(t) =Â † (t)Â(t) is the photon number operator.˜ p ,˜ sp and are the noise operators associated with the pump rate, the spontaneous emission rate into all non-lasing modes and the stimulated emission/absorption rates into a lasing mode, respectively.

Rate equations for the Ising machine based on an injection-locked laser network
Because of the simplification about the phase difference between the internal fields of the slave lasers and the injection signals, we can neglect the phase dynamics of the slave lasers (see the appendix) and describe the system performance only by the rate equations for the photon numbers and the electron numbers. In our previous paper [14], we showed that a general Ising Hamiltonian, is simulated by the experimental system in figure 2. The system is well described by the rate equations for the average number of right and left circularly polarized photons (n R and n L ) and the average electron number (N C ). Deriving these rate equations from the operator Langevin equations (2) and (6) requires the following procedure in summary [14,26]: (i) neglect the absorption loss asẼ C V −Ẽ V C Ẽ C V , (ii) derive the operator rate equation for the photon number operator with the identity (d/dt)n = [(d/dt)Â † ]Â +Â † [(d/dt)Â] and (iii) neglect the small quantum noise terms and approximate the operatorsÂ,N C ,Ẽ C V andn by the ensemble-averaged c-number values. Here, note that the Ising and Zeeman terms in the Hamiltonian (7) are implemented by the horizontally polarized component in the mutually injected slave signals and the master signal, respectively.
After these steps, the equations for the ith slave laser including injected signals from the master laser and other slave lasers are given by [14] d dt where a specific slave laser is denoted by an index i = 1, 2, . . . , M, and ω/Q = ω/Q e is the cavity photon decay rate. ζ and η i are the optical attenuation coefficients for a vertically polarized signal and a horizontally polarized signal from the master laser, and ξ i j is that for an injected signal from another slave laser j ( = i). E C V i = Ẽ C V i describes the average stimulated emission rate in the ith slave laser and is related to the average electron number N Ci by We define the common attenuation coefficient α, which relates the Ising and Zeeman coefficients in (7) to the optical attenuation coefficients in (8) and (9), as Here, our previous work [14] numerically showed that the extra factor 2 in (13) is necessary for simulating the Hamiltonian correctly. That paper also introduces the Ising spins expressed by the amplitude difference as σ i = ( √ n Ri − √ n Li )/ √ n Ri + n Li where √ n Ri + n Li is the total field amplitude in a slave laser. However, in the actual computation, we can read them out by simple majority vote after a steady-state condition is reached: and expect that the {σ i } gives the ground state of the given Hamiltonian (7). If a specific spin configuration is a ground state, a completely opposite spin configuration is also a ground state in (1). Thus the Zeeman term is introduced in (7) in order to lift this trivial degeneracy of the Ising Hamiltonian (1).

Locking bandwidths for lasing polarization modes
The rate equations used here treat the right and left circularly polarized modes separately. Thus, it is appropriate to define the locking bandwidths for the two modes separately. A locking bandwidth is expressed by the average injection signal amplitude F inj and the internal field amplitude A int in the slave laser as (see the appendix) Using (15) and the injection signal terms in (8) and (9), we can define the locking bandwidth for each polarization mode as Here, note that a locking bandwidth can be defined only for lasing states (n(t) > 1/β). Thus, we do not consider the locking bandwidths for the polarization states that are suppressed in the time evolution of the system later on.

Transient time versus problem size
In this section, we show the numerical results of a time-domain analysis of (23)- (25). Since these equations neglect the quantum noise sources, the results presented here represent the ensemble averaged dynamics. This means that, if we repeat the same computation process many times, then the ensemble averaged results are reduced to the results presented in this section. For a particular problem attacked here, however, the actual single-shot result including the quantum noise effect deviates only slightly from the ensemble averaged result. We used the standard fourth-order Runge-Kutta method for the numerical integration. The time step t = 10 −14 , 10 −13 or 10 −12 s is selected considering the transient time scale and accuracy necessary for each simulation. For the sake of initialization (t < 0), we assume that only the vertically polarized signal from the master laser is injected into the slave lasers. Thus, the variables of the initial state (t = 0) are given as the solutions of the steady-state equations The numerical parameters used here are as follows: β = 10 −4 , τ sp = 10 −9 s, ω/Q = ω/Q e = 10 12 s −1 (except for the (ω/Q) −1 dependence analysis), ω/Q 0 = 0 and ζ = 1/200. The former three are typical of semiconductor lasers. And the value of ζ ensures the locking bandwidth of about 1 GHz in the whole system with the photon lifetime (ω/Q) −1 = 10 −12 s. by the interaction term j 1 2 ξ i j ( √ n R j − √ n L j ), and the system ends up in the correct ground state {σ 1 , σ 2 } = {−1, +1}. In this way, we got the correct answers to all the problems treated in this paper. As shown in figure 3(b), the carrier numbers are all decreased as time elapses, which confirms the reduction of the total threshold gain i E C V i = i β N Ci /τ sp . Furthermore, the locking bandwidths of the lasing modes are also time varying due to the transition of the photon numbers in this system. Figure 3(c) shows that all of the lasing polarization modes have larger locking bandwidths in the final state than in the initial state. This stems from the increased injection signal power due to constructive interference among the injection signals from the master laser and the slave lasers. Note that the steady-state condition is reached at a time t T ∼ 1 ns, while the minimum injection locking bandwidth in the initial state is ω L /2π ∼ 1 GHz. Thus, the product of t T and ω L /2π has the order of one.

A symmetric case
In order to study the transient time for a very large problem size up to M = 1000, we now introduce a symmetric Ising model. We treat the Ising problems whose coefficients {J i j } and {λ i } are symmetric among the slave lasers with odd and even index numbers λ i = λ odd (i = odd number), λ even (i = even number).
Using (12), (13), (21) and (22), we can reduce the 3M rate equations (8)-(10) to the six rate equations for the slave lasers with even and odd number: where the parity sign {P,P} is either {odd, even} or {even, odd} depending on either i = odd or i = even. M P1 and M P2 describe the number of slave lasers with the same parity P andP, other than the slave i: where and mean floor and ceiling functions. Also, from (15), (23), (24), (26) and (27), the locking bandwidths for this case are given by Now we show the numerical simulation results on the symmetric Ising problems. Here we consider the two problem sets: .
Set 1 is a standard anti-ferromagnetic Ising problem set with constant Ising and Zeeman coefficients, while set 2 is a harder one in terms of the potential energy landscape (PEL). The energy difference between the ground state {σ 1 , σ 2 , σ 3 , . . .} = {−1, +1, −1, . . .} and the totally opposite spin configuration {+1, −1, +1, . . .}, which is a metastable state, is M|λ odd − λ even |. Therefore, in the problem set 1 the energy difference increases linearly with M, although it takes the constant value of 1/10 in the problem set 2. Thus, the energy density of metastable states is higher in problem set 2 than in problem set 1. Note that a local minimum exists with the farthest Hamming distance from the ground state in these two problem sets when i is an even number. When i is an odd number, there are M−1 2 local minima with the Hamming distance M − 1 from the ground state. From the experimental point of view, problem set 1 corresponds to the master laser with unlimited power while problem set 2 describes the master laser with fixed laser power.
Note that this machine starts with the coherent linear superposition of all polarization configurations in the PEL, with equal probability amplitudes. In addition, the Hamiltonian has a small Zeeman term to avoid degeneracy between the states with completely opposite spin configurations. Therefore, this machine can solve ferromagnetic problems very easily and it is therefore important whether it can find a correct ground state slightly below many metastable states, in anti-ferromagnetic problems.  figure 5(c), we see that the locking bandwidths in the steady state are larger than those of the initial states as in figure 3(c). This means that the system stabilizes itself by making the mutual injection signal maximum in each slave laser. We use the locking bandwidth of the initial state f as a reference later on.
The logarithmic plot of transient time t T versus photon lifetime (ω/Q) −1 is shown in figure 6. Here we added the inverse locking bandwidth 1/ f versus (ω/Q) −1 , where f (t = 0) is the initial locking bandwidth of the slowest mode. The transient time shows the linear dependence on ∼ 1/ f . Therefore, we conclude that the transient time of the Ising machine is mainly determined by the locking bandwidth. We used the fixed normalized pump rate I /I th ∼ 3, where I th is the threshold current and is inversely proportional to (ω/Q) −1 . In this case, the photon number in each slave laser is fixed at about 2 × 10 4 . We also did the same simulation with the fixed pumping current I = 4.8 mA; however, the curves we got are almost the same as figure 6, except that t T of set 1 gets slightly larger.

Small signal analysis of rate equations
In order to analyze the frequency-domain response of the injection-locked laser system, we introduce a small modulation into the common optical attenuation coefficient α as where α 0 and α(t) are the dc attenuation and the sinusoidal modulation signal at modulation frequency . Here, we assume that n Ri , n Li and E C V i can be expanded in the same form as α: We substitute (31)-(34) into the rate equations for the symmetric case (23)- (25) and keep the first-order terms. Note that any square-root variables are linearized by the first-order approximation . After subtracting all the dc terms, we have the linearized equations for the fluctuating variables where the coefficients A RP to C EP are defined by where we apply (d/dt) → i to the rhs of (38), (45) and (51) and have A RP ( ), B LP ( ) and C EP ( ).

Numerical results
Next, we discuss the numerical results of the frequency-domain analysis using (55)-(57) with the coefficients (38)-(51). Our goal is to investigate the cut-off frequency of the response functions (52)-(54). From the viewpoint of comparison with the overall transient time, it is appropriate to calculate the response functions at the initial state (t = 0). However, in order to examine the effect of the problem size and the Zeeman term, we need to include all injection signals in this analysis. Also, we have to get the dc terms {n RP0 , n LP0 , E C V P0 } to solve (55)-(57).
Here, we expect to obtain dc terms by solving (23)-(25) under the steady-state condition with a small α 0 . Quantitatively, this corresponds to reducing α 0 toward zero with its value kept finite. And qualitatively, this means that we change the values of α 0 and M under the condition that the photon numbers of both modes in the steady state are hardly varied from those for α 0 = 0. By the time-dependent analysis of (23)- (25), which is not shown here, we have found that 10 −6 α 0 10 −4 is a good range from the qualitative standpoint above. Figures 8(a) and (b) show an example of the normalized response function for the photon number of the right circularly polarized mode (H Ri ) and the spontaneous emission rate (H Ei ) as a function of modulation frequency /2π. We found that the response of the left circularly polarized mode almost overlaps that of the right circularly polarized mode in all the data given in this paper; thus we drop it. As shown previously by the experiment with a single injection-locked  laser [24], the photon number response ( figure 8(a)) has no relaxation oscillation peaks. This is consistent with the curves in figure 5(a) which do not display any oscillatory behaviors. On the other hand, that of spontaneous emission rate ( figure 8 (b)) features the relaxation oscillation peak and also a larger cut-off frequency. Figures 9 and 10 show the dependence of the cut-off frequency f C of the photon number and the spontaneous emission rate on M and α 0 , respectively. These graphs are plotted for problem set 1. Those for problem set 2 are almost identical to figures 9 and 10. In both figures, the cut-off frequency f C for the photon number response function almost equals the corresponding locking bandwidth when the variable (M or α 0 ) is low. This is also consistent with the previous single-laser experiment [24]. And the cut-off frequencies in these graphs show monotonically increasing behavior. Therefore, we reconfirm that the locking bandwidth determines the transient speed of the system. In addition, we can see that the Zeeman term is not related directly to the cut-off frequency f C and also expect that f C is dependent on α 0 M, meaning that the contributions of α and M to f C are equivalent in the range of the first-order approximation.
We finally show that the Zeeman term affects not f C but the magnitude of the response function. The low-frequency response function |H ( /2π = 1 Hz)| versus M for problem sets 1 and 2 is plotted in figures 11(a) and (b), respectively. Under the condition of problem set 1 (in figure 11(a)), the response functions exhibit flat behavior up to M ∼ 100. However, in problem set 2 ( figure 11(b)), the response functions of photon number are inversely proportional to M, which directly reflects the decreasing function of the Zeeman term on M. It follows that the Zeeman term plays a crucial role in driving the system at the initial stage. If this initial trigger is too weak, then a transient time must be considerably larger than 1/ f , or the system may even fail to ignite the significant change in the photon numbers n R and n L to find the correct ground state, in the range of this rate equation analysis.

Discussion
Here we discuss some possible limitations to the locking bandwidth (computational speed) and problem size (number of Ising sites) of the Ising machine based on an injection-locked laser network. In order to ensure a good value of the locking bandwidth, we need a certain amount of vertically polarized injection signal ζ √ n M . However, such an injection signal tends to equalize the photon numbers of two modes in each slave laser, thus hampering the switching between circular polarization modes, i.e. spin flip. This means ζ cannot be too small, in order to ensure a large enough locking bandwidth, and α has to be large enough to get the transient started.
α also determines the upper bound of the problem size M. Laser output from a slave laser is distributed into the other slave lasers as mutual injection signal. Here, the sum of distributed signal powers cannot exceed the total output from the slave laser; thus we can get the condition j =i (ξ i j ) 2 = j =i (2J i j α) 2 < 1. Note also that large α and M, i.e. a large injection signal power, can make the injection locking unstable [19]. Although we employed the systematic Ising model to study the problem size with M up to 1000, our time evolution analysis is restricted to a parameter range where stable solutions are obtained. The boundary between stable and unstable operations of this computing system will be discussed in a future publication.

Conclusion
We have theoretically studied the response time of the Ising machine based on an injectionlocked laser network. The response time is primarily determined by the inverse of the locking bandwidth of slave lasers. We have not observed a strong dependence of the response time on the problem size M, which suggests the interesting possibility that the proposed machine may solve an NP-complete Ising problem with a polynomial time. This statement must be examined by a further benchmarking test and our initial result for an NP-complete MAX-CUT problem in cubic graphs seems promising [27]. It is well known that all NP problems can be mapped to any NP-complete problem; thus this machine is expected to solve all NP problems in polynomial time eventually. The implementation of spin glass problems and general spin models, including quantum mechanical Heisenberg coupling, on this machine is also an interesting direction. the fractional spontaneous emission coupling efficiency β, the population inversion parameter n sp and the net gain coefficient G by The output field amplituder (t) is related to the internal fieldÂ(t) and the injection signal [F 0 +f (t)]e −iωt through [22] The output field is also linearized aŝ The steady-state solutions for the internal photon field are obtained from the nonfluctuating real and imaginary parts of (2): 14) The self-oscillation frequency ω 0 given by (5) is modified by using (A.7), (A.10) and (A.14) Here ω r0 is the self-oscillation frequency of the slave laser without an injection signal. The second term of (A.15) represents the shift of the oscillation frequency associated with the reduced threshold gain G − ω Q < 0 induced by an injection signal F 0 . Now we can introduce a frequency detuning parameter ω by ω ≡ ω − ω r0 = − (sin φ 0 + α cos φ 0 ) If we require the condition G − ω Q < 0 in (A.13), we note that the phase shift φ 0 must satisfy − π 2 < φ 0 < π 2 . When the linewidth enhancement factor α is zero, φ 0 = 0 at ω = 0 and φ 0 becomes ± π 2 at the edge of the locking bandwidth ω L defined by The above locking bandwidth is obtained by requiring |sin φ 0 | 1 in (A.16) together with α = 0. When the linewidth enhancement factor α is finite, the locking bandwidth becomes asymmetric (A. 18) In this case, φ 0 = 0 is satisfied at a finite frequency detuning ω = −α(F 0 /A 0 ) √ (ω/Q e ) from (A. 16). Figure A.1 shows the phase φ 0 versus normalized detuning parameter ω/ ω L for various values of α.
The steady-state solutions for the electron number can be obtained from the non-fluctuating part of (6): where N C0 is the electron number in the steady-state oscillation above threshold. The threshold pump rate p th of a free-running laser is defined by .20) and the normalized pump rate R is The steady-state solutions for the output field are obtained from the real and imaginary parts of the non-fluctuating component in (A.12) When φ 0 = 0, r 0 = (ω/Q e )A 0 − F 0 and ψ 0 = φ 0 = 0. However, when φ 0 = 0, the phase of the output field is not identical to that of the internal field but given by (A.23).
In our analysis, we assume that all slave lasers have zero linewidth enhancement factor α = 0 and there is no detuning ( ω = 0) for simplicity. In this case, the phases of the internal field and the output field are completely locked to the phase of the injection signal from the master laser. When α = 0, it is better to lock the slave laser at the edge of the negative detuning ω − √ 1 + α 2 ω L , where stable locking is possible and both the amplitude and phase noise are small [22]. It is also assumed that there is no internal photon loss ((ω/Q 0 ) = 0) and the population inversion is ideal (n sp = 1). We also assume that the injection signal F 0 is much smaller than the internal field so that the locking bandwidth is many orders of magnitude smaller than the cavity decay rate ω Q .