Polaritonic XY-Ising machine

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

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 θi. These phases can be mapped into continuous ‘spins’ si=(cosθi, sinθi) that can be further constrained to discrete values θi∈{0, π} 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], [22], [23], [24], in etched micropillars [25], [26], [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], [33], [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], [39], [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 Hamitonians 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], [47], [48], [49] by the generalised cGLE coupled to the exciton reservoir dynamics

where ψ(r, t) is the condensate wavefunction, m is the effective polariton mass, R(r, t) is the density distribution of the exciton reservoir, U0 and gR are the polariton-polariton and polariton-exciton interaction strengths respectively, RR is the rate of scattering from the reservoir into the condensate, η is the energy relaxation. The fRes 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 π [20]. Linear photon losses in the cavity and exciton losses are described by the condensate (γC) and the reservoir (γ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], [51], [52], [53], [54], [55]. The optimisation accuracy of either XY or Ising models does not change for small values of η that are asumed 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 ψ→ℏRR/2U0l02ψ, t→2l02t/RR, r→ℏl02/(mRR)r, R→R/l02, P→RRP/2l02, fRes→2l02fRes/ℏRR, g=2gR/RR, b0=2γRl02/RR, b1=ℏRR/U0, γ=γCl02/RR, l₀ = 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 [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,

where γout and γin are the dissipation rates outside and inside of a micropillar, respectively. Here, γout≫γin,ri denotes the center of the i-th micropillar, and nSG is the degree of a supergaussian that models micropillars as flat low-dissipative discs. The dramatically increased dissipation between the discs (γout=100γ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

where δγ,γin is the delta-function which is equal to one inside a micropillar, i. e. when γ(r)=γin, 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 ψ(r) in a micropillar i the respective values ψ(r+rj−ri) are added from the micropillar centred at j. For the relative coupling model we shall consider Jij=Jij while for the absolute coupling model we will use Jij=|ψi|Jij/|ψj|. 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 K(r)=Θ(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 ψ=∑iψi and R=∑iRi as N equations for the polariton condensates ψi=K(|r−ri|)ψ(r−ri, t) and N equations for reservoir densities Ri=K(|r−ri|)R(r−ri, t) noting that P(r, t)=∑iPi=∑iP(|r−ri|, t), fRes(r, t)=∑ifRes(i)=∑ifRes(|r−ri|, t):

The steady states of Eqs. (7) and (8) correspond to the minima of the XY, i. e. fRes(i)=0, and Ising, i. e. fRes(i)≠0, models as it becomes evident after we substitute ψi=ρiexp[iθi] into Eq. (8) and separate the real and imaginary parts. The equations read as

where θji=θj−θi. Here we considered the uniform phase distribution θi(|r−ri|,t)≈θi(t) which is a valid assumption near the micropillar’s centre, i. e. R<Rm with Rm 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 ρi(r)=ρj(r) stands. Such density equilibration can be robustly achieved by iteratively updating pumping intensity Pi so that ∫ρidr=ρ0 for all micropillars, where ρ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 Jij=|ψi|Jij/|ψ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:

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 fRes(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 gRi 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 Ri 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 |J0|≲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 π 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 |J0|≥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.

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 π 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 J0=−0.1. The expected anti-ferromagnetic state is observed for time-delays τ<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.

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 t−τ 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 J0=−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 τ up to 2%. The previously unstable local minimum, i. e. ferromagnetic state with 0 phase difference for J0=−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 τ>7%. 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 [41] 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 γout=100 outside nodes compared to low γ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 (fRes=0). In the former case, the densities of condensates are iteratively equalised over time by individually adjusting pumping intensities Pi. 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 HRelative|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, HGround 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.

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 J1=0.05 and J2=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 the XY Hamiltonian is verified by the gain-dissipative and the basin-hopping algorithms.

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 [62] and can be mapped to the Ising Hamiltonian by assigning antiferromagnetic couplings Jij=−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=−J0A with amplitudes J0 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 J0 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 J0 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.

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.

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 J0=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 J0 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.

3 Discussion

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 η=0, g=0.1, b0=1, b1=20, P(r, t)=∑iPiexp(−A⋅|r−ri|2) with Pi=P0=10 for all micropillars in case of the absolute coupling scheme and dynamically adjusted Pi to bring all the condensates to ρ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, nSG=10, the micropillar diameter is about dmicropillar=2, R≈0.25. In addition, the following parameters are used to simulate the resonant pumping: fRes(r, t)=50(tanh(0.1t−3)+1)∑i=1Nexp(−25|r−ri|) , where t∈[0, Tmax], Tmax is the time required to achieve a steady state.