Gain-dissipative systems of various physical origin have recently shown the ability to act as analogue minimisers of hard combinatorial optimisation problems. Whether or not these proposals will lead to any advantage in performance over the classical computations depends on the ability to establish controllable couplings for sufficiently dense short- and long-range interactions between the spins. Here, we propose a polaritonic XY-Ising machine based on a network of geometrically isolated polariton condensates capable of minimising discrete and continuous spin Hamiltonians. We elucidate the performance of the proposed computing platform for two types of couplings: relative and absolute. The interactions between the network nodes might be controlled by redirecting the emission between the condensates or by sending the phase information between nodes using resonant excitation. We discuss the conditions under which the proposed machine leads to a pure polariton simulator with pre-programmed couplings or results in a hybrid classical polariton simulator. We argue that the proposed architecture for the remote coupling control offers an improvement over geometrically coupled condensates in both accuracy and stability as well as increases versatility, range, and connectivity of spin Hamiltonians that can be simulated with polariton networks.
Various physical systems have recently emerged as unconventional computing platforms with a potential to successfully compete against the established state-of-the-art classical techniques in solving large-scale combinatorial optimisation problems. Among the newly proposed computational machines are the coherent Ising machine (CIM) based on the degenerate Optical Parametric Oscillators (OPOs) , , , a CMOS-based Fujitsu digital annealer , an all-electronic oscillator-based Ising machine , a simulated bifurcation Toshiba machine , a spatial light modulator (SLM) based photonic Ising machine , an optical simulator with the injection-locked multicore fibre lasers , and a molecular computing machine on a programmable microdroplet array . All these approaches aim to achieve a much faster, more efficient, and more accurate way of solving a particular class of optimisation problems, namely the quadratic unconstrained binary optimisation (QUBO) problem. In statistical physics QUBO problems appear when one seeks the global minimum (the ground state) of the Ising Hamiltonian. Indeed, the explicit polynomial mappings into the Ising Hamiltonian of many practically important problems of the discrete optimisation are available , such as the travelling salesman, graph colouring, warehouse inventory management, and low-risk portfolio optimisation. All these problems belong to the NP-hard computational complexity class meaning an exponential asymptotic growth of the number of operations with number of variables. It is no wonder that the advent of unconventional ways of finding the ground state of the Ising Hamiltonian was accompanied by the development of new classical algorithms as well as new methods to compare the performance of different platforms. The former resulted in newly proposed physically-inspired algorithms ,  while the latter facilitated the design of instances of interaction matrices with the controlled complexity and planted solutions , . In addition to the Ising model, other spin Hamiltonians have been proposed for solving with physical platforms including the XY spin Hamiltonian simulators based on the photon condensates , mode-locked fibre lasers , coupled laser , and nanolaser arrays .
A system of polariton condensates has attracted a considerable interest over the last few years by offering an alternative gain-dissipative system for tackling both discrete and continuous optimisation problems. Polaritons are the bosonic quasi-particles arising from strong coupling between excitons and photons in a semiconductor microcavity . Due to extremely low effective mass, polaritons can undergo Bose–Einstein condensation at temperatures higher than atomic condensates leading to two-dimensional networks of condensates operating at ambient conditions. Macroscopic coherence of such networks is characterised by a complex classical field with a well-defined condensates’ relative phases . These phases can be mapped into continuous ‘spins’ that can be further constrained to discrete values by employing the resonant excitation . The idea underlying the polariton simulator for solving optimisation problems originates from the belief that a huge combinatorial space of possible states can be sought in parallel near the condensation threshold, at which only the low-energy coherent states can form. At a coherent state, the phase-locked condensates have the same frequency, chemical potential, and constant relative phases. This state may correspond to the local or global minimum of a particular spin Hamiltonian, and since condensation occurs on a picosecond time scale, such polariton simulators may be potentially attractive for optimisation tasks.
Arbitrary networks of polariton condensates can be experimentally created in many different ways including optical imprinting , , , , in etched micropillars , , , in lead halide perovskite lattices , , in strain-induced traps , in hybrid air gap microcavities , by periodically etching the sample surface depositing metallic patterns , , , by surface acoustic waves , by direct fabrication with the gold deposition technique , in coupled mesas etching during the growth of the microcavity . The first of these, optical imprinting, is commonly realised with a liquid crystal SLM. The robustness of this technique has been demonstrated for a variety of lattice configurations and sizes , , , ,  proving the scalability of the polariton system for both trapped geometries and freely expanding condensates. The coupling between geometrically coupled condensates is generally of a complex nature  and consists of the dissipative (Heisenberg) and Josephson couplings. The latter prevents the system from achieving the minimum of a spin Hamiltonian as was previously elucidated . Even when the Josephson coupling is negligible in comparison with the dissipative coupling, the geometric coupling barely allows one to control the interactions beyond the nearest neighbours. This prevents the system from addressing complex, non-planar spin Hamiltonians. Nevertheless, the mapping of combinatorial problems into polariton networks was originally suggested by controlling the distance between geometrically coupled condensates ,  and later followed by propositions to control interactions in regular lattices of condensates with dissipative gates  or resonant pump barriers . The initial experimental demonstration of minimising the XY spin Hamiltonian with simple polarion lattice cells  has been shortly followed by a theoretical proposal extending the class of simulated spin Hamitonians from continuous to discrete .
Finding ways to dynamically control individual interactions between network nodes is a necessary step for addressing non-trivial spin Hamiltonians but not sufficient. In all proposed schemes the nearest neighbour interactions are attempted to be controlled while the beyond nearest neighbour interactions are assumed to be negligible, which is rarely the case. A recent study has shown the synchonisation between condensates across distances over 100 µm  noting that a typical lattice size constant is often in the range of 5–15 µm . This leads to a crucial and yet missing discussion of controlling the couplings beyond nearest neighbours for arbitrary graphs of polariton condensates. Moreover, spatially coupled polariton condensates are capable of representing different oscillator models  including the Kuramato, Sakaguchi–Kuramoto, Lang–Kobayashi, and Stuart–Landau models for different ranges of experimental parameters some of which are easier to adjust, e. g. exciton-polariton interactions, while others are harder, e. g. polariton lifetime. This apparent flexibility of a polariton system makes it harder to isolate a particular optimisation problem to address with polariton networks and, consequently, limits the optimisation accuracy of any objective function. To distinguish between many models that can be modelled with polariton networks, an instrumental calibration of experimental parameters may be required for spatially coupled polariton condensates even for nearest neighbour interactions.
In this work, we offer an alternative approach for simulating spin Hamiltonians with a network of spatially localised polariton condensates that do not interact with one another geometrically. For the network to become a spin Hamiltonian optimiser, we propose to couple any two condensates by redirecting the emission from one condensate to another or by exciting one condensate with an additional resonant pump tuned to the phase of that condensate. We emulate the polariton simulator with the two-dimensional complex Ginzburg–Landau equation (cGLE) coupled to the exciton reservoir equation and consider two possible coupling schemes, namely relative and absolute couplings. The performance of the emulated polariton simulator is demonstrated for discrete, i. e. Ising, and continuous, i. e. XY, spin Hamiltonians for sparse and dense interaction matrices J of various sizes from 9 to 49 condensates. This manuscript does not introduce a new optimisation algorithm and, therefore, the common metrics for comparing algorithms, e. g. time-to-solution, are intentionally omitted. Instead, we offer a proof-of-concept of a real XY-Ising polariton machine that can be realised within the vast range of experimental parameters. We outline the conditions under which the proposed polariton simulator becomes a “pure” or “hybrid-classical” physical optimiser.
where is the condensate wavefunction, m is the effective polariton mass, is the density distribution of the exciton reservoir, and are the polariton-polariton and polariton-exciton interaction strengths respectively, is the rate of scattering from the reservoir into the condensate, η is the energy relaxation. The term is an optional resonant pump at the double condensate frequency (second resonance) which forces phase differences between different condensates to be either 0 or π . Linear photon losses in the cavity and exciton losses are described by the condensate and the reservoir damping rates. The pumping intensity characterises the incoherent excitation. The process of Bose–Einstein condensation is essentially quantum, however, once a condensate is formed it can be accurately described by the cGLE as was shown in numerous experimental works , , , , , , , , . The optimisation accuracy of either XY or Ising models does not change for small values of η that are asumed in experiments  and, therefore, we will neglect it for the rest of the manuscript with a note that the bigger values would have had a negative effect on the accuracy. Introducing the dimensionless variables and parameters by , , , , , , , , , , l0 = 1 µm, we obtain
The cGLE (Eq. (3)) is a universal order parameter equation that describes non-linear phenomena in driven-dissipative systems ranging from non-linear waves and second-order phase transitions to lasers and superconductors. In this work, we use these equations to represent a network of isolated non-interacting polariton condensates which can be experimentally realised, for instance, with micropillars or with trapped polariton condensates. The former requires a lithographically modified sample and etching and leads to the formation of a polariton condensate, which coexists with the exciton reservoir density in each micropillar. The latter can be achieved without modifying the sample, e. g. by exciting each network’s element with a Gaussian ring pump, which would form a polariton condensate separated from the exciton reservoir. Although the following analysis can be readily applied to either experimental configuration, for ease of reading we will use an array of micropillars as our primary example of isolated condensates with occasional notes on the possible change in performance of the other. The position, shape, and size of each micropillar can be accurately controlled during their fabrication . Hundreds of coupled micropillars etched in a planar semiconductor microcavity have been used to study a wealth of phenomena from the Dirac cones in a honeycomb geometry  to the gap solitons in 1D Lieb lattices . To model the polariton condensation in a micropillar cavity, we introduce a spatially dependent dissipative profile,
where and are the dissipation rates outside and inside of a micropillar, respectively. Here, denotes the center of the i-th micropillar, and is the degree of a supergaussian that models micropillars as flat low-dissipative discs. The dramatically increased dissipation between the discs effectively blocks all the polariton outflows which leads to non-interacting condensates even for short separation distances of a few micrometres as would be expected for the system of micropillars. The condensates at different micropillars are noninteracting unless either relative or absolute remote couplings are introduced. In the former case, a part of the light emitted by the j-th micropillar condensate is re-injected into the i-th micropillar condensate at the amount proportional to the occupation of the j-th condensate. In the case of the absolute coupling, the same amount of light is exchanged between the i-th and the j-th condensates. Both coupling models can be represented by
where is the delta-function which is equal to one inside a micropillar, i. e. when , and zero outside, N is the number of micropillars, and τ represents a possible time-delay to supply couplings in an experimental setup. The coupling term represents the emission feedback when for each in a micropillar i the respective values are added from the micropillar centred at j. For the relative coupling model we shall consider while for the absolute coupling model we will use . The sign of the coupling strength can be made positive or negative by injecting the light with zero (ferromagnetic coupling) or π phase (anti-ferromagnetic coupling), respectively. For further derivations, we denote as the Heaviside function where R is the radius of the central part of the micropillar with a uniform phase distribution. In Eq. (6) we assume that the frequencies of each individual micropillar may be slightly different just below the condensation threshold. Nevertheless, the condensation process in presence of interpillar couplings locks these frequencies of different condensates resulting in a single energy condensate level . Recent experimental reports on two coupled micropillar lasers have demonstrated such frequency locking for detunings of up to 1 GHz in the few photons regime , . For negligible time-delay and geometric couplings between condensates, we can rewrite Eqs. (4) and (6) for each micropillar i using and as N equations for the polariton condensates and N equations for reservoir densities noting that , :
The steady states of Eqs. (7) and (8) correspond to the minima of the XY, i. e. , and Ising, i. e. , models as it becomes evident after we substitute into Eq. (8) and separate the real and imaginary parts. The equations read as
where . Here we considered the uniform phase distribution which is a valid assumption near the micropillar’s centre, i. e. with being the micropillar’s radius. In case of the relative coupling scheme, the fixed points of Eqs. (9)–(11) represent the minima of the XY or Ising spin Hamiltonians only for the equal polariton densities  across all micropillars, that is, when the condition stands. Such density equilibration can be robustly achieved by iteratively updating pumping intensity so that for all micropillars, where is the predefined integral luminosity. In contrast, the absolute coupling model naturally optimises the XY and Ising models and doesn’t require the equalised polariton densities due to the coupling coefficients that represent the exchange of a fixed number of photons between sites. The steady state solution of Eqs. (9)–(11) is given for both coupling models by equations:
where µ is the global oscillation frequency shared between all condensates at a coherent state.
One can see from the Eq. (14) that for a fixed point solution the maximised total polariton density corresponds to the minimum of the total reservoir density, which together with Eq. (12) leads to the minimisation of the spin Hamiltonians:
where denotes the plane of the microcavity. The resonant force term acts as a penalty in the objective function and leads to optimisation of the Ising model while the XY Hamiltonian is optimised for zero penalty term. We note that the term has a destabilising effect on the steady states solutions corresponding to minima of spin Hamiltonians meaning that small exciton-polariton interactions and/or small exciton reservoirs could possibly improve the optimisation accuracy. In experiments, a small reservoir density can be achieved for a high conversion rate of excitons into polaritons or by spatially separating polaritons from the reservoir by considering, for example, trapped condensates.
The validity of the proposed relative and absolute coupling models is verified by applying the two-dimensional Eqs. (7) and (8) for optimisation of the XY and Ising Hamiltonians on various coupling matrices. Firstly, we determine the minimum value of the coupling strength required for phase-locking of two condensates. Figure 1(top) shows the phase difference for a polariton dyad in the case of different interaction strengths with a zero time-delay. For each coupling strength, we simulate 50 random initial conditions and calculate the phase difference between the condensates in a final steady state. The region of decoupled condensates can be identified for coupling strengths by observing random phase differences between the condensates in Figure 1(top). For bigger coupling strengths, the condensates become phase locked and can reach ferromagnetic ground state (with zero phase difference between the condensates) for the positive couplings or antiferromagnetic ground state (with π phase difference) for the negative couplings. The local minima become unstable for coupling strengths bigger than 0.05 and the system finds the ground state regardless of the initial conditions. The demonstrated minimum coupling strengths for phase-locking of two condensates are similar for both relative and absolute coupling models in case of the XY Hamiltonian. For the Ising Hamiltonian, the destabilisation of excited states (local minima) happens for bigger coupling strengths of about . This is therefore the minimum coupling strength needed for the system to find the dyad’s ground state independently of the initial conditions. We note that the presence of intrinsic noise has a positive effect on destabilising such local minima.
In an experimental implementation of interactions, a possible time-delay τ may appear in constructing couplings between the network elements due to multiple reasons including the phase readout time, the time required to re-route photons, or the time for adjusting an SLM. As a result, the delayed phase information of condensates at time will be used for creating couplings between the condensates at time t whose phases will be shifted due to the global oscillation frequency. To demonstrate this effect of a time-delay in realising coupling strengths between different micropillars, we consider the absolute coupling model in optimising the XY Hamiltonian. Figure 1(bottom) shows the phase difference dependence on the time-delay for the polariton dyad with the coupling strength . The percentage time-delay is defined as a ratio to the time T that is required for the dyad to reach a steady state in the absence of the delay. For each time-delay value, we simulate 50 random initial conditions and show the resulting phase difference with scatter points of varied sizes proportional to the fraction of initial conditions that lead to this phase. The anticipated anti-ferromagnetic ground state is observed for time-delays τ up to . The previously unstable local minimum, i. e. ferromagnetic state with 0 phase difference for , becomes now stable in the presence of time-delay. Interestingly, the subsequent de-synchronisation area is followed by a clear ferromagnetic coupling between condensates which is in turn followed by another anti-ferromagnetic area for . This peculiar synchronization behaviour can be attributed to the global phase rotation with frequency µ of each condensate which can lead to phase-locking of condensates with an additional π phase difference for large time-delay values. This time-delay effect is similar for both coupling schemes in simulating either spin Hamiltonian. Although for networks of condensates, the presence of a time-delay would result in a phase lag  in Eqs. (14)–(12) which for significant τ can decrease the optimisation accuracy of the XY Hamiltonian, but not Ising. For simplicity, in the following investigations, we will not consider any time-delay in the couplings.
Having established the minimum coupling strength for phase-locking of two condensates, we now consider nine fully-connected polariton condensates. Each condensate is created with a non-resonant Gaussian pump in a lattice of three by three condensates (see Figure 2(a)). To realise spatially non-interacting polariton condensates we introduce a dissipative profile as shown in Figure 2(b) where the absence of particle outflows is ensured by the high value of outside nodes compared to low values inside nodes. A random interaction matrix is constructed of positive and negative couplings of amplitude as shown in Figure 2(c). As an illustrative example, we apply the relative and absolute coupling models described by Eqs. (7) and (8) for optimising the XY Hamiltonian (). In the former case, the densities of condensates are iteratively equalised over time by individually adjusting pumping intensities . The absolute coupling model does not require equal polariton densities at the steady state and, consequently, non-equal densities can be realised in a final state. The phase configurations and corresponding density profiles are shown in Figure 2(d–e) for the lowest energy states out of 10 runs for both models. To quantify the optimisation performance of coupling models, we consider the median accuracy that is defined by a proximity to the ground state:
where is the spin Hamiltonian energy for the phase configurations obtained with the mean-field approach (Eqs. (7) and (8)) in case of the relative or absolute coupling schemes, is the ground state energy found by the classical optimisation algorithms. In Figure 2(d–e), the found minima are within and from the ground state of the XY Hamiltonian that was verified with the gain-dissipative  and the basin-hopping  algorithms (Figure 2(f)). The median accuracy over 100 random fully-connected matrices of size generalises to and for the XY Hamiltonian in case of the relative and absolute coupling models, respectively.
To investigate the performance of the proposed polaritonic XY-Ising machine on the bigger size problems, an analysis of the optimal range of coupling values and edge density effects is required. In what follows we study the relative and absolute coupling models on the random unweighted MaxCut problems for the XY and Ising spin Hamiltonians. For the unweighted MaxCut problem, one seeks to divide the graph into two subgraphs with the maximised number of edges between them. This problem is known to be NP-hard  and can be mapped to the Ising Hamiltonian by assigning antiferromagnetic couplings to the graph edges. We construct three such random adjacency matrices of size of degree 5, 9, and 13. Both coupling models are simulated on matrices with amplitudes in the range . For each coupling strength amplitude, the Eqs. (7) and (8) are simulated for 20 random initial conditions. Figure 3(a–b) shows the ground state proximity as a function of amplitude for the XY Hamiltonian. The optimal range of couplings with the median accuracy over can be identified for the amplitudes in the range for the relative coupling model and slightly smaller range of for the absolute coupling model. For the Ising model, a smaller batch of coupling amplitudes allows one to achieve the median accuracy greater than (see Figure 3(c–d)). Such difference between the optimal coupling ranges can be possibly anticipated since the hard problems for the Ising Hamiltonian are not necessarily hard for the XY Hamiltonian optimisation. The clear shift to bigger optimal couplings for bigger edge densities (>0.8) is especially pronounced for the Ising Hamiltonian. This analysis confirms the lower bound and provides the upper bound of the coupling strength for achieving higher optimisation accuracies for both models. We note that the ground states of the Ising Hamiltonians were verified with the gain-dissipative  and CIM  algorithms.
With the identified optimal range of coupling amplitudes, we apply the relative and absolute coupling models to bigger spin Hamiltonian problems. Table 1 shows the median accuracy for both coupling models simulated on 20 unweighted MaxCut instances of size 25 and 49 with edge density of . For such connectivity, we pick the amplitude strength of from the optimal range. The number of initial conditions is fixed to 20 per each coupling matrix. We say that the coupling matrix J is globally optimised if the actual ground state is found at least once out of 20 runs for the Ising Hamiltonian. In case of the XY Hamiltonian, we require at least one phase configuration that is closer than to the ground state for claiming global optimisation. This number of globally optimised interaction matrices is indicated in parentheses in Table 1. The relative coupling model shows a consistently better performance on both the Ising and XY Hamiltonians than the absolute coupling model. The less accurate results for the Ising Hamiltonian, which are even more pronounced for the absolute coupling model, may be due to the greater hardness of generated interaction matrices for discrete optimisation than continuous. The drastic difference between coupling models could be possibly mitigated with a better choice of or may be a signal of a better local minima escape mechanism of the relative scheme. Nevertheless, the demonstration of the optimal performance of either of the proposed coupling methods is not the focus of this manuscript since both methods can be easily outperformed by standard heuristic algorithms. Instead, the achieved results clearly demonstrate a proof-of-principle for using polariton condensates, modelled with the mean-field approach equations (7) and (8), as the XY-Ising computing machine.
|25 (5 × 5 lattice)||99.3% (20)||87.8% (20)||96.8% (20)||72.9% (14)|
|49 (5 × 5 lattice)||98.2% (20)||81.7% (4)||93.3% (3)||52.3% (0)|
3.1 Experimental implementation
The spatially non-interacting condensates can be experimentally realised using lithographically etched micropillars or with trapped polariton condensates. The couplings are established remotely according to the elements of the coupling matrix . We envision two types of remote couplings. In the first scheme, the couplings are constructed by redirecting the emission of each condensate with either free-space optics or optical fibres to an SLM. At the SLM, the signal from each node is multiplexed and redirected to other nodes with the desired coupling strength allowing one, in principle, to create an all-to-all coupled network. Each matrix of couplings J can be programmed on the SLM in advance. We refer to this implementation as all-optical implementation. In the second approach, the frequency and phase of the condensate emission are read out and fed forward to an additional resonant excitation. Such resonant excitation will have to be iteratively updated based on the phase and energy of the emission until the polariton network synchronises. Consequently, the time-performance of the second scheme would be dependent on the operational frequency of the reading system and the SLM, which could be on the order of a few kHz .
The comparable or better time-performance can be possibly achieved with the digital micro mirror devices which have a similar millisecond operational time-scale or with electro-optical modulators which can operate at nanosecond scale. We will refer to this implementation as hybrid-classical implementation, since the condensate must first form to acquire a well-defined phase that is read out and passed to other nodes. Note that in both implementations we consider symmetric interactions, i. e. for any two condensates in a network, though directional interactions can be readily constructed, e. g. by using an optical isolator.
In addition to two possible experimental implementations of the remote coupling control, we propose two kinds of couplings: absolute and relative. The absolute coupling scheme implies the exchange of equal amounts of photons (equal signals’ intensities) between i-th and j-th nodes and guarantees the use of the correct coupling matrix for the spin Hamiltonian minimisation. In the relative coupling scheme, the condensates are coupled at the rate defined by relative intensities of emission and, therefore, a further density adjustment is required . This adjustment is crucial for the operation of nonequilibrium condensates, lasers or Degenerate Optical Parametric Oscillators (DOPOs) as the density heterogeneity changes the values of the coupling strengths . Since the equilibration of densities will be done at the operation frequency of the SLM, the relative coupling model shares the same limitations as the hybrid-classical implementation.
Thus, the absolute coupling scheme with the all-optical implementation may lead to a pure polaritonic XY-Ising machine for optimising spin Hamiltonians since it doesn’t require any external control: all couplings of a given spin Hamiltonian can be programmed on the SLM in advance. By approaching the condensation threshold from below, the polariton network will condense at one of the lowest energy states corresponding to a local or global minimum of the spin Hamiltonian. The term ”pure” indicates that the system can operate at its own physical time-scale, i. e. picosecond scale for the polariton condensation. Among other pure physical simulators are the time-delay CIM  and the recently proposed pure molecular simulator . The absolute coupling scheme with the hybrid-classical implementation as well as the relative coupling scheme with either of the proposed implementations would lead to the classical hybrid polariton simulators with an operational time limited by the frequency of the SLM. These approaches would be reminiscent of the CIM with a measurement feedback via FPGAs  or hybrid molecular simulator .
Similar to other analogue optimisers with optical feedback, one of the fundamental bottlenecks of the proposed polaritonic machine may be the appearance of phase errors when creating couplings between the condensates. These phase errors may be present in the phase pattern imprinted on the micropillar array due to aberrations and misalignment of the optics interfacing the SLM with the array. Such errors can be minimised by a direct imaging of phase pattern of a resonant laser going through the SLM and reflecting from the array of micropillars. The measured pattern would be matched with the desired one that includes the engineered hoppings via slight adjustments of the SLM pattern. Assuming the worst case scenario when the phase errors may disturb the final phase configuration of the low energy state, we emphasise a potentially non-trivial nature of the observed state compared to minima that can be found by standard classical techniques. Consequently, the classical optimisation algorithms can benefit from a warm-start, i. e. when the solution of an analogue machine is supplied as an initial condition for a classical algorithm. Such possible hybrid optimisation approach has been recently considered for the D-Wave machine . We also note that the implementation of interactions in polariton machine can be potentially much simpler to realise comparing to the feedback loops in OPOs or pure SLM implementations. For example, the couplings can be robustly controlled by adjusting the resonant pump intensities in the relative coupling scheme irrespective of the hopping pattern of the chosen spin Hamiltonian. Many existing approaches including the recent XY  and Ising  simulators suffer from couplings that are dependent on the final steady state which significantly limits the optimisation accuracy for non-trivial coupling matrices. For the polaritonic optimiser, this effect is considered and further mitigated by adjusting non-resonant pumping intensities so that the output power is uniform throughout the network.
The potentially advantageous performance of polaritonic machine in optimising spin Hamiltonians stems from the nature of polariton quasiparticles. Polariton condensates have a much stronger nonlinearity (coming from self-interactions between polaritons) than any of the purely photonic or laser-based optimisers. The stronger interactions should allow easier and faster exploration of phase configurations during the condensation process and narrower linewidth for the final measurement. In addition, the Bose-condensation process itself may facilitate the efficient low-energy sampling of spin Hamiltonians in a polariton simulator thanks to quantum effects present during the coherence formation.
3.2 Polaritonic XY-Ising machine
In this work, we introduce a new approach for simulating discrete and continuous spin Hamiltonians, e. g. Ising and XY, with polariton networks. We propose two experimental implementations for realising remote phase locking of any two condensates in a micropillar array or in a lattice of trapped condensates with a potential to have fully-connected coupling matrices. The first scheme could possibly result in a pure optical polariton simulator in which the interactions are organised by redirecting the leaking photons from one condensate to another, therefore forming photonic feedback mechanism. The second leads to a hybrid-classical polariton simulator in which the interactions are realised with additional resonant injections. Both methods can be a viable option for building a real polaritonic XY-Ising machine. We verify the performance of the proposed machine by simulating polariton networks with the mean-field approach for two types of couplings between condensates: relative and absolute. Both methods clearly demonstrate the ability to optimise spin Hamiltonians of various sizes, up to 49 condensates, and various connectivities, up to 24 connections per element. Moreover, the possibility to simulate spin Hamiltonians with beyond nearest neighbour couplings is proposed for the first time in polaritonic networks. Thus, the proposed polaritonic machine possesses such essential qualities of an analogue optimiser as robust programmability of interactions via SLMs, the ability to simulate sparse and possibly fully-connected matrices, and the implementation of arrays up to thousand condensates with existing experimental techniques which have great potential for further scale-up. The strong-coupling regime of polariton quasi-particles should be advantageous for the bottom-up optimisation approach and facilitate the achievement of low-energy states by a parallel-scanning through all phase configurations near the condensation threshold. The real physical machine would further benefit from low noise to signal ratio, ultra-fast operational time-scale, high energy-efficiency with a milliwatt excitation power per condensate, and potential room-temperature operation.
4 Materials and methods
The numerical evolution of Eqs. (7) and (8) is performed with the fourth-order Runge-Kutta time integration scheme and fourth order spatial finite difference scheme. The simulation parameters are , , , , with for all micropillars in case of the absolute coupling scheme and dynamically adjusted to bring all the condensates to in case of the relative coupling scheme for the XY Hamiltonian, , the distance between micropillars is , , , , , the micropillar diameter is about , . In addition, the following parameters are used to simulate the resonant pumping: , where , is the time required to achieve a steady state.
Funding source: Engineering and Physical Sciences Research Council
Funding source: Cambridge Trust
Funding source: Huawei
Conflict of interest statement: The authors declare that they have no competing interests.
KPK acknowledges the support from EPSRC and Cambridge Trust. NGB acknowledges the support from Huawei.
 Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, “Coherent Ising machine based on degenerate optical parametric oscillators,” Phys. Rev. A, vol. 88, no. 6, 2013 Dec 30, Art no. 063853, https://doi.org/10.1103/physreva.88.063853.Search in Google Scholar
 P. L. McMahon, A. Marandi, Y. Haribara, et al., “A fully programmable 100-spin coherent Ising machine with all-to-all connections,” Science, vol. 354, no. 6312, pp. 614–617, 2016, https://doi.org/10.1126/science.aah5178.Search in Google Scholar
 T. Inagaki, Y. Haribara, K. Igarashi, et al., “A coherent Ising machine for 2000-node optimization problems,” Science, vol. 354, no. 6312, pp. 603–606, 2016 Nov 4, https://doi.org/10.1126/science.aah4243.Search in Google Scholar
 M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. Katzgrabeer, “Physics-inspired optimization for quadratic unconstrained problems using a digital annealer,” Front. Phys., vol. 7, p. 48, 2019, https://doi.org/10.3389/fphy.2019.00048.Search in Google Scholar
 J. Chou, S. Bramhavar, S. Ghosh, and W. Herzog, “Analog coupled oscillator based weighted Ising machine,” Sci. Rep., vol. 9, no. 1, pp. 1–10, 2019 Oct 15, https://doi.org/10.1038/s41598-019-49699-5.Search in Google Scholar
 H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,” Sci. Adv., vol. 5, no. 4, 2019 Apr 1, Art no. eaav2372, https://doi.org/10.1126/sciadv.aav2372.Search in Google Scholar
 D. Pierangeli, G. Marcucci, and C. Conti, “Large-scale photonic Ising machine by spatial light modulation,” Phys. Rev. Lett., vol. 122, no. 21, 2019 May 31, Art no. 213902, https://doi.org/10.1103/physrevlett.122.213902.Search in Google Scholar
 M. Babaeian, D. T. Nguyen, V. Demir, et al., “A single shot coherent Ising machine based on a network of injection-locked multicore fiber lasers,” Nat. commun., vol. 10, no. 1, pp. 1–11, 2019 Aug 6, https://doi.org/10.1038/s41467-019-11548-4.Search in Google Scholar
 S. Y. Guo, P. Friederich, Y. Cao, et al. “A molecular computing approach to solving optimization problems via programmable microdroplet arrays,” ChemRxiv, 2018.10.26434/chemrxiv.10250897.v1Search in Google Scholar
 K. P. Kalinin and N. G. Berloff, “Global optimization of spin Hamiltonians with gain-dissipative systems,” Sci. Rep., vol. 8, no. 1, pp. 1–9, 2018 Dec 12, https://doi.org/10.1038/s41598-018-35416-1.Search in Google Scholar
 T. Leleu, Y. Yamamoto, P. L. McMahon, and K. Aihara, “Destabilization of local minima in analog spin systems by correction of amplitude heterogeneity,” Phys. Rev. Lett., vol. 122, no. 4, 2019 Feb 1, Art no. 040607, https://doi.org/10.1103/physrevlett.122.040607.Search in Google Scholar
 Y. Pang, C. Coffrin, A. Y. Lokhov, and M. Vuffray. The Potential of Quantum Annealing for Rapid Solution Structure Identification. arXiv preprint arXiv:1912.01759, 2019 Dec 4.Search in Google Scholar
 F. Hamze, J. Raymond, C. A. Pattison, K. Biswas, and H. G. Katzgraber. “Wishart planted ensemble: A tunably-rugged pairwise ising model with a first-order phase transition,” Physical Review E, vol. 101, no. 5, pp. 052102, 2020 Jun 1.Search in Google Scholar
 B. Kassenberg, M. Vretenar, S. Bissesar, and J. Klaers. Controllable Josephson Junction for Photon Bose–Einstein Condensates. arXiv preprint arXiv:2001.09828, 2020 Jan 27.Search in Google Scholar
 S. Tamate, Y. Yamamoto, A. Marandi, P. McMahon, and S. Utsunomiya. Simulating the Classical XY Model with a Laser Network. arXiv preprint arXiv:1608.00358, 2016 Aug 1.Search in Google Scholar
 V. Pal, S. Mahler, C. Tradonsky, A. A. Friesem, and N. Davidson. Rapid Fair Sampling of XY Spin Hamiltonian with a Laser Simulator. arXiv preprint arXiv:1912.10689, 2019 Dec 23.10.1103/PhysRevResearch.2.033008Search in Google Scholar
 M. Parto, W. Hayenga, A. Marandi, D. N. Christodoulides, and M. Khajavikhan. “Realizing spin-hamiltonians in nanoscale active photonic lattices,” Nat. Mater., pp. 1–7, 2020 Dec 27.10.1038/s41563-020-0635-6Search in Google Scholar PubMed
 J. Kasprzak, M. Richard, S. Kundermann, et al., “BoseEinstein condensation of exciton polaritons,” Nature, vol. 443, no. 7110, pp. 409–414, 2006 Sep, https://doi.org/10.1038/nature05131.Search in Google Scholar
 K. P. Kalinin and N. G. Berloff, “Simulating Ising and n-state planar Potts models and external fields with nonequilibrium condensates,” Phys. Rev. Lett., vol. 121, no. 23, 2018 Dec 4, Art no. 235302, https://doi.org/10.1103/physrevlett.121.235302.Search in Google Scholar
 E. Wertz, L. Ferrier, D. D. Solnyshkov, et al., “Spontaneous formation and optical manipulation of extended polariton condensates,” Nat. Phys., vol. 6, no. 11, pp. 860–864, 2010 Nov, https://doi.org/10.1038/nphys1750.Search in Google Scholar
 F. Manni, K. G. Lagoudakis, T. C. Liew, R. Andr, and B. Deveaud-Pldran, “Spontaneous pattern formation in a polariton condensate,” Phys. Rev. Lett., vol. 107, no. 10, 2011 Sep 2, Art no. 106401, https://doi.org/10.1103/physrevlett.107.106401.Search in Google Scholar
 G. Tosi, G. Christmann, N. G. Berloff, et al., “Sculpting oscillators with light within a nonlinear quantum fluid,” Nat. Phys., vol. 8, no. 3, pp. 190–194, 2012 Mar, https://doi.org/10.1038/nphys2182.Search in Google Scholar
 G. Tosi, G. Christmann, N. G. Berloff, et al., “Geometrically locked vortex lattices in semiconductor quantum fluids,” Nat. Commun., vol. 3, no. 1, pp. 1–5, 2012 Dec 4, https://doi.org/10.1038/ncomms2255.Search in Google Scholar
 D. Bajoni, P. Senellart, E. Wertz, et al., “Polariton laser using single micropillar GaAs GaAlAs semiconductor cavities,” Phys. Rev. Lett., vol. 100, no. 4, 2008 Jan 28, Art no. 047401, https://doi.org/10.1103/physrevlett.100.047401.Search in Google Scholar
 D. Sanvitto, F. M. Marchetti, M. H. Szymaska, et al., “Persistent currents and quantized vortices in a polariton superfluid,” Nat. Phys., vol. 6, no. 7, pp. 527–533, 2010 Jul, https://doi.org/10.1038/nphys1668.Search in Google Scholar
 T. Boulier, M. Bamba, A. Amo, et al., “Polariton-generated intensity squeezing in semiconductor micropillars,” Nat. Commun., vol. 5, no. 1, pp. 1–7, 2014 Feb 12, https://doi.org/10.1038/ncomms4260.Search in Google Scholar
 Y. Su, S. Z. Phua, Y. Li, et al., “Ultralong room temperature phosphorescence from amorphous organic materials toward confidential information encryption and decryption,” Sci. Adv., vol. 4, no. 5, 2018 May 1, Art no. eaas9732, https://doi.org/10.1126/sciadv.aas9732.Search in Google Scholar
 R. Su, S. Ghosh, J. Wang, et al., “Observation of exciton polariton condensation in a perovskite lattice at room temperature,” Nat. Phys., vol. 16, pp. 301–306, 2020 Jan 13.10.1038/s41567-019-0764-5Search in Google Scholar
 R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, “Bose–Einstein condensation of microcavity polaritons in a trap,” Science, vol. 316, no. 5827, pp. 1007–1010, 2007 May 18, https://doi.org/10.1126/science.1140990.Search in Google Scholar
 S. Dufferwiel, F. Fras, A. Trichet, et al., “Strong exciton-photon coupling in open semiconductor microcavities,” Appl. Phys. Lett., vol. 104, no. 19, 2014 May 12, Art no. 192107, https://doi.org/10.1063/1.4878504.Search in Google Scholar
 C. W. Lai, N. Y. Kim, S. Utsunomiya, et al., “Coherent zero-state and -state in an excitonpolariton condensate array,” Nature, vol. 450, no. 7169, pp. 529–532, 2007 Nov, https://doi.org/10.1038/nature06334.Search in Google Scholar
 A. Mischok, V. G. Lyssenko, R. Brückner, et al., “Zeroand πStates in a periodic array of deep photonic wires,” Adv. Opt. Mater., vol. 2, no. 8, pp. 746–750, 2014 Aug, https://doi.org/10.1002/adom.201400126.Search in Google Scholar
 L. Zhang, W. Xie, J. Wang, et al., “Weak lasing in one-dimensional polariton superlattices,” Proc. Natl. Acad. Sci., vol. 112, no. 13, pp. E1516–E1519, 2015 Mar 31, https://doi.org/10.1073/pnas.1502666112.Search in Google Scholar
 E. A. Cerda-Mndez, D. N. Krizhanovskii, K. Biermann, R. Hey, M. S. Skolnick, and P. V Santos, “Dynamic exciton-polariton macroscopic coherent phases in a tunable dot lattice,” Phys. Rev. B, vol. 86, no. 10, 2012 Sep 25, Art no. 100301, https://doi.org/10.1103/physrevb.86.100301.Search in Google Scholar
 Y. Kim, J. Ihm, E. Yoon, and G. D. Lee, “Dynamics and stability of divacancy defects in graphene,” Phys. Rev. B, vol. 84, no. 7, 2011 Aug 10, Art no. 075445, https://doi.org/10.1103/physrevb.84.075445.Search in Google Scholar
 K. Winkler, J. Fischer, A. Schade, et al., “A polariton condensate in a photonic crystal potential landscape,” New J. Phys., vol. 17, no. 2, 2015 Jan 30, Art no. 023001, https://doi.org/10.1088/1367-2630/17/2/023001.Search in Google Scholar
 H. Ohadi, A. J. Ramsay, H. Sigurdsson, et al., “Spin order and phase transitions in chains of polariton condensates,” Phys. Rev. Lett., vol. 119, no. 6, 2017 Aug 8, Art no. 067401, https://doi.org/10.1103/physrevlett.119.067401.Search in Google Scholar
 H. Ohadi, Y. D. Redondo, A. J Ramsay, et al., “Synchronization crossover of polariton condensates in weakly disordered lattices,” Phys. Rev. B, vol. 97, no. 19, 2018 May 8, Art no. 195109, https://doi.org/10.1103/physrevb.97.195109.Search in Google Scholar
 N. G. Berloff, M. Silva, K. Kalinin, et al., “Realizing the classical XY Hamiltonian in polariton simulators,” Nat. Mater., vol. 16, no. 11, 2017 Nov, pp. 1120–1126, https://doi.org/10.1038/nmat4971.Search in Google Scholar
 K. P. Kalinin and N. G. Berloff, “Polaritonic network as a paradigm for dynamics of coupled oscillators,” Phys. Rev. B, vol. 100, no. 24, 2019 Dec 26, Art no. 245306, https://doi.org/10.1103/physrevb.100.245306.Search in Google Scholar
 K. P. Kalinin, P. G. Lagoudakis, and N. G. Berloff, “Matter wave coupling of spatially separated and unequally pumped polariton condensates,” Phys. Rev. B, vol. 97, no. 9, 2018 Mar 23, Art no. 094512.10.1103/PhysRevB.97.094512Search in Google Scholar
 K. P. Kalinin and N. G. Berloff, “Toward arbitrary control of lattice interactions in nonequilibrium condensates,” Adv. Quant. Technol., vol. 3, no. 2, 2020 Jan 7, Art no. 1900065, https://doi.org/10.1002/qute.201900065.Search in Google Scholar
 S. Alyatkin, J. D. Topfer, A. Askitopoulos, H. Sigurdsson, and P. G. Lagoudakis. Optical Control of Synchronous Phases in a Programmable Polariton Cell. arXiv preprint arXiv:1907.08580, 2019 Jul 19.Search in Google Scholar
 J. D. Topfer, H. Sigurdsson, L. Pickup, and P. G. Lagoudakis, “Time-delay polaritonics,” Commun. Phys., vol. 3, no. 1, pp. 1–8, 2020 Jan 7, https://doi.org/10.1038/s42005-019-0271-0.Search in Google Scholar
 J. Keeling and N. G. Berloff, “Spontaneous rotating vortex lattices in a pumped decaying condensate,” Phys. Rev. Lett., vol. 100, no. 25, 2008 Jun 23, Art no. 250401, https://doi.org/10.1103/physrevlett.100.250401.Search in Google Scholar
 M. Wouters and I. Carusotto, “Excitations in a nonequilibrium Bose–Einstein condensate of exciton polaritons,” Phys. Rev. Lett., vol. 99, no. 14, 2007 Oct 3, Art no. 140402, https://doi.org/10.1103/physrevlett.99.140402.Search in Google Scholar
 P. Cristofolini, A. Dreismann, G. Christmann, et al., “Optical superfluid phase transitions and trapping of polariton condensates,” Phys. Rev. Lett., vol. 110, no. 18, 2013 May 1, Art no. 186403, https://doi.org/10.1103/physrevlett.110.186403.Search in Google Scholar
 H. Ohadi, A. Dreismann, Y. G. Rubo, et al., “Spontaneous spin bifurcations and ferromagnetic phase transitions in a spinor exciton-polariton condensate,” Phys. Rev. X, vol. 5, no. 3, 2015 Jul 8, Art no. 031002, https://doi.org/10.1103/physrevx.5.031002.Search in Google Scholar
 K. G. Lagoudakis, M. Wouters, M. Richard, et al., “Quantized vortices in an excitonpolariton condensate,” Nat. Phys., vol. 4, no. 9, pp. 706–710, 2008 Sep, https://doi.org/10.1038/nphys1051.Search in Google Scholar
 K. G. Lagoudakis, B. Pietka, M. Wouters, R. Andre, and B. Deveaud-Pledran, “Coherent oscillations in an exciton-polariton Josephson junction,” Phys. Rev. Lett., vol. 105, no. 12, 2010 Sep 15, Art no. 120403, https://doi.org/10.1103/physrevlett.105.120403.Search in Google Scholar
 T. Gao, E. Estrecho, K. Y. Bliokh, et al., “Observation of non-Hermitian degeneracies in a chaotic exciton-polariton billiard,” Nature, vol. 526, no. 7574, pp. 554–558, 2015 Oct, https://doi.org/10.1038/nature15522.Search in Google Scholar
 M. D. Fraser, G. Roumpos, and Y. Yamamoto, “Vortexantivortex pair dynamics in an excitonpolariton condensate,” New J. Phys., vol. 11, no. 11, 2009 Nov 25, Art no. 113048, https://doi.org/10.1088/1367-2630/11/11/113048.Search in Google Scholar
 T. Jacqmin, I. Carusotto, I. Sagnes, et al., “Direct observation of Dirac cones and a flatband in a honeycomb lattice for polaritons,” Phys. Rev. Lett., vol. 112, no. 11, 2014 Mar 18, Art no. 116402, https://doi.org/10.1103/physrevlett.112.116402.Search in Google Scholar
 V. Goblot, B. Rauer, F. Vicentini, et al., “Nonlinear polariton fluids in a flatband reveal discrete gap solitons,” Phys. Rev. Lett., vol. 123, no. 11, 2019 Sep 13, Art no. 113901, https://doi.org/10.1103/physrevlett.123.113901.Search in Google Scholar
 K. P. Kalinin and N. G. Berloff, “Networks of non-equilibrium condensates for global optimization,” New J. Phys., vol. 20, no. 11, 2018 Nov 14, Art no. 113023, https://doi.org/10.1088/1367-2630/aae8ae.Search in Google Scholar
 E. Schlottmann, S. Holzinger, B. Lingnau, et al., “Injection locking of quantum-dot microlasers operating in the few-photon regime,” Phys. Rev. Appl., vol. 6, no. 4, 2016 Oct 31, Art no. 044023, https://doi.org/10.1103/physrevapplied.6.044023.Search in Google Scholar
 S. Kreinberg, X. Porte, D. Schicke, et al., “Mutual coupling and synchronization of optically coupled quantum-dot micropillar lasers at ultra-low light levels,” Nat. Commun., vol. 10, no. 1, pp. 1–11, 2019 Apr 4, https://doi.org/10.1038/s41467-019-09559-2.Search in Google Scholar
 D. J. Wales and J. P. Doye, “Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms,” J. Phys. Chem. A, vol. 101, no. 28, pp. 5111–5116, 1997 Jul 10, https://doi.org/10.1021/jp970984n.Search in Google Scholar
 G. Lazarev, A. Hermerschmidt, S. Kruger, and S. Osten, “LCOS spatial light modulators: trends and applications,” Opt. Imag. Metrol. Adv. Technol., pp. 1–29, 2012 Aug 22, https://doi.org/10.1002/9783527648443.ch1.Search in Google Scholar
 A. Marandi, Z. Wang, K. Takata, R. L. Byer, and Y. Yamamoto, “Network of time-multiplexed optical parametric oscillators as a coherent Ising machine,” Nat. Photon., vol. 8, no. 12, p. 937, 2014 Dec, https://doi.org/10.1038/nphoton.2014.249.Search in Google Scholar
© 2020 Kirill P. Kalinin et al., published by De Gruyter, Berlin/Boston
This work is licensed under the Creative Commons Attribution 4.0 International License.