A Quantum Model for Coherent Ising Machines: Discrete-time Measurement Feedback Formulation

TL;DR

This paper introduces a quantum theoretical model for coherent Ising machines using discrete-time measurement feedback, providing insights into quantum coherence during optimization.

Contribution

It presents a novel discrete-time measurement feedback model for CIMs and performs numerical simulations on a simple two-oscillator network.

Findings

Quantum coherence persists during the optimization process.

The model effectively captures the dynamics of CIMs with feedback.

Numerical results demonstrate the feasibility of the approach.

Abstract

Recently, the coherent Ising machine (CIM) as a degenerate optical parametric oscillator (DOPO) network has been researched to solve Ising combinatorial optimization problems. We formulate a theoretical model for the CIM with discrete-time measurement feedback processes, and perform numerical simulations for the simplest network, composed of two degenerate optical parametric oscillator pulses with the anti-ferromagnetic mutual coupling. We evaluate the extent to which quantum coherence exists during the optimization process.

Tables

Table 1. Table 1: The parameters for numerical simulation for the coherent Ising machine with two DOPOs.

physical meaning	name	value
Net Gain in One loop	$G_{tot}$	1.05
Background Loss Rate	$1-T^{\prime }$	0
Feedback Rate	$R$	0.005
Two Photon Loss Rate	$L$	0.002

Table 2. Table 2: Parameters for the numerically simulation in Figure 6 .

physical meaning	name	value
Measurement Strength	$T$	0.99
Feedback Rate	R	0.005
Two Photon Loss Rate	L	0.002

Equations

$$\mathcal{H}$$

$$B=i\sum _{k}g({\omega }_k)e^{-i{\omega }_{k}t}b({\omega }_k),$$

$$b(2\omega )\to b(2\omega )+\beta ,$$

$$\mathcal{H}$$

$${\mathcal{H}}_{int}^I$$

$$\frac{d{\rho }^I}{dt}$$

$${\Gamma }_k$$

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

$$\frac{d}{dt}⟨x|\rho |x^{\prime }⟩$$

$$U^{†}aU=\sin \theta a_{vac}+\cos \theta a,$$

$$U^{†}{a}_{vac}U=\sin \theta a+\cos \theta a_{vac}.$$

$$M_{x_m}$$

Full text

I introduction

There is significant interest in finding alternatives to modern von-Neuman computers, in particular for solving combinatorial optimization problems in NP-hard and NP-complete classes. Artificial neural networks can implement NP-hard Ising problems Hopfield and Tank (1985) and NP-complete k-SAT problems Ercsey-Ravasz and Toroczkai (2011). Other efforts include adiabatic quantum computation Farhi et al. (2001), and quantum annealing Kadowaki and Nishimori (1998) with superconducting quantum circuits. However, providing dense connectivity between qubits in such physical systems remains a major challenge in achieving a satisfying efficiency with these alternative approaches Rieffel et al. (2015).

Recently, many types of Coherent Ising machines (CIM) have been studied to solve Ising-type combinatorial optimization problems Nigg et al. (2017); Puri et al. (2017); Goto et al. (2017). Among these systems, the CIMs using denegerate optical parametric oscillator (DOPO) networks are also being discussed Wang et al. (2013); Takata et al. (2015); Maruo et al. (2016). The first generation of this CIM implements spin-spin coupling through optical delay lines Marandi et al. (2014); Roslund et al. (2014); Takata et al. (2016); Inagaki et al. (2016a); Treps et al. (2016). This type of direct coupling Ising machine can implement $O(N^{2})$ spin-spin connections with $N-1$ optical delay lines, which remains still a technical challenge for a large spin size ($N\gg 1$). Recently, the DOPO networks with the measurement feedback circuit were implemented at Stanford University McMahon et al. (2016) and NTT Inagaki et al. (2016b). In these machines, the oscillators are coupled indirectly with the discrete-time quantum measurement and feedback processes, which can generate classical correlations between oscillators. In this paper, we formulate a quantum model with completely positive trace preserving (CPTP) maps of the measurement feedback process for the first time, and numerically evaluate the extent to which quantum coherence exists during the optimization process. This paper is organized as follows. In Section II, we introduce our theoretical model of the DOPO network with the discrete-time measurement feedback processes. In Section III, we present numerical simulation results for a simple DOPO network. Finally, in Section IV we conclude with a brief summary.

II The Theoretical Model

Our theoretical model consists of four components: a PPLN waveguide as a phase sensitive amplifier; two output couplers for simulating the measurement loss and background loss; and a feedback circuit consisting of optical homodyne detectors, an analog-to-digital converter (ADC), a field-programmable gate array (FPGA), a digital-to-analog converter (DAC), and an optical amplitude/phase modulator (Figure 1). The first output coupler represents all of the background loss in the ring cavity. The signal fields extracted by the second output coupler are used to measure the in-phase amplitudes from the homodyne detectors. In the feedback process, the feedback pulses are generated as coherent states with an average excitation amplitude $x_{i}=\sum_{j}J_{ij}\tilde{x_{j}}$, where $\tilde{x_{j}}$ is the measurement result for the j th pulse and $J_{ij}$ is the Ising coupling constant. Each signal pulse undergoes these four processes each time it completes one round trip along the ring cavity. To simulate the way in which the states of signal pulses evolve, we calculate the CPTP maps of measurement feedback processes. Then, we numerically simulate the system by expanding the field density operators in terms of the eigenvectors of the in-phase amplitude operator $x=(a+a^{\dagger})/\sqrt{2}$ 111When we define the in-phase and quadrature-phase amplitude operators as $x=(a+a^{\dagger})/\sqrt{2}$ and $p=(a-a^{\dagger})/\sqrt{2}$, they satisfy the commutation relation $[x,p]=i$ and the uncertainty principle $\langle\Delta x^{2}\rangle\langle\Delta p^{2}\rangle=1/4$ ., where $a$ and $a^{\dagger}$ are the annihilation/creation operators of the signal pulse. Because these processes along the ring cavity consists of local operators and classical communications, the states of the signal pulses are not entangled. This is in sharp contrast to the direct optical coupling DOPO network Takata et al. (2015); Maruo et al. (2016). We calculate conditional density matrices governed by the randomly determined measurement results $x_{m}$ of the in-phase amplitudes of the signal pulses extracted by the second output coupler.

In the following subsections, we derive the time evolution equation for the density operators and CPTP maps for the four processes in the ring cavity. For simplicity, we consider a rotating coordinate and ignore the free field Hamiltonian.

II.1 Phase sensitive amplifier (PSA)

In the PSA, the signal pulse with an initial state $\rho$, frequency $\omega$, and annihilation operator $a$, interacts with the pump pulse, whose initial state is a coherent state $|\beta\rangle$ ($\beta$ is real) with a frequency $2\omega$. The pump pulse after the interaction is thrown away and in the next round trip, new pump pulse $|\beta\rangle$ is generated.

The amplitude $\beta$ of the pump pulse is very large and the interaction with the signal pulse is weak.

The Hamiltonian for the PSA is

[TABLE]

Here, $B$ is defined as follows:

[TABLE]

where $b(\omega_{k})$ is the annihilation operator of the pump field of frequency $\omega_{k}$, and $g(\omega_{k})$ is a parametric coupling constant. Initially, the pump field at $2\omega$ is in a coherent state $|\beta\rangle$ and all the other modes are assumed to be in vacuum state $|0\rangle$. When we displace the pump field as follows

[TABLE]

the initial state of the pump field is displaced to a vacuum state. In this case, the Hamiltonian (1) is rewritten as

[TABLE]

We shall interpret that the first term $\mathcal{H}_{s}$ represents the linear phase sensitive amplification/ deamplification (squeezing effect) imposed by the pump, while the second term represents up conversion process, which results in an effective two photon loss of the signal mode. The first part of the Hamiltonian can be absorbed as an internal Hamiltonian of the signal. Then, the second part is expressed as

[TABLE]

in the interaction picture. We take the Born-Markov and the rotating wave approximations to eliminate the Hilbert space of the pump pulse. Then, we can find the master equation corresponding to the interaction Hamiltonian as follows Breuer and Petruccione (2002) :

[TABLE]

where

[TABLE]

When we go back to the Schrodinger picture, we have

[TABLE]

The first term in Eq.(7) represents a standard unitary (squeezing) process, while the second term is a Lindblad form representing the two photon loss process associated with the parametric pump photon generation. We defined the squeezing rate $S=g(2\omega)\beta t$ and the two photon loss rate $L=g(2\omega)^{2}t$, where $t$ is the time duration of this interaction in the PSA. The linear power-gain $G$ can be represented as $G=\exp(2S)=\exp(2g(2\omega)\beta t)$. One of the important assumption leading to Eq.(7) is that the gain saturation is relatively weak, i.e. the signal pulse intensity grows and depletes the pump power only slightly, instead of decaying to zero, as is the case for a traveling-wave PSA with strong signal-pump interaction Hamerly et al. (2016). The other important assumption is that the pump field is dissipated into external reservoirs each time of PSA and a fresh coherent state $|\beta\rangle$ is always prepared as a new pump field for the next round of PSA. The equation for our simulation is obtained by expanding the density operator $\rho$ in terms of the in-phase amplitude eigenstates $|x\rangle$ as follows:

[TABLE]

where $z=x+x^{\prime}$ and $w=x-x^{\prime}$.

All coefficients which appears in the equation above are real-valued because annihilation and creation operators, which appear in the equation (7), are summation of in-phase amplitude operator $x$ and its derivative $\partial_{x}$ with real coefficients, which allows $\langle x|\rho|x^{\prime}\rangle$ are real-valued in CIM.

II.2 Output couplers and homodyne detectors

A portion of the signal pulse is extracted from the ring cavity by the two beam splitters. At the first output coupler, the extracted signal-field is simply dissipated in external reservoirs, which represents the background loss in the ring cavity. At the second output coupler, the in-phase amplitude $x=(a+a^{\dagger})/\sqrt{2}$ of the extracted field is projectively measured by the homodyne detectors. We define the transmittance of the first and second splitters as $T^{\prime}=\sin^{2}\theta^{\prime}$ and $T=\sin^{2}\theta$. When the signal pulse goes into the beam splitter, it is combined with the incident vacuum state from the external environments. Thus, the measurement performed by the homodyne detectors has a finite measurement error, which stems from the vacuum fluctuation. We define the annihilation operators of the signal and the vaccum field as $a$ and $a_{vac}$. Then, the output field annihilation operators can be written in terms of the unitary operator $U$ of the beam splitter with a parameter $\theta$ as follows:

[TABLE]

[TABLE]

From these equations, the Kraus operator corresponding to the measured value of $x_{m}$ can be expressed by

[TABLE]

To calculate the conditional density matrix for the post-measurement state, we generate a random number and determine a measured value $x_{m}$ with the probability of $Tr(M_{x_{m}}\rho M_{x_{m}}^{\dagger})$. For the operation for the first output coupler, we use the same Kraus operators $\{M_{x_{m}}\}$ and ensemble many conditional density matrices governed by the probabilistically determined $x_{m}$.

II.3 Feedback process

In the feedback injection process, the signal pulse and feedback pulse, the latter of which is prepared in a coherent state $|\alpha\rangle$, are combined with a third beam splitter. The transmittance rate of the third beam splitter, defined as $T_{f}=\sin^{2}\theta_{f}$, is very high ($T_{f}\approx 1$). In this parameter region, the quantum flucutation due to the injected coherent state is much smaller than the flutucation in the signal pulse, and this feedback process can be described with a simple unitary displacement operator $D(\alpha\theta_{f})=\exp(\alpha\theta_{f}a^{\dagger}-\alpha^{*}\theta_{f}a)$. In the Heisenberg picture, the in-phase amplitude operator $x$ will be translated as $D(\alpha\theta_{f})xD^{\dagger}(\alpha\theta_{f})=x+\alpha\theta_{f}/\sqrt{2}$.

The amplitude $\alpha$ of the feedback pulse is determined by the measured values of the homodyne detectors. We define the feedback rate $R$ as the ratio of $\alpha\theta_{f}/\sqrt{2}$ to the in-phase amplitude $x$ of the signal pulse, estimated using the value $x_{m}$ measured by the homodyne detectors.

II.4 Summary of the modeling

In summary, the signal pulses experience four processes described above during each loop around the ring cavity. The system is described by five physical parameters, the gain $G$ (or squeezing rate $S$), two photon loss rate $L$, background loss rate $1-T^{\prime}$ , measurement loss rate $1-T$, and feedback rate $R$. The total net linear amplitude-gain, before the gain saturation is switchwed on, during one round trip of the cavity is $G_{tot}=\sqrt{G\times T\times T^{\prime}}=\exp(g(2\omega)\beta t)\sin\theta\sin\theta^{\prime}$, where we assume that $T_{f}\approx 1$.

For numerical simulations, we expand the conditional density matrices of pulses in terms of $x$-eigenvectors and calculate the elements of the density matrices. Note that the elements of density matrices $\langle x|\rho|x^{\prime}\rangle$ are real numbers in this system. In this paper, we will visualize density matrices as functions of $x+x^{\prime}$ and $x-x^{\prime}$. While this function is equivalent to the Wigner function because Wigner function can be obtained by the Frourier transformation of $\langle x|\rho|x^{\prime}\rangle$ along the $x-x^{\prime}$ axis, it is easier to see the quantum coherence than Wigner function as shown in Figure 2.

Figure 2 shows the contour maps of the functions $\langle x|\rho|x^{\prime}\rangle$ corresponding to typical quantum states. The line of $x-x^{\prime}=0$ is the diagonal line of a density matrix and represents the probability distribution on $x$, while the line of $x+x^{\prime}=0$ represents quantum coherence between the element of $|x\rangle$ and the one of $|-x\rangle$. The vacuum state $|0\rangle$ and a one photon state $|1\rangle$ can be represented with a simple gaussian function whose variance is $\langle\Delta x^{2}\rangle=0.5$ and a phase reversed Hermite Gaussian function (Figure 2(a),(b)). When it is anti-squeezed along the $x$-axis, it becomes broader gaussian distribution (Figure 2(c)). On the other hand, a thermal state is a gaussian whose width is large along the $x+x^{\prime}$ axis but small along the $x-x^{\prime}$ axis (Figure 2(d)). This clearly shows the low qunatum coherence in the thermal state. While non-zero values appear along the vertical line of $x+x^{\prime}=0$ when two coherent states $|\alpha\rangle$ and $|-\alpha\rangle$ are superposed with quantum coherence, they disappear when the two coherent states are classically ensembled (Figure 2(e)(f)). These states shown in Figure 2 and states in our CIM model don’t have imaginary elements, so the functions can be described with contour maps of real values.

III Numerical Simulation Results

We simulated the time evolution of the simplest DOPO network, consisting of two oscillators interacting with out-of-phase coupling (anti-ferromagnetic coupling). We consider that a [math] phase pulse is an up spin, and that a $\pi$ phase pulse is a down spin. In this case, the two degenerate ground states are [math] phase-$\pi$ phase (up-down $|\uparrow\downarrow\rangle$) and $\pi$ phase-[math] phase (down-up $|\downarrow\uparrow\rangle$) states.

III.1 Time evolution of typical conditional density matrices

In this subsection, we assume that the background loss in the cavity is zero ($T^{\prime}=1$), and the transmittance of the output coupler is $T=0.99$. Thus, in this case, the ratio of power extracted from the cavity for the measurement is 0.01. Other numerical parameters are given in Table 1. Figure 3 illustrates the time-development of typical conditional density matrices of signal pulses governed by the sequence of measured values $x_{m}$. The initial states of the two pulses are vacuum states ( with the number of round trips $N=0$ in Figure 3 (b)). At that time, $\langle x\rangle=0$ and $\langle\Delta x^{2}\rangle=0.5$. The optimization process consists of three stages. In the first stage, the in-phase amplitudes of two pulses are anti-squeezed by the phase sensitive amplifier, and $\langle\Delta x^{2}\rangle$ become larger ($N=30$ in Figure 3 (b) ). Note that both the diagonal x-distribution along the horizontal axis and the off-diagonal quantum coherence plotted along the vertical axis become broad. In the second stage, because of the gain saturation and feedback processes, the expectation value $\langle x\rangle$ moves to either negative or positive value ($N=60$ in Figure 3 (b)). The gain saturation and the linear photon loss are responsible for the spontaneous symmetry breaking of DOPO, while the feedback process makes the system to select an anti-ferromagnetic order instead of a ferromagnetic order. Finally, in the third stage the state becomes close to the highly excited coherent state and $\langle\Delta x^{2}\rangle$ is reduced to $0.5$ ($N=150$ in Figure 3 (b) ). At this stage, the DOPO state is already in a classical level and the optimization process of the CIM is completed.

III.2 Various measurement strengths

Next, we present the simulation results for various measurement strengths 1-$T$ under a condition of no background loss, $T^{\prime}=1$. The parameters for the numerical simulations are shown in Table 1. Figure 4 shows typical conditional density matrices $<x|\rho|x^{\prime}>$ governed by the sequence of measurement results $x_{m}$. It can be seen that the anti-squeezing effect at the early period of the optimization process is more significant when the transmittance of the output coupler is larger or the measurement strength is weaker so that the wavepacket reduction is not so significant. When $T=0.9$ or $T=0.5$, the states undergo mild anti-squeezing and they are quickly displaced. On the other hand, when $T=0.999$, the states maintains the quantum coherence between the macroscopically separated ”up state” and ”down state” ($N=60$ in Figure 4 (a)). Here, the probability distributions are not localized but the centers of the two wavepackets are negatively correlated. As shown in Figure 4(e), the Wigner function features an oscillatory behavior with negative amplitudes that manifests the quantum interference effect between the macroscopically separated ”up-state” and ”down-state”.

Figure 5 illustrates the relationship between $\langle\Delta x^{2}\rangle$ and $\langle\Delta p^{2}\rangle$ for a typical conditional density matrix and the full density matrix of a one signal pulse for various values of transmittance rate $T$. Here, $x,p$ are the in-phase and quadrature-phase amplitudes, respectively, defined by $x=\frac{a+a^{\dagger}}{\sqrt{2}},ip=\frac{a-a^{\dagger}}{\sqrt{2}}$. Initially at N=0, the state satisfies $\langle\Delta x^{2}\rangle=\langle\Delta p^{2}\rangle=0.5$. As the signal pulses complete many round trips around the ring cavity, the value of $\langle\Delta x^{2}\rangle$ for a conditional density matrix first becomes larger, and then decreases to $\langle\Delta x^{2}\rangle=0.5$ at well above DOPO threshold. These lines form loops. The dashed curve represents the minimum uncertainty product. The Heisenberg uncertainty principle dictates that $\langle\Delta x^{2}\rangle\langle\Delta p^{2}\rangle\geq 1/4$. Because the state is squeezed vacuum state in the early stage of the optimization process, it lies on this dashed curve when $\langle\Delta x^{2}\rangle$ is relatively small. The results depicted in Figure 5 clearly demonstrate the two facts. First, when the measurement strength is weaker, the states are more anti-squeezed and more quantum coherence is present between the ”up state” (the region $x>0$) and the ”down state” (the region $x<0$). Second, as $T$ becomes smaller or the measurement strength increases, the conditional density matrix becomes close to a squeezed coherent state, because the density matrix moves almost onto the curve of the minimum uncertainty product. In this case, the density matrices will be efficiently simulatable with the displaced squeezing basis method discussed in Tezak et al. (2017).

The degree of squeezing is more than 3dB in Figure 5 . It is well-known that a continuous wave (CW) pumped DOPO features only 3dB squeezing inside a cavity. However, such a limit does not exist for a pulsed DOPO. The numerical results shown in Figure 5(a)(b) are the uncertainty product at the input to the PSA, while that at the output of the PSA is shown in Figure 5(c).

III.3 Probability of success and the effect of background loss

After the many round trips of signal pulses, we projectively measure the in-phase amplitude of the pulses. If the amplitude is positive $x>0$, we treat it as a up-spin state, and otherwise, a down-spin state. Thus, in the case of two signals with out-of-phase coulping, the optimization by CIM is successful if the projectively measured in-phase amplitude of a signal is positive and the one of the other signal is negative. Thus, the definition of the success probability is $\int_{x_{1}x_{2}<0}\langle x_{1}|\langle x_{2}|\rho|x_{1}\rangle|x_{2}\rangle dx_{1}dx_{2}$, where $x_{1}$ and $x_{2}$ denote the in-phase amplitudes of the two signal pulses and $\rho$ is the full density matrix for the two signal DOPO pulses.

We produced many conditional density matrices, in order to calculate how the probability of success $P$ of the optimization depends on the background loss rate $1-T^{\prime}$, with the numerical parameters presented in Table 2 and three different pump schedules of net linear amplitude-gain $G_{tot}$. The results are shown in Figure 6.

The initial success rate for the two vacuum states is $0.5$ and as $N$ increases it becomes higher. To see the dependance of the probability of success on the background loss rate and the time schedule of the net linear gain, we set the feedback rate $R$ not to be sufficiently strong and $P$ does not reach 100%. Of course, when $R$ is sufficiently strong, $P$ reaches to 100%. In the case of a low background loss ($T^{\prime}=1.0$ and $0.9$), the initial increasing rate of $P$ is smaller than the case of a larger background loss. When the loss rate is low, the state is largely anti-squeezed, and this makes the signal-to-noise ratio of measurement is poor at the early stage, as shown in Figure 4. Thus, the in-phase amplitudes of the feedback pulses generated from the measured values suffer from small signal-to-noise ratio, and this leads to a lower increase rate for the probability of success at the early stage.

However, once $P$ starts to increase, it suddenly goes up and reaches to the final constant value. The final value greatly depends on the time shedule of the net gain $G_{tot}$. When $G_{tot}$ increases rapidly, the final probability of success for the low-loss case becomes lower than that in the case of slowly increasing $G_{tot}$. This is because the small fluctuation leads to the states of the DOPO pulses being easily trapped into the potential wells of DOPO (up-spin state or down-spin state) before the correlation between the two pulses permanently forms. On the other hand, when the background loss is large or $T^{\prime}$ is small ($T^{\prime}=0.5$ and $0.7$), the state fluctuates strongly and can move from the up (down) state to the down (up) state even after $G_{tot}$ becomes larger than one (above the threshold)Kinsler and Drummond (1991) . Thus, the success rate continues to increase also in the later stage of the optimization process.

IV Conclusions

We developed the quantum theory for the coherent Ising machines with the discrete-time measurement feedback process. In this optical system, signal pulses in a ring cavity are weakly measured and the mutual coupling via feedback pulses allows them to interact with each other. The CIM can detect the ground state of the Ising Hamiltonian with a certain probability. In this paper, we showed the simulation results of the simplest DOPO network, which is composed of two anti-ferromagnetically coupled signal pulses.

We showed that there are three stages in the optimization process: anti-squeezing of an initial vacuum state by phase-sensitive amplification, spontaneous symmetry breaking by gain saturation and linear photon loss with dispacement by the feedback process, and excess noise supression to approach the final coherent states.

We discussed quantum coherence during the optimization process. When the loss rate is small, the signal pulses undergo strong anti-squeezing along the in-phase amplitude, and quantum coherence between different amplitudes is maintained. On the other hand, in the case of a large loss rate, the state is close to the squeezed coherent states, and close to the Heisenberg limit $\langle\Delta x^{2}\rangle\langle\Delta p^{2}\rangle\approx 1/4$ during the optimization process.

We also calculated the success probabilities for different values of background loss rates $1-T^{\prime}$. In the case of a small background loss, the higher success probability is obtained with a slower increase in the net gain. On the other hand, in the case of a large background loss, the success probability is not as sensitive to the net gain increase schedule.

Acknowledgement

This research was funded by the Impulsing Paradigm Change through Disruptive Technologies (ImPACT) Program of the Council of Science, Technology and Innovation (Cabinet Office, Government of Japan). The authors wish to thank Taime Shoji, Ryan Harmerly, Peter Drummond and Hideo Mabuchi for their critical discussions.