Quantum Correlations in the Kerr Ising Model

In this article we present a full description of the quantum Kerr Ising model---a linear optical network of parametrically pumped Kerr non-linearities. We consider the non-dissapative Kerr Ising model and, using variational techniques, show that the energy spectrum is primarily determined by the adjacency matrix in the Ising model and exhibits highly non-classical cat like eigenstates. We then introduce dissipation to give a quantum mechanical treatment of the measurement process based on homodyne detection via the conditional stochastic Schrodinger equation. Finally we identify a quantum advantage in comparison to the classical analogue for the example of two anti-ferromagnetic cavities.


INTRODUCTION
In an increasingly connected and dynamic society, the demand to solve complex problems requiring optimal configurations of large systems is growing. Many of these optimisation problems fall into the NP-hard or NPcomplete complexity class that are practically impossible to solve on a classical digital computer. These types of problems can be mapped onto the Ising model where the ground state yields the optimal solution. Attempts to find the ground state has lead to the development of many classical and quantum approaches including adiabatic quantum computing [1] and quantum annealing [2,3]. However, a significant hurdle facing these architectures is the connectivity of individual physical qubits, an essential requirement for solving large optimisation problems.
A new approach dubbed the Coherent Ising Machine (CIM) overcomes this issue by implementing an optical analog of the Ising model [4][5][6][7]. Recent results have compared the performance of the CIM against semi-definite programs and simulated annealing [8] as well as neural network architectures [9] and quantum annealers [10]. Further theoretical results have shown comparable results can be obtained when modelled using Gaussian optics [11]. In the CIM, a series of degenerate parametric oscillators (DOPO's) form a coupled network of potential spins. At sufficiently high pumping, the DOPO experiences a pitchfork bifurcation into two steady state solutions. These two solutions are coherent states which are π out of phase with one another and play the role of a 'spin' in the Ising model. When coupled, each DOPO in the network bifurcates in accordance with an Ising Hamiltonian. The lowest energy configuration of coherent states in the cavities therefore corresponds to the ground state. The network of DOPO's become non separable during bifurcation [5,12] which may lead to a quantum advantage [13], although it is not yet clearly identified. Recently, several impressive classical experimental realisations of the CIM have been published demonstrating up to 2000 connected spins [6,10,13]. Furthermore, off-the-shelf electronics have been used to build robust opto-electronic CIMs [14].
Here we will consider a Kerr Ising Model (KIM) based on a network of parametrically pumped cavities containing a Kerr non-linearity. This interaction has been studied extensively with seminal results showing the onset of chaos [15], single and multi-photon blockade [16][17][18][19], as well as qubit construction [20], quantum gates [21][22][23] and computation [24,25]. Furthermore, experimental realisations of parametrically driven of Kerr non-linearities and their subsequent quantum behaviour have been discussed [26] and realised in superconducting circuits [27]. There is now a small but growing literature on driven networks of Kerr non-linearities in cascaded systems [28], under adiabatic evolution [24,29,30] and dissipative single photon loss [31,32]. Notably, the results presented in [29,30] describe the implementation of an Ising machine constructed from a network of Kerr non-linearities but do not consider the nature or effects of any quantum correlations; an analysis we present here.
We start by presenting a comprehensive theoretical description of the KIM and show that the network exhibits non-classical correlations at the corresponding classical bifurcation -indicating a quantum phase transition [33]. Beyond the bifurcation the ground state corresponds to an entangled state of minimum energy spin configurations of the Ising model. We then introduce low temperature dissipation using the quantum optics master equation and show that the steady state corresponds to a probabilistic mixture of these highly entangled ground states. This enables us to describe the conditional dynamics corresponding to continuous homodyne detection of the output fields showing the stochastic approach to the desired ground state. Lastly, we compare the fully quantum model against a classical analogue at finite temperatures for two coupled cavities and clearly identify the advantage arising from quantum correlations.

THE KERR ISING MODEL
In our analysis, we will start by assuming that the pump mode has been adiabatically eliminated. In the interaction picture, each cavity in the network is described by the Hamiltonian where χ is proportional to the third-order nonlinear susceptibility, ∆ = ω a − ω p /2 is the detuning of the pump frequency, ω p , from cavity resonance frequency, ω a , and is pump field. Due to the parity symmetry of Eq.(1), on resonance the ground state is a doubly degenerate subspace with energy E 0 = 0 spanned by two coherent states |α i = | ± α 0 where α 0 = /χ. The introduction of a small detuning ∆ creates a linear shift in the energy levels E 0 = ∆|α 0 | 2 .
Likewise, for a network of N uncoupled Kerr-cavities the ground state is 2N -fold degenerate spanned by product state | α m = | ± α 0 , ±α 0 , ... with a perturbed energy of If we now introduce a linear coherent coupling scheme in the form of a network of connected beam-splitters and phase shifters -we obtain the Hamiltonian [34] where η is fixed coupling strength for each cavity, S is the adjacency matrix between cavities and a = (a 1 , a 2 , ..., a N ) is a vector representation of annihilation operators for each cavity. On closer inspection, H bears considerable resemblance to the classical Ising model of N coupled spins σ i = ±1 in a uniform magnetic field B where H B = −µB i σ i is the magnetic dipole Hamiltonian, J is the uniform coupling between spins andŜ is again, the adjacency matrix between spins. In the KIM, the parametric driving of each Kerr non-linearity creates a subspace of pseudo-spins coupled through a linear optical network. This model mimics the magnetic field driving and the magnetic coupling between spins in the traditional Ising model. As a result, the energy spectrum in the KIM will reflect similar characteristics to the Ising model in Eq.(3). Using degenerate perturbation theory for a weak coupling η and large driving > 0, the nth shifted energy level of the network is determined by where α n = | |/χ σ n corresponds to the configuration vector of cavity coherent states for the nth eigenstate | α n [29]. The eigenstates of the new spectrum can be found by solving the perturbed eigenvector equation where | α = 2N i c ni | α i is written as a linear combination of the uncoupled coherent eigenstates | α i . Multiplying through by α k | and using the annihilation operator eigenvalue equation a| α = α| α we obtain the matrix equation As c ik = 0, this relationship is satisfied so long as α † k .Ŝ. α k = α † n .Ŝ. α n . This ensures any linear superposition of spin configurations α n which have the same energy will form a degenerate subspace. Given the parity symmetry in the Hamiltonian, the subspace will be spanned by the orthogonal cat-like superpositions of the Ising solutions.
To illustrate this, we consider the simplest 2-spin Ising model for zero detuning ∆ = 0 shown in Fig.(1). In the absence of coupling η = 0, the ground state is 4fold degenerate spanned by every unique configuration of spins |α i . Now introducing a weak coupling, the degeneracy is broken into a 2-fold degenerate ground state spanned by the anti-ferromagnetic states |λ − ∝ |α 0 , −α 0 ± | − α 0 , α 0 and a 2-fold degenerate excited state spanned by the ferromagnetic states |λ + ∝ |α 0 , α 0 ± | − α 0 , −α 0 . The new ground state in the high driving limit has an energy E g = −2 η| |/χ whereas the excited state has energy E e = 2 η| |/χ. A notable feature is the limit of low driving plotted in Fig.(1b) where the system is not degenerate. In this limit, the coupling η perturbation is strong compared to the driving and breaks the degeneracy.
The evolution of the KIM can be described by the Von-Neumann time development equation This may be converted to a Fokker-Planck equation for the positive P-representation of ρ using an expansion of the density operator over the off-diagonal projectors | µ ν|, where the diffusion term is Φ( µ, ν) = −i(χ µ 2 + ) and the drift is Ω( µ, ν) = −i(2χ µ 2 ν − 2 ν + ∆ µ + ηŜ. µ).
Replacing µ = α and ν = α * and ignoring the quantum diffusion terms yields the semi-classical equations of motion whereσ x,y,z are the Pauli spin matrices. The structure ofÂ is analogous to the spin model coupled alongσ z and driven alongσ x where the Hilbert space of the cavities constructs the pseudo-spin structure and interactions. The bistability of the system can then be computed when a sign change occurs in the trace tr(Â) or the determinant det(Â) with the variation of the driving [35].

ENTANGLEMENT IN THE KIM Ground state entanglement
The KIM has an energy spectrum analogous to the Ising model of a spin network immersed in a constant magnetic field. Earlier work on the quantum Ising model has shown the ground state undergoing a quantum phase transition, through which the degenerate ground states become highly entangled [33]. Likewise, previous dissipative CIM models have also shown that the entanglement and quantum discord occur in the pseudo-spin network [5,12].
There are many different tests and criterion for quantum correlations and entanglement. One such measure of non-separability of Gaussian states was introduced by Duan et al [36] and used in the dissiptive CIM models [5,12]. This criterion is a suitable choice in dissipative models if the steady state solution is Gaussian. The eigenstates of the KIM are highly non-Gaussian cat like states seen in the joint photon number distribution of n 1 ,n 2 |λ 0 plotted in Fig.(2)a which exhibits significant interference suggesting the presence of strong nonclassicality. Given the non-Gaussian nature of the eigenstates, using the Duan inequality is not a necessary condition of inseparability [36]. A more robust metric to compare between Hamiltonian and dissapative models is the logarithmic negativity (LN) defined as where ρ T A 1 = 2N + 1 is the trace norm defined in terms of the negativity N which is the absolute sum of all negative eigenvalues of ρ T A [37]. In the limit of maximally entangled pure states, the LN reaches equality with the Von-Neumann entropy metric of entanglement.
Plotting the LN in Fig.(2)b for a large non-zero detuning ∆ = 1 we see that the non-resonant pumping of the cavity pushes the bifurcation away from = 0 and plays an analogous role to the damping considered in the dissapative models-which we will consider in the next section-but retains the unitarity of the system. The ground state smoothly transitions to a highly entangled state and remains so for increased driving and coupling η.

Steady state entanglement
In an experimental setting, one needs to measure the values of each pseudo-spin in the model and thus the closed Hamiltonian system must be made open in some way. This could be an impulsive coupling to the external apparatus after some period of unitary evolution. In Ref. [29] the solution to the model was obtained under adiabatic evolution followed by a final read out. In quantum optics however a more conventional approach would be to open the system by letting a small amount of light leak out of the cavity and subject that to continuous measurement. Once the cavity is opened in this way the system becomes a damped non linear system and the elliptic fixed points of the Hamiltonian now become stable fixed points of a dissipative dynamical system. As for what to measure, we note that the key feature of the KIM solutions we seek are in the phase of the intracavity field so we will need to consider a phase-dependent measurement scheme such as homodyne detection. Given a measurement record we can now ask for the conditional dynamics of the intracavity field given a given stochastic homodyne current record.
At a non-zero temperature-where the heat bath is treated as white-noise on the input field-and for unit quantum efficiency homodyne detection, the conditional state obeys the conditional Stochastic Schrodinger equation (SSE) [38] , where all photons emitted from the cavity at rate γ are detected and dW i (t) is the classical Weiner process. The mean-photon number due to thermal excitation is determined by the Bose-Einstein statistics n = (exp [ ω/k b T ] − 1) −1 where T is the temperature and ω is the cavity resonance frequency (assumed to be identical). Here, D[a i ]ρ is the single photon loss channel and H[a]ρ = aρ + ρa † i − tr a i ρ + ρa † i ρ is the conditional stochastic term of the ith cavity. If we only seek the unconditional dissipative evolution, we can average over the stochastic term in Eq.(12) which then vanishes giving a master equation for the state of the damped cavity field. This master equation may now be converted to an equivalent Fokker-Planck equation as Eq. (8) with the addition of a decay term introduced into the drift Ω( µ, ν) → Ω( µ, ν) − γ µ/2 and an additional noise term γn ∂ 2 ν µ P ( µ, ν). Instead of asking about the structure of the ground state our attention now turns to the steady state solution to the master equation or equivalently to the Fokker Planck equation.
The LN of the steady state solution of two cavities as a function of pumping and coupling η is depicted in Fig. (2)c at zero temperaturen = 0. The entanglement is reduced due to the introduction of cavity dissipation but is non zero at the bifurcation. The two cavities become non-separable around parameter values corresponding to the classical bifurcation in the semi-classical dissipative model indicating a dissipative quantum phase transition [33]. The entanglement then decreases to zero in the limit of large driving. A similar peak in non-separability was also observed in the CIM [5,12].
In the case of the conditional dynamics based on homodyne detection, the cavities stochastically evolve into one particular spin configuration | ± α 0 , ∓α 0 . One can show the steady state solutions correspond to the different possible spin configuration sin the Ising model determined by the adjacency matrixŜ. This is achieved by finding the roots of Eq.9 when ∂ α/∂t = 0 [7,32]. As we can regard the unconditional steady state as arising from ensemble averaging over the measurement records and thus the conditional states, we conclude in the large driving limit, the steady state can be regarded as a mixture of the possible spin-state configurations that satisfy Eq.(6), slightly shifted due to the damping.

COMPARISON BETWEEN QUANTUM AND CLASSICAL MODELS
In the previous section we outlined the quantum dynamics of the 2 spin anti-ferromagnetic KIM subject to weak damping so as to allow constant homodyne detection through the conditional SSE Eq. (12). By modelling the conditional SSE we provide a numerical experiment of results that would likely be observed in a physical realisation of the KIM. We seek to understand more carefully what advantage may arise at the quantum limit. This can be achieved by making a comparison between the semi-classical model and the full quantum model at finite temperature as the temperature is lowered. Here we will compare the conditional field amplitude averages computed by solving the conditional SSE model with stochastic solutions to the SDE described by Eq. (9) with the addition of thermal noise where ξ 1 (t) and ξ 2 (t) are delta correlated stochastic forces.
The measured homodyne currents from either cavity, in the quantum case is J i (t) = X i (t) c + γ(2n + 1)ξ i (t) where the first term is a quantum average that is computed from the conditional state. In the classical case the analogous quantity is simply j i (t) = α i (t) + α * i (t). There are several quantities determined directly from the homodyne current capable of informing us about the performance of the KIM. Firstly, consider the difference squared of the measured currents of either cavity; ignoring the white noise terms this gives, (J 1 (t) − J 2 (t)) 2 = X 1 −X 2 2 c . (Note the square of a conditional average). In the anti-ferromagnetic example, if the two cavities converge on the same signed steady state then the mean difference will be zero. Conversely, if the two cavities settle into the anti-ferroamgnetic solution then this quantity will reach 2|α| 2 since the cavities bifurcate alongX. Notably, this quantity will also give us a measure of the intracavity dynamics.
Secondly the correlations between the individual homodyne currents J 1 (t) and J 2 (t) is fundamental to understanding the consequences of the quantum correlations ascertained in Fig.(2)c. Here we measure the normalised cross-correlation function between the two cavities defined where M is the number of sampled currents and σ 2 Ji is the variance in the current of the ith cavity. The correlation function is bounded −1 ≤ R X1,X2 ≤ 1 where equality signifies perfect (anti-)correlation.
Finally, the 'pseudo-spin' of the KIM is determined by the sign of the current J i (t). An error arises when both cavities have the same sign sign[J 1 (t)] = sign[J 2 (t)] in the anti-ferromagnetic example. Thus we can calculate the probability of obtaining an error as a function of time and measure the success probability of the KIM.
Here we sample both the quantum Eq. (12) and the classical Eq. (13) trajectories 2000 times-each equivalent to an experimental trial. From this sample we compute the average of the quantities depicted in Fig. (3). Both simulations have identical cavity parameters of = 1, ∆ = 0, χ = 0.5, γ = 1, η = 1 and are varied over several temperatures from 0 ≤n ≤ 1.
The results exhibit several interesting features. Firstly, increasing the temperature leads to expected agreement between the classical and quantum models as they approach the steady state. The system experiences bifurcation approximately where the cross-correlation peaks R X1,X2 . Before bifurcation, the two cavities exhibit positive correlations present in both classical and quantum models. This leads to an increase in detecting an error but rapidly decays once the system has bifurcated. In the limit of zero temperature (blue) the classical model stalls in the absence of thermal fluctuations: the KIM remains stuck at the unstable fixed point centred at the origin. The quantum mechanical model on the other hand still continues to function due to the spontaneous emission of the cavities. In-fact at zero temperature, it is less-likely to lead to an error in the steady state P r(error). Also, at zero temperature the quantum model exhibits extremely large positive correlations prior to bifurcation as a consequence of the squeezing. This quantum effect is rapidly suppressed in the presence of small amounts of thermal noise.
Another noticeable feature is the time taken to bifurcate. At all temperatures, the quantum mechanical model appears to peak in correlations at t ∼ 0.5. The mean difference in X 1 −X 2 2 appears to plateau between 1 < t < 2 before increasing again for larger times. The classical model on the other hand at low temperaturen = 0.25-is noticeably slower to begin convergence on the Ising solution but begins to resemble the quantum model for t > 2. In the classical model Fig. (3), this suggests a critical slowing in the convergence rate to a stable fixed point at lower temperatures. This tension introduces a trade off between speed to undergo bifurcation and accuracy; lower thermal fluctuations lead to higher accuracy but also result in critical slowing of the dynamics. At low temperatures, this critical slowing could be overcome by thermal annealing. In the quantum mechanical model at times below t < 2, the phase transition occurs sooner. This is likely due to the quantum entanglement between the cavities also known as quantum activation [39]. For larger times-after bifurcation-the cavities transition over into the classical regime as noted by the agreement between the results at finite tempera- In the classical model at zero temperaturen = 0, the system remains fixed at the unstable fixed point at the origin. The middle figure is a measurement of the normalised cross correlation function RX 1 ,X 2 . At zero the quantum trajectories are strongly positively correlated likely due to the non-classical squeezing in the cavities. The peak in positive correlations occurs when the system bifurcates and is present in both the classical and quantum models at finite temperature. The bottom figure shows the probability of measuring an error Pr(error) from the homodyne current. The positive correlations lead to a peak in errors but are then rapidly suppressed as the cavities settle into opposing 'pseudo-spin' states. The zero temperature limit again exhibits large probability of detecting an error prior to bifurcation due to the non-classical correlations.
This indicates one significant advantage over traditional quantum annealing architectures.
At zero temperature-a reasonable assumption in optical systems and superconducting circuits-with an initial state in the vacuum, the system is situated at an unstable stationary fixed point. Classically, a mechanism to induce this symmetry breaking is thermal annealing. Quantum mechanically, the cavities can quantum tunnel and begin converging on the stationary solutions determined by the Ising model [7]. Therefore no annealing is required in the KIM when finding the lowest energy solution.

DISCUSSION
In this article we have focused on the two-cavity KIMboth classically and the full quantum description. We have shown that the quantum model exhibits a highly entangled ground state but constantly measuring the cavities under homodyne detection forces the system converge on one particular spin-configuration solution. Consequently this leads to a classically correlated steadystate solution whereby quantum correlations exist only during the process of bifurcation. Regardless, the na-ture of this quantum phase transition elicits a quantum advantage in the trade-off between speed and accuracy of the model as we have shown comparing the classical model to the full quantum description.
We have shown conclusively that there exists an advantage to using a quantum mechanical KIM at the zero temperature limit as shown in Fig. (3). The quantum fluctuations due to spontaneous emission allow the KIM to operate in the limit where there are no thermal fluctuations. Furthermore, the non-classical correlations arising due to the quantum mechanical squeezing of the joint cavity state lead to an apparent decrease in the time taken for the system to reach bifurcation, likely resulting from quantum activation [39]. However it must be noted that this advantage rapidly decays with increasing temperature. Furthermore, we have considered the simplest two cavity model; further investigation is required to determine whether the advantage still exists at a larger number of cavities. It would also be interesting to consider how this advantage manifests itself with more complex adjacency matricesŜ and determine whether or not it is capable of outperforming other classical heuristic Ising solvers.