Networks of non-equilibrium condensates for global optimization

Recently several gain-dissipative platforms based on the networks of optical parametric oscillators, lasers, and various non-equilibrium Bose-Einstein condensates have been proposed and realised as analogue Hamiltonian simulators for solving large-scale hard optimisation problems. However, in these realisations the parameters of the problem depend on the node occupancies that are not known a priori, which limits the applicability of the gain-dissipative simulators to the classes of problems easily solvable by classical computations. We show how to overcome this difficulty and formulate the principles of operation of such simulators for solving the NP-hard large-scale optimisation problems such as constant modulus continuous quadratic optimisation and quadratic binary optimisation for any general matrix. To solve such problems any gain-dissipative simulator has to implement a feedback mechanism for the dynamical adjustment of the gain and coupling strengths.

In the last five years we have seen the rapid emergence of a new field at the intersection of laser and condensed matter physics, engineering and complexity theories which aims to develop quantum devices to simulate classical spin problems faster than on classical von Neumann architecture using a gain-dissipative principle of their operation. If the simulated spin problem is NP-hard, then solving it efficiently opens a new route for solving many practically relevant NP problems [1]. Among such strongly NP-hard problems are the Ising and XY models [2], for which even an approximate solution is hard to find [3]. Several platforms were proposed and demonstrated the proof-of-principle for finding the global minimum of such spin Hamiltonians: injection-locked laser systems [4], the networks of optical parametric oscillators, [5][6][7][8], coupled lasers [9], polariton condensates [10], and photon condensates [11]. The main principle of such simulators is based on a gain process that raises the system above the threshold for a phase transition to a coherent state (a Bose-Einstein condensate, laser coherence, synchronisation of oscillators, etc.). Since the threshold is first reached at the state that maximises the occupation for a given pumping intensity, it is natural to expect that this state is related to a ground state of a particular spin Hamiltonian.
Such Hamiltonians are written for N classical "spins" s j = (cos θ j , sin θ j ). In the XY model "spins" are continuous, whereas for the Ising model the "spins" take discrete values as the phases θ j are restricted to 0 or π. The "spins" are coupled to each other with the "strengths" J ij that can be positive (ferromagnetic) or negative (antiferromagnetic), so that at the ground state "spins" arrange their orientation as to minimise Finding the global minimum of the Ising or XY Hamiltonians are also known as quadratic binary optimisation (QUBO) and the constant modulus quadratic continuous optimisation (QCO) problems, respectively.
In the gain-dissipative simulators the "spin" -the node (the "bit") of the simulators -is represented by the phase of a condensate at a particular spatial position [10,11] or by the phase of a coherent state generated in a laser cavity [8,9], so by the phase of the so-called coherent center (CC). All of the proposals to use the gain-dissipative simulators to find the absolute minimum of the spin Hamiltonians suffer from a serious limitation. As we show below, the coupling strengths J ij in such systems are modified by the occupations (number densities) of CCs i and j. The densities, however, are not a priory known and, for a general matrix J J J will be different from one CC to another. The density inhomogeneity may also lead to the system fragmentation when different subsystems acquire their own mode of oscillation, and the coherent steady state across all CCs will never be reached. In all the previous experimental realisations of the gain-dissipative simulators, an explicit or implicit assumption was made that the coupling terms are small so that each laser or condensate is stabilised independently at the same steady-state amplitude [4,5,[8][9][10]. Such assumption is justified only for the simplest structures of the coupling matrix J J J, where all CCs are almost equally connected with about the same coupling strengths. The found solution for a more general matrix is bound to be either only approximate or invalid. The problem of unequal densities has been recognised before [12], but the proposed method to reduce the heterogeneity in densities was to drive the system using randomly modulated signals and can yield only very modest improvement and only if the densities are quite close for the unmodulated signals.
In this paper we formulate the technological requirements for gain-dissipative platforms to be used as analog Hamiltonian optimisers. We develop a general framework for the operation of the gain-dissipative analogue simulators based on the Langevin gain-dissipative equations written for a set of CCs. We derive the rate equations for the geometrically coupled CCs such as in polariton or photon condensates. We show that by establishing a feedback connection between the gain mechanism and the density of the CC we can drive the system to the coherent ground state of the XY model, while the minimisers of this ground state will give the true minimum precisely for the externally provided coupling strengths. This framework allows us to formulate the hardware requirements for a physical realisation of any such simulator to achieve such a minimum and argue that such requirements are within the recent technological capabilities.
The operation of a gain-dissipative simulator consists of two stages: bosonic stimulation below the threshold and the coherence of operations at and above the threshold. As one ramps up the power of the gain mechanism (e.g. laser intensity) the gain overcomes the linear losses and is stabilised by the nonlinear gain saturation. The emergent coherent state maximises the total number of particles (minimises losses) and, therefore, minimises the Ising or XY Hamiltonian depending on whether the phases of the wavefunctions describing the CCs are constrained to 0 and π or not, respectively. To derive the governing equations for the system one can describe each CC at a position r = r i by a classical complex function Ψ i (t). Depending on the system, the couplings K ij between CCs can have a different origin: they can be geometrically induced by the particle outflow from other CCs [10] or induced by the mutual injection rate between lasers [8,13] or spatially separated condensates; they can be controlled by external potentials [11]. In what follows we derive the rate equations for the geometrically induced flow for the network of non-equilibrium condensates (such as e.g. polariton or photon condensates) starting from the mean-field description by the Ginzburg-Landau equations [14]. Such derivation has been done before using the infinite quantum well orthogonal basis [15] which is not natural for the localized condensates in the lattice. As the result, the phases of the time-dependent modes are not directly relevant to the phases of the condensates. Instead, we will use the natural basis associated with the wavefunction of an individual condensate. The Ginzburg-Landau equation that gives a phenomenological description of the wavefunction of a gain-dissipative condensate [14] is where ψ(r, t) is the wavefunction of the system,Ũ is the strength of the interaction (pseudo) potential, the sum on the right-hand side represents the rate of adding particles in N spatial is a given spatially localised pumping profile, that creates the condensate with a wavefunction φ i (r) = φ(|r − r i |) centred at r = r i with a normalised number density so that 2π |φ(r)| 2 r dr = 1. Furthermore, f i is the timedependent part of the pumping at r = r i , γ c is the rate of linear losses andσ is the rate of the density-dependent losses. In writing Eq. (1) we let = 1, m = 1/2. If the distances between CCs are larger than the width of P (r), the wavefunction of the condensate can be is the time-dependent complex amplitude. We substitute the expression for ψ into Eq. (1), multiply by φ * j for j = 1, ..., N , and integrate in space to yield the rate equations for a j (t) where a a a = {a i }, An asymptotics of φ(r) for the Gaussian pumping profile have been developed in [16]. Function φ can be approximated by φ(r) = β 2 exp[−βr + ik c r], where k c is the outflow velocity and β is the inverse characteristic width of the condensate [16]. The integrals χ ij = φ i φ * j dr for i = j can be evaluated in elliptical coordinates in terms of the Bessel functions [17] where l ij = |r i −r j |. We assumed that the CCs are well separated: l ij β 1, so that for i = j we have χ ij 1 as follows from Eq. (3). The integrals in the off-diagonal terms in Q Q Q are of the same order in smallness as χ ij . Since at the threshold |a i (t)| 2 are small, the off-diagonal terms in Q Q Q are quadratic in small quantities and can be neglected to the linear order in Eq.
(2), so that Q Q Q is the diagonal matrix with elements |a i | 2 a i q, where q = 2π |φ| 4 r dr.
To the linear order terms in small quantities |a i (t)| 2 , p ij and χ ij for i = j Eqs. (2) become where we used the quadratic smallness of ) which represents intrinsic vacuum fluctuations and classical noise with a diffusion coefficient D which disappears at the threshold. In other platforms such as the Coherent Ising Machines [5] the injection does not have to be symmetric between the nodes -this will be modelled by introducing an asymmetry parameter δ so that ∆ ij = γ inj i (t) + (1 − δ)γ inj j (t), where δ = 0 for symmetrically coupled CCs. Equations (4) are the rate equations on the CCs coupled with the strengths ∆ ij K ij . Note, that by writing the coupling strength in this form we separated the effect of what is not a priory known (pumping intensity, energy at the threshold etc.) from K ij that depends on the characteristics of the system that are known and for geometrically coupled condensate, for instance, depend on the distance between CCs i and j. To show how Eqs. (4) lead to the XY model minimization we rewrite them in terms of the number densities ρ i and phases θ i using the Madelung transformation Ψ i = √ ρ i exp[iθ i ] while suppressing the noise: where θ ij = θ i − θ j . The first term on the right-hand side of Eq. (6) tends to provide θ i with its own frequency of oscillations whereas the second term tends to synchronise the phases.
Phase synchronisation in such a system has been extensively studied especially in the context of semiconductor laser arrays [18,19]. To guarantee that the steady state of the system is reached and coincides with a minimum of the XY model, however, one needs to make sure that the gain mechanism is chosen so that all densities ρ i are the same. Only under this condition Eq. (6) describes the gradient decent to the minimum of the XY model. Previously, any realisation of the XY model using the laser systems (or non-equilibrium condensates) was based on the assumption that all lasers (condensates) have the same steady-state photon (particle) number [8,18,19]. This limits the problems that can be addressed by such a framework to trivial ones where all CCs have an almost equal number of connections with almost the same pumping rate. To implement any couplings and connectivities, one needs to be able to control the pumping rate of the individual CCs and bring all the densities to the same value at the threshold. We schematically illustrate the operational principle of such control mechanism in Fig. 1a. Starting from below the threshold, all CCs are equally pumped as shown at some t = t 1 . Depending on the node connectivity, the non-zero densities emerge at different rates for each CC as the pumping intensity increases and takes some CCs above the specified threshold ρ = ρ th as illustrated on Fig. 1a at some later time t = t 2 .
The pumping mechanism must be adjusted for each CC to enable the saturation at the same density: decreased for CCs with densities above the threshold and increased for the CCs below the threshold. The feedback mechanism can be implemented via optical delay lines (in a NOPO system), by adjusting the injection via the spatial light modulator (in the polariton and photon condensates) or by electrical injection (e.g. in the polariton lattices [25]). The mathematical description of such a feedback mechanism is where is a parameter that can be tuned to control the speed of approaching the threshold.
The fixed point of Eqs. (5-7) is ρ i = ρ th = (γ inj i − γ c + j;j =i ∆ inj ij K ij cos θ ij )/σ, with the total particle number given by M = ( i γ inj i − N γ c + i,j;j =i ∆ inj ij K ij cos θ ij )/σ. The condensation will first take place at the minimum of i γ inj i , therefore, at the minimum of the XY Hamiltonian H XY = − i,j;j =i ∆ inj ij K ij cos θ ij . The minimum is approached from below by gradually raising the pumping above the threshold which facilitates the achievement of the global minimum as Fig. 1a at t = t 3 illustrates.
In the context of the cGLE with a saturable nonlinearity [22] appropriate for the description of non-equilibrium condensates the nonlinear dissipation σ can be proportional to the pumping intensity σ = cγ inj i , where c is a system dependent parameter that for polariton condensates for instance depends on the decay rate of the particles in the hot exciton reservoir.

The steady state of Eqs. (5-7) becomes
In this case, at the threshold i 1/γ inj i is maximized, so again the system reaches the minimum of the XY Hamiltonian.
If one removes the density heterogeneity the global minimum of the XY model will be achieved, but the coupling terms ∆ inj ij (t)K ij now depend on the particle injection rates γ inj i that would not be known a priori if one requested the equal densities at the threshold.
Therefore, not only γ inj i (t) has to be adjusted in time to equalise the densities using Eq. (7) but also the coupling parameters K ij have to be modified in time to bring the required couplings J ij at the steady state by whereˆ controls the rate of the coupling strengths adjustment. Sinceˆ such adjustments do not significantly slow down the operation of the simulator as they have to be performed much more rarely then adjustments of the gain. Equation (8) indicates that the couplings need to be reconfigured depending on the injection rate: if the coupling strength scaled by the gain at time t is lower (higher) then the objective coupling J ij , it has to be increased To illustrate the implementation of the density and coupling adjustments and, therefore, realisation of the global minimum of the XY model by the actual physical system, we apply the developed algorithm to the square lattice of polariton condensates experimentally achieved in our previous work [10]. Figure 2a shows the density profile of 45 polariton condensates that interact by the outflow of the particles from neighbouring CCs. All CCs are equally pumped (with, say, γ inj ≡ γ inj i ) and, therefore, the CCs away from the margins have the largest occupation -they are fed by the particles coming from the eight neighbours.
The CCs at the margins have the lowest occupation as they interact with only four or five neighbours. Such density heterogeneity between the lattice sites is clearly seen on Fig. 2a with the condensates on the margins being barely visible. The resulting configuration realises the global minimum of the XY model, but with the coupling strengths between i−th and j−th condensates given by 2γ inj K ij √ ρ i ρ j with number densities ρ i and ρ j that are not known before the system reaches the configuration shown in Fig. 2a. The numerical simulation of the 7 × 7 polariton lattice shows not only the density variation between the sites in agreement with the experimental result but indicates the formation of a spin wave state ( Fig. 2b) which manifests the presence of various couplings in the lattice. To realise the global minimum of the XY model for the given couplings (J ij ) we need to implement the feedback mechanisms described above and that we illustrate step by step. First, we remove the density heterogeneity by adjusting the gain mechanism described by Eq. (7). The resulting pumping profile is shown in Fig. 2c with the corresponding steady state number densities and phases in Fig. 2d. The spin wave in the presence of equal densities between the lattice sites is due to the different pumping intensities, and therefore, different couplings ∆ ij K ij between the CCs across the lattice. We adjust K ij according to Eq. (8) by changing the distances between the sites as Fig. 2e illustrates. The final steady state has equal densities and equal antiferromagnetic coupling strength between the nearest neighbours with phases alternating between 0 and π to give the expected global minimum of the XY model.
The developed procedure for the dynamical adjustment of the gain and coupling strengths, which is reflected in Eqs. (4,7,8) could be simulated by a classical computer and lead to a new class of global optimisation algorithms as we explore elsewhere [23] in particular, analysing if or when the actual physical gain-dissipative simulator can outperform the classical computer algorithm due to inherent parallelism of spanning various phases before the configuration that enables coherence at the threshold is reached and whether quantum superpositions contribute to processing of the phase configurations. These results can also be extended to other spin Hamiltonians such as the Ising and Potts [24].
In this paper we discussed gain-dissipative systems such as OPOs, optical cavities, lasers, and non-equilibrium condensates. The physics of the gain-dissipative oscillators may cover a much large variety of systems such as electronic circuits, voltage controlled oscillators, microelectromechanical systems, spin-torque nano-oscillators, oxide-based oscillators, etc.
Our analysis is applicable to a broad family of oscillatory networks regardless of the nature of each oscillator as long as it is governed by the generic Eqs. (4), can be operated stably at the threshold and allows for the feedback mechanism of adjusting the pumping intensity of +10% +10% +12% FIG. 2: (a) Contour plot of the number density of the polariton condensates formed by nonresonantly pumping polaritons with equal intensities (reproduced with permission from [10]). The distance between the neighbours is such that the coupling strength is antiferromagnetic for a polariton dyad with such separation at the condensation threshold. Red lines show the particle fluxes between the sites: the central site experiences the inflow of the particles from eight neighbouring sites whereas the sites on margins have only four or five neighbours. Dashed figures embrace the condensates with densities lower than the central condensates. (b,d,f) Contour plots of the steady state number density function ρ = |ψ| 2 dr obtained by the numerical integration of the full dynamical governing equations for 7 × 7 lattice and for the parameters used previously [10]. The lattice constant is 7.5µm. (c,e) Contour plots of the pumping profiles at the steady state. Panels (c,d) were obtained by applying the density adjustments according to Eq. (7). Panels (e,f) were obtained by applying both the density adjustments and coupling adjustments according to Eq. (7) and Eq. (8). The resulting pumping intensities at the lattice sites are indicated for the top right quarter of (c,e) as the factor of the pumping at the lattice centre. The coupling strength adjustments are achieved by shifting the lattice sites as shown in red for the bottom left corner only in (e).