# Polaritonic XY-Ising machine

Kirill Kalinin, Alberto Amo, Jacqueline Bloch, Natalia Berloff

► To cite this version:

Kirill Kalinin, Alberto Amo, Jacqueline Bloch, Natalia Berloff. Polaritonic XY-Ising machine. Nanophotonics, 2020, 9 (13), pp.4127-4138. 10.1515/nanoph-2020-0162. hal-03037847

HAL Id: hal-03037847

https://hal.science/hal-03037847v1

Submitted on 16 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research article

Kirill P. Kalinin, Alberto Amo, Jacqueline Bloch and Natalia G. Berloff*

# Polaritonic XY-Ising machine

https://doi.org/10.1515/nanoph-2020-0162

Received March 2, 2020; accepted April 15, 2020; published online

June 18, 2020

**Abstract:** 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.

**Keywords:** analogue computing machine; exciton-polaritons; Ising Hamiltonian; polaritonic XY-Ising machine; XY Hamiltonian; unconventional computing.

## 1 Introduction

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) [1–3], a CMOS-based Fujitsu digital annealer [4], an all-electronic oscillator-based Ising machine [5], a simulated bifurcation Toshiba machine [6], a spatial light modulator (SLM) based photonic Ising machine [7], an optical simulator with the injection-locked multicore fibre lasers [8], and a molecular computing machine on a programmable microdroplet array [9]. 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 [10], 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 [11, 12] while the latter facilitated the design of instances of interaction matrices with the controlled complexity and planted solutions [13, 14]. 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 [15], mode-locked fibre lasers [16], coupled laser [17], and nanolaser arrays [18].

*Corresponding author: Natalia G. Berloff, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, CB3 0WA, UK; and Skolkovo Institute of Science and Technology, Bolshoy Boulevard 30, Build. 1, Moscow, 121205, Russian Federation, E-mail: N.G.Berloff@damtp.cam.ac.uk. https://orcid.org/0000-0003-2114-4321

Kirill P. Kalinin: Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, CB3 0WA, UK

Alberto Amo: Univ. Lille, CNRS, UMR 8523 – PhLAM – Physique des Lasers Atomes et Molécules, Lille, F-59000, France

Jacqueline Bloch: Université Paris-Saclay, CNRS, Centre de Nanosciences et de Nanotechnologies, Palaiseau, 91120, France

<figure><img src="image_3.png" /></figure>
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 [19]. 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 $\theta_i$. These phases can be mapped into continuous 'spins' $s_i = (\cos\theta_i, \sin\theta_i)$ that can be further constrained to discrete values $\theta_i \in \{0, \pi\}$ by employing the resonant excitation [20]. 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 [21–24], in etched micropillars [25–27], in lead halide perovskite lattices [28, 29], in strain-induced traps [30], in hybrid air gap microcavities [31], by periodically etching the sample surface depositing metallic patterns [32–34], by surface acoustic waves [35], by direct fabrication with the gold deposition technique [36], in coupled mesas etching during the growth of the microcavity [37]. 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 [21, 22, 38–40] 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 [41] 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 [41]. 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 [40, 42] and later followed by propositions to control interactions in regular lattices of condensates with dissipative gates [43] or resonant pump barriers [44]. The initial experimental demonstration of minimising the XY spin Hamiltonian with simple polarion lattice cells [40] has been shortly followed by a theoretical proposal extending the class of simulated spin Hamiltonians from continuous to discrete [20].

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 [45] noting that a typical lattice size constant is often in the range of 5–15 µm [40]. 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 [41] 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.

## 2 Results

Over the last decade, polariton condensates have been successfully modelled [46–49] by the generalised cGLE coupled to the exciton reservoir dynamics

$$ i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} (1 - i\eta \mathcal{R}) \nabla^2 \psi + U_0 |\psi|^2 \psi + g_R \mathcal{R} \psi \\ + \frac{i\hbar}{2} [R_R \mathcal{R} - \gamma_C] \psi + i f_{\text{Res}}(r, t) \psi^* \quad (1) $$

$$ \frac{\partial \mathcal{R}}{\partial t} = -(\gamma_R + R_R |\psi|^2) \mathcal{R} + P(r, t), \quad (2) $$

where $\psi(r, t)$ is the condensate wavefunction, $m$ is the effective polariton mass, $\mathcal{R}(r, t)$ is the density distribution of the exciton reservoir, $U_0$ and $g_R$ are the polariton-polariton and polariton-exciton interaction strengths respectively, $R_R$ is the rate of scattering from the reservoir into the condensate, $\eta$ is the energy relaxation. The $f_{\text{Res}}$ 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 $\pi$ [20]. Linear photon losses in the cavity and exciton losses are described by the condensate ($\gamma_C$) and the reservoir ($\gamma_R$) damping rates. The pumping intensity $P(r, t)$ 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 [23, 24, 26, 50–55]. The optimisation accuracy of either XY or Ising models does not change for small values of $\eta$ that are assumed in experiments [20] 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 $\psi \to \sqrt{\hbar R_R / 2U_0 l_0^2} \psi$, $t \to 2l_0^2 t / R_R$, $r \to \sqrt{\hbar l_0^2 / (mR_R)} r$,

$\mathcal{R} \to \mathcal{R} / l_0^2$, $P \to R_R P / 2l_0^2$, $f_{\text{Res}} \to 2l_0^2 f_{\text{Res}} / \hbar R_R$, $g = 2g_R / R_R$, $b_0 = 2\gamma_R l_0^2 / R_R$, $b_1 = \hbar R_R / U_0$, $\gamma = \gamma_C l_0^2 / R_R$, $l_0 = 1 \mu\text{m}$, we obtain

$$ i \frac{\partial \psi}{\partial t} = - \nabla^2 \psi + |\psi|^2 \psi + g \mathcal{R} \psi + i (\mathcal{R} - \gamma) \psi + i f_{\text{Res}}(r, t) \psi^*, \quad (3) $$

$$ \frac{\partial \mathcal{R}}{\partial t} = -(b_0 + b_1 |\psi|^2) \mathcal{R} + P(r, t). \quad (4) $$

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 [27]. 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 [56] to the gap solitons in 1D Lieb lattices [57]. To model the polariton condensation in a micropillar cavity, we introduce a spatially dependent dissipative profile,

$$ \gamma(r) = \gamma_{\text{out}} - (\gamma_{\text{out}} - \gamma_{\text{in}}) \sum_i \exp(-\alpha|r - r_i|^{2n_{\text{SG}}}), \quad (5) $$

where $\gamma_{\text{out}}$ and $\gamma_{\text{in}}$ are the dissipation rates outside and inside of a micropillar, respectively. Here, $\gamma_{\text{out}} \gg \gamma_{\text{in}}$, $r_i$ denotes the center of the i-th micropillar, and $n_{\text{SG}}$ is the degree of a supergaussian that models micropillars as flat low-dissipative discs. The dramatically increased dissipation between the discs ($\gamma_{\text{out}} = 100\gamma_{\text{in}}$) 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

$$i\psi_t = -\nabla^2\psi + |\psi|^2\psi + g\mathcal{R}\psi + i(\mathcal{R} - \gamma)\psi + if_{\text{Res}}\psi^* \\ +i\delta_{\gamma, y_{\text{in}}} \sum_{j=1, j \neq i}^{N} \mathcal{J}_{ij}\psi(r + r_j - r_i, t - \tau) \quad (6)$$

where $\delta_{\gamma, y_{\text{in}}}$ is the delta-function which is equal to one inside a micropillar, i. e. when $\gamma(r) = \gamma_{\text{in}}$, and zero outside, $N$ is the number of micropillars, and $\tau$ represents a possible time-delay to supply couplings in an experimental setup. The coupling term represents the emission feedback when for each $\psi(r)$ in a micropillar $i$ the respective values $\psi(r + r_j - r_i)$ are added from the micropillar centred at $j$. For the *relative* coupling model we shall consider $\mathcal{J}_{ij} = J_{ij}$ while for the *absolute* coupling model we will use $\mathcal{J}_{ij} = |\psi_i|J_{ij}/|\psi_j|$. The sign of the coupling strength can be made positive or negative by injecting the light with zero (ferromagnetic coupling) or $\pi$ phase (anti-ferromagnetic coupling), respectively. For further derivations, we denote $K(r) = \Theta(R - r)$ 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 [58]. 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 [59, 60]. For negligible time-delay and geometric couplings between condensates, we can rewrite Eqs. (4) and (6) for each micropillar $i$ using $\psi = \sum_i \psi_i$ and $\mathcal{R} = \sum_i \mathcal{R}_i$ as N equations for the polariton condensates $\psi_i = K(|r - r_i|)\psi(r - r_i, t)$ and $N$ equations for reservoir densities $\mathcal{R}_i = K(|r - r_i|)\mathcal{R}(r - r_i, t)$ noting that $P(r, t) = \sum_i P_i = \sum_i P(|r - r_i|, t)$, $f_{\text{Res}}(r, t) = \sum_i f_{\text{Res}}^{(i)} = \sum_i f_{\text{Res}}(|r - r_i|, t)$:

$$i\partial_t \psi_i = -\nabla^2 \psi_i + |\psi_i|^2 \psi_i + g \mathcal{R}_i \psi_i + i(\mathcal{R}_i - \gamma_{\text{in}}) \psi_i + if_{\text{Res}}^{(i)} \psi_i^* \\ + i \sum_{j=1, j \neq i}^{N} \mathcal{J}_{ij} \psi_j \quad (7)$$

$$\partial_t \mathcal{R}_i = -(b_0 + b_1 |\psi_i|^2) \mathcal{R}_i + P_i. \quad (8)$$

The steady states of Eqs. (7) and (8) correspond to the minima of the XY, i. e. $f_{\text{Res}}^{(i)} = 0$, and Ising, i. e. $f_{\text{Res}}^{(i)} \neq 0$, models as it becomes evident after we substitute

$\psi_i = \sqrt{\rho_i} \exp[i\theta_i]$ into Eq. (8) and separate the real and imaginary parts. The equations read as

$$\frac{1}{2} \partial_t \rho_i = (\mathcal{R}_i - \gamma_{\text{in}}) \rho_i + \sum_{j=1, j \neq i}^{N} \sqrt{\rho_j \rho_i} \mathcal{J}_{ij} \cos(\theta_{ji}) \\ + \rho_i f_{\text{Res}}(r, t) \cos(2\theta_i), \quad (9)$$

$$\partial_t \theta_i = \frac{\nabla^2 \sqrt{\rho_i}}{\sqrt{\rho_i}} - \rho_i - g \mathcal{R}_i + \sum_{j=1, j \neq i}^{N} \sqrt{\frac{\rho_j}{\rho_i}} \mathcal{J}_{ij} \sin(\theta_{ji}) \\ - f_{\text{Res}}(r, t) \sin(2\theta_i), \quad (10)$$

$$\partial_t \mathcal{R}_i = -(b_0 + b_1 \rho_i) \mathcal{R}_i + P_i (|r - r_i|, t), \quad (11)$$

where $\theta_{ji} = \theta_j - \theta_i$. Here we considered the uniform phase distribution $\theta_i(|r - r_i|, t) \approx \theta_i(t)$ which is a valid assumption near the micropillar's centre, i. e. $R < R_m$ with $R_m$ 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 [58] across all micropillars, that is, when the condition $\rho_i(r) = \rho_j(r)$ stands. Such density equilibration can be robustly achieved by iteratively updating pumping intensity $P_i$ so that $\int \rho_i dr = \rho_0$ for all micropillars, where $\rho_0$ 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 $\mathcal{J}_{ij} = |\psi_i| J_{ij} / |\psi_j|$ 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:

$$(R_i - \gamma_{\text{in}}) = - \sum_{j=1, j \neq i}^{N} J_{ij} \cos(\theta_{ji}) - f_{\text{Res}}^{(i)} \cos(2\theta_i), \quad (12)$$

$$\mu - \frac{\nabla^2 \sqrt{\rho_i}}{\sqrt{\rho_i}} + \rho_i + g \mathcal{R}_i = \sum_{j=1, j \neq i}^{N} J_{ij} \sin(\theta_{ji}) - f_{\text{Res}}^{(i)} \sin(2\theta_i), \quad (13)$$

$$R_i = P_i (b_0 + b_1 \rho_i)^{-1}, \quad (14)$$

where $\mu$ 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:

$$\max \sum_{i=1}^{N} \int_{\Omega} \rho_i dr \Leftrightarrow \min \sum_{i=1}^{N} \int_{\Omega} \mathcal{R}_i dr \Leftrightarrow \min H_{XY|Ising}$$

$$H_{XY|Ising} = -\frac{1}{2} \sum_{i,j=1}^{N} J_{ij} \cos(\theta_{ij}) - \sum_{i=1}^{N} \left( \int_{\Omega} f_{\text{Res}}^{(i)} dr \right) \cos(2\theta_i),$$
<figure><img src="image_4.png" /></figure>

**Figure 1:** Top: Phase difference as a function of coupling strength for a polariton dyad. The Eqs. (7) and (8) are simulated for 50 random initial conditions for each coupling value. The coherence occurs for the absolute values of strengths greater than 0.02 leading to ferromagnetic state with 0 phase difference for positive couplings and to antiferromagnetic state with $\pi$ phase difference for negative couplings. The slowly decaying unstable solutions are shown in grey. Bottom: Phase difference as a function of time-delay for a polariton dyad. The time-delay percentage is defined with respect to the time required to reach a steady state in the absence of delay. The scatter point size indicates how many states out of 50 initial conditions end with a particular phase difference. The coupling strength between condensates is chosen to be $J_0 = -0.1$. The expected anti-ferromagnetic state is observed for time-delays $\tau < 2\%$ and followed by the region with decoupled condensates. The further phase-locking of condensates becomes possible for bigger time-delays due to the global phase presence.

where $\Omega$ denotes the plane of the microcavity. The resonant force term $f_{\text{Res}}^{(i)}$ 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 $gR_i$ 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 $R_i$ 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 $|J_0| \le 0.02$ 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 $\pi$ 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 $|J_0| \ge 0.07$. 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 $\tau$ 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 $t - \tau$ 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 $J_0 = -0.1$. 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 $\tau$ up to 2%. The previously unstable local minimum, i. e. ferromagnetic
state with 0 phase difference for $J_0 = -0.1$, 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 $\tau > 7\%$. This peculiar synchronization behaviour can be attributed to the global phase rotation with frequency $\mu$ of each condensate which can lead to phase-locking of condensates with an additional $\pi$ 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 [41] in Eqs. (14)–(12) which for significant $\tau$ 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 $\gamma_{out} = 100$ outside nodes compared to low $\gamma_{in} = 1$ values inside nodes. A random interaction matrix is constructed of positive and negative couplings of amplitude $\{0.05, 0.1\}$ 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

$(f_{Res} = 0)$. In the former case, the densities of condensates are iteratively equalised over time by individually adjusting pumping intensities $P_i$. 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:

$$ \text{Median Accuracy} = \left\langle \frac{H_{\text{Relative|Absolute}}}{H_{\text{Ground State}}} \right\rangle. \qquad (15) $$

where $H_{Relative|Absolute}$ 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, $H_{Ground State}$ is the ground state energy found by the classical optimisation algorithms. In Figure 2(d–e), the found minima are within 1% and 0.4% from the ground state of the XY Hamiltonian that was verified with the gain-dissipative [11] and the basin-hopping [61] algorithms (Figure 2(f)). The median accuracy over 100 random fully-connected matrices of size $N = 9$ generalises to 99.2% and 99.5% 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

<figure><img src="image_5.png" /><figcaption>the XY Hamiltonian is verified by the gain-dissipative and the basin-hopping algorithms.</figcaption></figure>

<caption>**Figure 2:** Finding the global minimum of the XY Hamiltonian of size $N = 9$ with a 3 × 3 polariton lattice by simulating Eqs. (7) and (8). (a) The intensity distribution of the incoherent pumping profile $P(r)$. The condensate emissions within the black circles are used for couplings between condensates. (b) The dissipative profile for realising spatially isolated polariton pillars. (c) The fully-connected coupling matrix $J$ which is randomly constructed from $J_1 = 0.05$ and $J_2 = 0.1$ of random signs. (d–e) The polariton density profiles and phase configurations are plotted for relative and absolute coupling models, respectively. The white arrows represent the phase difference with respect to the central condensate (with respect to the vertical arrow). The corresponding energy values of the XY Hamiltonian are shown in the top-right corner. (f) The ground state solution of</caption>
<figure><img src="image_6.png" /></figure>

**Figure 3:** Optimal amplitude range study for *relative* and *absolute* coupling models on the unweighted MaxCut problems of size *N* = 16 with degrees 5, 9, and 13. The median accuracy is shown for the XY Hamiltonian in (a–b) and the Ising Hamiltonian in (c–d). Both models are simulated with Eqs. (7) and (8) for 20 random initial conditions per each coupling strength. Shading indicates 25th and 75th percentile range of instances.

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 [62] and can be mapped to the Ising Hamiltonian by assigning antiferromagnetic couplings $J_{ij} = -1$ to the graph edges. We construct three such random adjacency matrices $A$ of size $N = 16$ of degree 5, 9, and 13. Both coupling models are simulated on matrices $J = -J_0 A$ with amplitudes $J_0$ in the range [0.001, 0.3]. 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 $J_0$ amplitude for the XY Hamiltonian. The optimal range of couplings with the median accuracy over 90% can be identified for the amplitudes in the range [0.01, 0.11] for the *relative* coupling model and slightly smaller range of [0.02, 0.1] for the *absolute* coupling model. For the Ising model, a smaller batch of coupling amplitudes allows one to achieve the median accuracy greater than 90% (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 $J_0$ for achieving higher optimisation accuracies for both models. We note that the ground states of the Ising Hamiltonians were verified with the gain-dissipative [11] and CIM [12] 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 50%. For such connectivity, we pick the amplitude strength of $J_0 = 0.04$ 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 98% 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 $J_0$ 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.
<table><caption>**Table 1:** Optimisation of the Ising and XY spin Hamiltonians with *relative* and *absolute* coupling models on unweighted MaxCut problems of size 25 and 49 with edge density 0.5. The median accuracy of both models is calculated for 20 random initial conditions per each coupling matrix which was further averaged over 20 random coupling matrices with coupling strength $J_0 = 0.04$. The number in parentheses indicates how many problems with different coupling matrices were globally optimised. The ground state solutions are calculated with the gain-dissipative and the basin-hopping algorithms for the XY Hamiltonians and the gain-dissipative and CIM algorithms for the Ising Hamiltonians.</caption><thead><tr><th rowspan="2">Problem Size</th><th colspan="2">Relative</th><th colspan="2">Absolute</th></tr><tr><th>XY</th><th>Ising</th><th>XY</th><th>Ising</th></tr></thead><tbody><tr><td>25 (5 × 5 lattice)</td><td>99.3% (20)</td><td>87.8% (20)</td><td>96.8% (20)</td><td>72.9% (14)</td></tr><tr><td>49 (5 × 5 lattice)</td><td>98.2% (20)</td><td>81.7% (4)</td><td>93.3% (3)</td><td>52.3% (0)</td></tr></tbody></table>

# 3 Discussion

## 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 $J_{ij}$. 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 $J_{ij}$ 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 [63].

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. $J_{ij} = J_{ji}$ 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 [11]. 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 [58]. 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 [64] and the recently proposed pure molecular simulator [9]. 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 [2] or hybrid molecular simulator [9].
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 [13]. 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 [16] and Ising [8] 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 $\eta = 0$, $g = 0.1$, $b_0 = 1$, $b_1 = 20$, $P(r, t) = \sum_i P_i \exp(-A \cdot |r - r_i|^2)$ with $P_i = P_0 = 10$ for all micropillars in case of the absolute coupling scheme and dynamically adjusted $P_i$ to bring all the condensates to $\rho_0 = 1$ in case of the relative coupling scheme for the XY Hamiltonian, $A = 5$, the
distance between micropillars is $d = 2.4$, $γ_{in} = 1$, $γ_{out} = 100$, $α = 10$, $n_{SG} = 10$, the micropillar diameter is about $d_{micropillar} = 2$, $R ≈ 0.25$. In addition, the following parameters are used to simulate the resonant pumping: $f_{Res}$ ($r, t$) = 50 (tanh(0.1$t$ - 3) + 1) $\sum_{i=1}^{N} \exp(-25|r - r_i|)$, where $t \in [0, T_{max}]$, $T_{max}$ is the time required to achieve a steady state.

**Conflict of interest statement:** The authors declare that they have no competing interests.

**Acknowledgements:** KPK acknowledges the support from EPSRC and Cambridge Trust. NGB acknowledges the support from Huawei.

## References

[1] 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.

[2] 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.

[3] 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.

[4] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H. Katzgraber, "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.

[5] 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.

[6] 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.

[7] 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.

[8] 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.

[9] S. Y. Guo, P. Friederich, Y. Cao, et al. "A molecular computing approach to solving optimization problems via programmable microdroplet arrays," *ChemRxiv*, 2018.

[10] A. Lucas, "Ising formulations of many NP problems," *Front. Phys.*, vol. 2, p. 5, 2014 Feb 12, https://doi.org/10.3389/fphy.2014.00005.

[11] 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.

[12] 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.

[13] 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.

[14] 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.

[15] 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.

[16] 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.

[17] 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.

[18] 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.

[19] J. Kasprzak, M. Richard, S. Kundermann, et al., "Bose-Einstein condensation of exciton polaritons," *Nature*, vol. 443, no. 710, pp. 409–414, 2006 Sep, https://doi.org/10.1038/nature05131.

[20] 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.

[21] 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.

[22] 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.

[23] 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.

[24] 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.

[25] 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.

[26] 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.

[27] 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.

[28] 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.

[29] 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.

[30] 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.

[31] 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.

[32] 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.

[33] A. Mischok, V. G. Lyssenko, R. Brückner, et al., "Zero and $n$ 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.

[34] 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.

[35] 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.

[36] 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.

[37] 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.

[38] 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.

[39] 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.

[40] 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.

[41] 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.

[42] 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.

[43] 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.

[44] 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.

[45] 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.

[46] 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.

[47] 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.

[48] J. Keeling and N. G. Berloff, "Excitonpolariton condensation," *Contemp. Phys.*, vol. 52, no. 2, pp. 131–151, 2011 Mar 1, https://doi.org/10.1080/00107514.2010.550120.

[49] I. Carusotto and C. Ciuti, "Quantum fluids of light," *Rev. Mod. Phys.*, vol. 85, no. 1, p. 299, 2013 Feb 21, https://doi.org/10.1103/revmodphys.85.299.

[50] 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.

[51] 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.

[52] 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.

[53] 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.

[54] 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.

[55] 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.

[56] 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.

[57] 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.

[58] 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.
[59] 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.

[60] 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.

[61] 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.

[62] R. M. Karp, "Complexity of computer computations," *Reducibility among Comb. Probl.*, vol. 23, no. 1, pp. 85-103, 1972.

[63] 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.

[64] 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.