Antiferromagnetic Ising model in a triangular vortex lattice of quantum fluids of light

Vortices are topologically distinctive objects appearing as phase twists in coherent fields of optical beams and Bose-Einstein condensates. Structured networks and artificial lattices of coupled vortices could offer a powerful platform to study and simulate interaction mechanisms between constituents of condensed matter systems, such as antiferromagnetic interactions, by replacement of spin angular momentum with orbital angular momentum. Here, we realize such a platform using a macroscopic quantum fluid of light based on exciton-polariton condensates. We imprint all-optical hexagonal lattice that results into a triangular vortex lattice, with each cell having a vortex of charge l = ±1. We reveal that pairs of coupled condensates spontaneously arrange their orbital angular momentum antiparallel, implying a form of artificial orbital “antiferromagnetism.” We discover that correlation exists between the emergent vortex patterns in triangular condensate lattices and the low-energy solutions of the corresponding antiferromagnetic Ising system. Our study offers a path toward spontaneously ordered vortex arrays with nearly arbitrary configurations and controlled couplings.


Introduction
Arrays of quantized vortices, appearing as phasesingular twists in wavefunctions, exemplify fascinating self-organising phenomena, originally studied in macroscopic interacting Bose gases such as rotating superfluids and Bose-Einstein condensates displaying nucleation and ordering of single-charge vortices 1,2 .More recently, optical vortex arrays in lasers have also gained increased interest [3][4][5][6] due to the promising applications of optical vortices 7 and phase-singular optics 8,9 in particle manipulation, high resolution imaging, and quantum communication protocols 10,11 .While optical systems offer superior spatial programmability over said vortex arrays 6 they do not possess the large interaction strengths inherent to quantum fluids, which play a crucial role in phase transitions and pattern formation 12 .For this reason, numerous studies have been focused on hybrid systems, which combine the best properties of light and matter to explore vortices 13 .Here, we investigate such a light-matter platform in the strong-coupling regime based on cavity exciton-polariton (hereafter polariton) condensates 14 offering a compromise between optical programmability 15 and intrinsic evolution of vortex arrays in driven-dissipative quantum fluids [16][17][18][19][20] .
Polariton condensates 14 are bosonic quantum fluids, characterised by large interaction strengths and light effective mass, that can be all-optically driven into artificial lattices using structured pump patterns [21][22][23] .A notable feature of polariton condensates is their nonequilib-rium superfluid character [24][25][26] , and ability to form vortices that are pinned by sample disorder 27 or through optical trapping 19,28,29 , or spun-up 30 .To date, polariton vortex arrays have been observed in the interference pattern of multiple condensate modes 17,20 or in specially patterned cavities 18 .However, large-scale programmable vortex lattices with tunable coupling strengths, where each vortex site can be individually addressed, remain unexplored in polariton fluids.
As a particular example, we explore here the concept of driven-dissipative geometric frustration by designing a triangular lattice of "antiferromagnetically" (AFM) coupled polariton vortices.Here, AFM refers to polariton orbital-angular-momentum (OAM) instead of spinangular-momentum. Conventionally, geometric frustration can be described as an inability of a system to minimize its real energy through reduction of pairwise interactions between constituent elements 31 .A prominent example are Ising spins arranged into an AFM triangular graph 32 , as shown schematically by the black arrows in Fig. 1A.Only two pairwise interactions can be minimized at any time leaving the third at a higher energy.So far, artificial spin ice systems are popular candidates to explore frustration at large scale 33 but recently lattices in bosonic systems have steered towards this direction using ultracold atoms 34 , coupled lasers 35 and excitonpolariton condensates [36][37][38][39] .However, instead of minimizing their real energy, polariton condensates minimize their losses (i.e.maximize their imaginary energy component) when they ballistically couple over mutually pumped region [37][38][39] .From this perspective, one can define frustration for polaritons as their inability to minimize losses through non-Hermitian interactions between lattice nodes 40 .
In this article, we report on the evidence of frustration in a driven-dissipative triangular lattice of polariton condensates.To realize this, we optically imprint a lattice of trapped polariton condensates that populate a vortex state at each lattice node, as shown schematically in Fig. 1B.Instead of spin-angular-momentum, we work with condensates carrying OAM, with topological charges l = ±1, defined by a superposition of the traps first excited dipole modes |P x ⟩ ± i|P y ⟩.We show that we can tune the coupling between synchronized vortex pairs to favour either parallel or antiparallel OAM.The latter constitutes an effective AFM coupling mechanism which results in frustration when triangular geometry of condensates is introduced, due to incommensurate symmetries between the C 6 (sixfold rotational symmetry of the triangular lattice) and SO(4) (rotation of the twocomponent OAM condensate) operators.Scaling up to a large 22-vortex system we observe spontaneous formation of non-periodic vortex configurations above threshold, changing from realization to realization, implying existence of multiple available modes.We provide experimental evidence that observed vortex patterns correlate significantly with the ground state solutions of the corresponding Ising system.The large nonlinearities inherent to exciton-polariton systems are crucial for the stability of the vortices, facilitating the readout of their OAM configurations even over long excitation pulses (≈ 10 7 longer than the polariton lifetime).

Results
We use a planar 2λ GaAs-based microcavity with embedded InGaAs quantum wells 41 held at ≈4 K in a closed-cycle helium cryostat.Using a reflective phaseonly spatial light modulator (SLM), we transform an incident nonresonant single-mode continuous-wave laser beam into an ordered honeycomb array (the dual of the triangular lattice) of small Gaussian beams focused onto the microcavity plane (see schematic in Fig. 1B) 15 .The nonresonant circularly polarized pump is tuned at the first Bragg minimum of the microcavity reflectivity stopband to avoid heating of the sample (1.5578 eV).As we show in the Supplementary Material, at the pump power used in experiments, trapped polariton condensate inherits the polarization state from the pump even under nonresonant excitation 42 .This allows us to excite a single energy state corresponding to polaritons with a predominant pseudo-spin and make therefore, a "clean" system without energy splitting (between σ + and σ − components).
The short distance between the pump spots (d = 6.6 µm) is chosen so as to obtain trapped polariton condensates in the centre of each cell (i.e. in the region between six pumping spots) 22 .Physically, polaritons become trapped in this region due to their strong repulsive interactions with the pump-injected background exciton reservoir [43][44][45] .In Supplementary Note 1 we show that for larger honeycomb cells (d = 8.8 µm) the polaritons mostly condense on top of the pump spots instead of becoming trapped, underlining the importance of tuning the pump parameters.

Vorticity onset in a single trap
In order to construct spatially extended triangular structure with coupled vortex states, we first implement and analyze their formation in a single, a pair, and a triangle of cells.In particular, the single cell experiment allows us to confirm that the formation of vorticity is truly a spontaneous symmetry-breaking event with a random sign of vortex charge from realization-to-realization.luminescence (PL) with a donut shape intensity profile for a single cell pumped at P = 1.2P thr .
The formation of a macroscopic wavefunction Ψ implies irrotational flow everywhere except at the vortex core, where the phase singularity is located and particle density vanishes.Phase winding around the core is quantized in units of 2πl, where l ∈ Z.In order to verify this, we apply a homodyne interferometric technique 46 to extract the phase of the condensate wavefunction under single-shot excitation conditions (i.e. the sample is excited by a single pulse of 50 µs duration).Our measurements of the single-shot phase maps across different realizations, where vortices form, reveal that 45% of the shots resulted in formation of a vortex (l = +1, Fig. 1D), and 55% in an antivortex (l = −1, Fig. 1E).Due to finite negligible pump anisotropy and sample disorder, we also find that 18 realizations out of 100 result in distributions with step-like π phase jumps, indicating the formation of a dipolar state, rather than of single vortex state, as it occurs in other 82 realizations, as described in Supplementary Note 2. The absence of a predominant sign of the observed vortex charges means that imprinted excitation pattern does not break chiral symmetry and that vorticity emerges spontaneously during the polariton condensation process.
In the Materials and Methods section we show, using a simple variational generalised Gross-Pitaevskii model, that l = ±1 vortices in an ideal cylindrically symmetric optical trap (approximating the actual trap induced by six pump spots) are the only stable condensate solutions, in qualitative agreement with our observations.We emphasise that our nonresonant excitation pattern (shown in fig.S3) does not carry any OAM, nor does it imprint any phase onto the condensate, in contrast to resonant excitation schemes 47 .

Two coupled vortices
Next, we focus on two neighbouring cells, shown in Fig. 2, which convincingly demonstrate signatures of coupling and synchronization.Figures 2A-D show measured time-averaged real-space polariton PL as a function of pump power.Clearly below (Fig. 2A) and at condensation threshold (Fig. 2B) the polariton PL does not reveal univocal signatures of vortex formation.However, at slightly higher power (P = 1.06P thr ) the condensate collapses inside the traps forming dipoles (Fig. 2C), oriented head-to-tails (with parallel nodal lines) in a σ-bonding fashion as predicted recently 48 .Then, above some vortex-threshold power the condensates in both cells reveal a donut-shaped intensity distributions (Fig. 2D), indicating the onset of vorticity.
In order to verify that the coupling in this configuration is effectively AFM, we extract the phase maps of the vortex pairs at P = 1.2P thr (see PL in Fig. 2E).Figures 2F-H show three typical examples of corresponding phase maps observed in a hundred independent singleshots.Remarkably, we obtained that 93% of the condensate realizations formed with opposite OAM between the cells.The remaining 7% correspond to vortex-dipole (or antivortex-dipole) states.Therefore, our data analysis conclusively proves that the vortex coupling between cells strongly favours AFM order.Physically, this can be understood as an optimal constructive interference condition between the traps which maximizes the gain of the condensate pair 48 .
The advantage of our optical technique is that by changing some parameters of the pump profile we can force the system to occupy vortex-vortex state (i.e., OAM "ferromagnetic" coupling) instead of vortex-antivortex states.To implement this, we optically imprinted weak potential barriers in the middle of the 2-cell structure in the same spirit as in demonstrated earlier coupling control approach 46 .As a result, we managed to demonstrate reversible switching the system from vortex-antivortex to vortex-vortex state as described in Supplementary Note 4.

Three vortices in a triangle
Figure 3A shows the cavity PL after we optically construct an equilateral triangle of trapped condensates.The system's preference for the AFM ordering as well as signatures of the expected topological vortex charge frustration is observed through homodyne interferometric measurements.Details on the interferometric technique are given in Supplementary Note 3 with an example of interference pattern given in fig.S5.In Fig. 3B-D we provide few examples of extracted phase maps corresponding to different configurations even though the excitation conditions (geometry and pump power) are kept the same.For example, in Fig. 3B one can see π phase jump lines implying formation of dipole states inside the cells in this particular realization.Dipole states can appear because of the sensitivity of polaritons to inherent cavity disorder or slight pump inhomogeneities which pin the dipoles along a specific direction 44 .We note that when a condensate is in a perfect superposition of vortex and antivortex modes then the resulting dipole state has zero net OAM and therefore, cannot be assigned to discrete variable of +1 or -1 following the Ising model.However, our further analysis indicates that when there is finite OAM in individual condensate sites, projecting it on classical binary variables shows pattern formation reminiscent of a frustrated triangular Ising system.None-the-less, for a major part of individual system's realizations (47% out of one hundred realizations), we find vortex-antivortex (antivortex-vortex) pair accompanied by a dipole state in the third cell, as shown in Fig. 3C.In 43% of condensate realizations we reveal vortex-antivortex-vortex (or antivortex-vortexantivortex) states in all cells, as depicted in Fig. 3D.We stress that we never observed vortices with the same OAM at the same time in all three cells.We also stress that Fig. 3C,D show examples over multiple possible arrangements of these type of patterns.These results, in the same vein as those of two cell structure, clearly indicate the tendency for development of the AFM configuration in each pair of neighbouring cells.The fact that multiple realizations display a cell with an indeterminate vortex (see top right cell in Fig. 3C) implies dynamic indecisiveness reminiscent of frustration.

Mean field modelling
Our experimental observations can be reproduced through numerical mean-field simulations, based on the generalised two-dimensional Gross-Pitaevskii equation (2DGPE) describing a macroscopic polariton wavefunction, Ψ(r, t), coupled to the rate equation for the exciton reservoir, n X (r, t) (see Materials and Methods).First, we analyze the vorticity onset in the cells of the lattice.In order to obtain families of possible solutions (modes) corresponding to one-, two-and three-cell pump configurations, we apply the Newton method and perform linear stability analysis.Families of stationary states Ψ(r, t) = w(r)e −iεt with ∂ t n X = 0 (here ε is the energy detuning) that emerge as stable attractors around threshold correspond to the coloured curves in Fig. 4A.In this graph we show the dependence of the scaled peak amplitude of the wavefunction (g c m eff r 2 0 /ℏ 2 ) 1/2 |Ψ| on pump power P 0 .Here, r 0 = 1 µm is the characteristic spatial scale, m eff is the polariton mass, and g c is the polaritonpolariton interaction strength.One can see that the condensation threshold power (starting point of the curves) decreases with increasing number of pumped cells as expected 49 .
On each solid curve, the onset of vorticity with pump power (i.e., finite net clockwise or anticlockwise rotation) is indicated with a blue star.For a single honeycomb pump cell (black curve), a stable vortex emerges almost at condensation threshold in good agreement with experiments.For the two-cell configuration (red curve), we find that the stable branch emerging from threshold corresponds indeed to vortex-antivortex pairs forming at a critical power given by the blue star, whilst below this power vorticity is absent and field distribution exhibits dipole structure.This critical pump power represents a vorticity threshold, observed experimentally between Fig. 2C and 2D.When we add a third pump cell, the first stationary family emanating from threshold is associated with a symmetric vortex-antivortex-dipole configuration appearing at the blue star on the green curve in Fig. 4A, in agreement with observations in Fig. 3A,C.For lower powers (below blue star) the condensate resembles a triple dipole state, whereas for high powers it gradually transforms into an asymmetric vortex-antivortexdipole solution (solid green curve with dots).Even further increase of pump power results in unstable behaviour (denoted with green dashed curve in Fig. 4A) and the emergence of different stationary states w(r), such as the antivortex-vortex-antivortex shown in Fig. 3D.Corresponding solutions around the blue stars on these three branches are shown in Fig. 4B-J obtained from direct numerical integration of the 2DGPE using white noise initial conditions and damped boundary conditions.
In the Materials and Methods section we additionally develop a simplified model describing the coupling between vortices by treating them as localised orbitals in each trap (i.e., tight binding approach).This reduces the problem of modelling a 2+1 nonlinear partial differential equation into a set of 2N ordinary differential equations (6).Here N is the number of traps and the factor "2" appears because of clockwise and anticlockwise currents.This effective "spinor" model allows us to analytically prove the stability of vortices in isolated traps through construction of a Lyapunov function.It also correctly predicts increased power of the vortex threshold when two traps are coupled together, in agreement with experiment and 2DGPE simulations.One advantage of the model (10) is that it possesses analogies with conventional spinor polaritons in patterned cavities with effective spin-orbit-coupling (subject to future work).

Triangular lattice of AFM coupled vortices: evidence of low energy Ising order
Finally, we study the vorticity patterns appearing in much larger spatially extended structures, consisting of 22 cells in a triangular geometry.Figure 5A shows polariton condensation into a macroscopic coherent state containing multiple vortices across the lattice for a pump power above condensation threshold unlike previous observations in square 15,37,50 , triangular 15 , and Lieb 22 opti-cal lattices.Similar to our observations in smaller structures, each cell possesses a donut shaped PL intensity profile.
Figure 5B shows extracted from experiment a singleshot phase map with multiple single-charge (l = ±1) phase singularities across the lattice, schematically marked by blue and red arrows.We underline that the excitation pulse width (50 µs) is much longer than the cavity lifetime (≈ 5.5 ps), indicating the topological robustness of the vortices.We observe that across multiple realizations (provided in Supplementary Note 5) the pattern of vortices changes as one might expect due to the spontaneous symmetry breaking upon condensation.
The apparently random configurations of vortices from realization-to-realization would suggest that the system is completely stochastic, which would be the case if the traps are very weakly coupled.However, as we know from observations on the two-cell (Fig. 2) and three-cell (Fig. 3) experiments, the traps are AFM coupled which should manifest in some form of AFM order in the lattice OAM with a highly degenerate "ground state" in analogy to Ising systems 32 .
To address this question, we analyzed the spatial OAM order of the vortices appearing in 25 single-shot realizations of the 22-cell structure, and checked for correlations with the low energy configurations of the Ising Hamiltonian, see Supplementary Note 8 for derivation of its role in picking optimal vortex arrangements, where J < 0 is the AFM coupling energy, the variables σ n ∈ {±1} represent Ising spins and the sum is taken over nearest neighbours in the triangular lattice.
The up/down OAM of each lattice node is then projected onto its corresponding Ising spin and classified into one of three groups 32 : (1) spin has more co-aligned neighbours ("FM" spins, E is raised); (2) spin has more counter-aligned neighbours ("AFM" spins, E is lowered); (3) spin has equal co-and counter-aligned neighbours ("free" spins, E is unchanged).Here, we refer the reader to small insets in Fig. 5C, on the right.The ground state of ( 1) can be organised in many different ways using conditionally ordered AFM spins and free spins.All 2 22 possible spin configurations and associated energies from equation ( 1) can be easily found by brute force.Binning each spin of each configuration into one of the three possible categories described above, we can plot the mean and standard deviation (SD) of the normalized number of spins in each category as a function of Ising energy E (see curves and shaded areas in Fig. 5C).As expected, the highest energy state has all spins parallel (only FM spins) with zero SD whereas the ground state has only AFM and free spins with finite SD (implying degeneracy).We then apply the same classification to the vortices in our experimental data (see fig.S8), and extract corresponding statistical occurrence of FM, AFM, and free vortices normalised against the number of detected vortices and plot on top of Fig. 5C.The three data points (as black, green and red empty markers) are placed where their error is minimum from the solid curves which puts them at E = −22J.Remarkably, our data points are within the lower ≈6% of Ising spin configurations from the ground state.From this we conclude, that our vortex lattice does not display stochastic arrangement of vortices, which would give E = 0, but instead tends towards AFM order.
Figures 5D and 5E show a representative example of the simulated condensate density and phase distributions emerging from random initial conditions, evidently illustrating the formation of a vortex lattice.Excitation from random noise typically yields dynamical, breathing, structures that nevertheless are characterised by persistence of vortices, once they form in the cells.In agreement with the experiment, the arrangement of vortices changes from simulation to simulation (see more examples in fig.S9), illustrating the highly multi-modal nature of the lattice, yet maintaining a dominant AFM order across the 22-cell structure.

Discussion
We note that self-arranging vortices are well-known in rotating superfluids 51 , atomic BECs 2 , and in nonlinear optical systems 3,4,6,7 .Here, we provide evidence that our all-optical lattice supports large scale polaritonic modes (or super-modes), which arise due to the inter-cell interactions in our nonlinear system, rather than to a collection of uncoupled individual single-site states.As such, system states are prone to develop strong instabilities that could demolish their robustness and, hence, their observation.For this reason, the observation of a large number of coupled vortices constitutes in itself a remarkable finding, that could not be anticipated a priori.With 22 cells (Fig. 5A,B) a large number of stationary states coexist for pump powers just above the condensation threshold and the system exhibits a large degree of multi-stability.As a consequence, the different experimental realizations with 22 cells reveal the formation of a wide diversity of states where individual vortex locations appear stochastically from one realization to the next.
Future investigations could verify whether it is possible to excite a macroscopic polariton system in a superposition of eigenstates carrying the vorticity and provide a clear answer to a fundamental question: is stochasticity of observed vortex states in such a system a solid proof of its multi-stability?In general, interferometry cannot unambiguously prove that the system is not in a quantum superposition of states within the same energy manifold.Therefore, to verify the hypothesis of a superposition of eigenstates, one would need to implement quantum vortex tomography, which requires experimental single-shot realization of a vortex sorter in free-space and tempo-ral correlation measurements with high resolution.The systematic study of extended lattices of coupled vortices, although challenging, it offers a promising avenue of research with an outlook to quantum simulators that requires thorough theoretical consideration.
In conclusion, we demonstrate evidence of an alloptically imprinted lattice of interacting vortices in nonresonantly driven-dissipative quantum fluids of light.Whereas in a single cell, we observe an intrinsically stochastic behaviour of an optically trapped polariton condensate, in the presence of the second neighbouring cell we reveal synchronization into a vortex-antivortex state, resembling the OAM analogue of AFM coupling.We also analytically predicted and experimentally demonstrated optical control of the vortex interactions from AFM to FM.Therefore, by tuning the interactions between ordered trapped polariton condensates, one can expect observation of various artificial magnetic phases between the nodes, emulating complex magnetic systems, and the study of other exotic coherent states in expanded structures.Our findings address the challenge of realizing the node-to-node coupling in spatially confined quantum fluids of light possessing OAM.

Materials and Methods
Generalised Gross-Pitaevskii model.The condensate dynamics is assumed to be described by the 2D generalised Gross-Pitaevskii equation (2DGPE) for the polariton wavefunction Ψ(r, t) coupled to a rate equation describing the exciton reservoir feeding the condensate n X (r, t) 52 : Our parameters are based on the sample properties and past experiments 41 .Here, m eff ≈ 5.63 × 10 −5 m e is the effective mass of the lower-branch polaritons where m e is the free electron rest mass, g c = 2.4 µeVµm 2 is the polariton-polariton and g r = 2g c is the polaritonreservoir interaction strengths typical for GaAs-based systems, R = 0.021 µm 2 ps −1 is fitted rate of stimulated scattering of polaritons from active reservoir, γ c ≈ 0.182 ps −1 and γ X ≈ 0.05 ps −1 , are the decay rates for polariton condensate and reservoir excitons, respectively, ηγ X ≈ 0.49 quantifies an additional blueshift coming from a background of dark reservoir excitons.The function P (r) = (P 0 γ X γ c /R) n Q(r − r n ) describes spatial pump profile consisting of identical Gaussian spots Q(r) = e −r 2 /d 2 of width d (corresponding to 2.5 µmwide pump spots used in the experiment) placed in the points with coordinates r n defining particular pump configuration (one or several honeycomb cells or large-scale honeycomb lattice with period D).Dimensionless pump amplitude P 0 is defined here in units of γ X γ c /R corresponding to condensation threshold for spatially uniform pump.
The system (2)-( 3) was solved using a variant of split-step fast Fourier transform method combined with fourth-order Runge-Kutta method to account for interaction between the condensate and the reservoir.To uncover all existing stationary states in the system and possible types of evolution, including essentially dynamical regimes, when excitation of stationary states does not occur, the modeling was initiated using different realizations of small-scale noise for initial Ψ, n X distributions.To identify the majority of all possible stationary solutions of the system or dynamical evolution regimes, multiple (up to several hundreds) input noise realizations were implemented for each pump power and pump-configuration.To account for considerable ballistic expansion of polaritons from the pumped regions leading to nonzero polariton density even far from them, we used in modeling sufficiently large spatial domain 300 × 300 µm 2 greatly exceeding the pumped area.Stationary states of the system (2)-( 3) corresponding to ∂ t n X = 0 were found in the form Ψ(r, t) = w(r)e −iεt , where w(r) is a complex function describing stationary polariton distribution in the cavity plane, ε is the energy detuning, using Newton iteration method.Due to the dissipative nature of our system, the energy detuning ε of stationary states, some of which appear as stable attractors, is also determined by the pump amplitude P 0 .Stability of such stationary states, whose families are presented in Fig. 4 for different pump configurations, were tested by modeling their evolution in the system (2)-( 3) and by performing linear stability analysis on them.
Single condensate vorticity onset.Here, we will write a variational form of the generalised Gross-Pitaevskii model describing only the dynamics of the first excited degenerate pair of l = ±1 OAM modes in the condensate.We will then proceed to prove, for an ideal cylindrically symmetric trap, that the only (and equally) stable condensate solutions correspond to these OAM modes.
We will assume for simplicity, that the cw driven excitonic reservoir follows the condensate adiabatically so that ∂ t n X ≃ 0 which implies, n X = P (r) . In the last step we have Taylor expanded the reservoir solution around small R|Ψ| 2 /γ X which is valid at pump powers not too far from threshold.This allows us to write a simpler 2DGPE, We are interested in a condensate that occupies the degenerate pair of clockwise and anticlockwise OAM states.
Assuming an ideal cylindrically symmetric trap P (r) = P (r) we can write our truncated basis as, Here, ξ(r) is the radial steady state (∂ t |Ψ| 2 = 0) part of the condensate in a single optical trap by solving the 2DGPE.Plugging in this truncated solution into (4) and integrating over the space, exploiting the orthogonality of the states, we arrive at two coupled equation of motion describing the vortex phase and occupation, up to an overall energy shift.Note the factor 2 in the counter-rotating nonlinearity (cross-Kerr term) which is orders of magnitude larger than the singlet interaction strength in spinor polariton condensates.The new coefficients are given by, p = P0γc 6) can be written in terms of the amplitude and phase of each mode We see that the change in the phase of the modes is trivially determined by the dynamics of their amplitudes.We therefore only need to focus on solutions of equation (7a) which has three equilibrium points: N ± = 0 and N ∓ = p/ R; and N + = N − = p/3 R. It is easy to show that the only (and equally) stable equilibrium points are the former through the eigenvalues λ of the Jacobian of equations (7a).This can also be visulized by plotting the system Lyapunov potential with extrema corresponding to these solutions, The Lyapunov potential satisfies the condition dL/dt ≤ 0 which means that all points in phase space flow towards the two minima indicated by the white dots in The colorscale depicts the Lyapunov potential (8) of the system with two minima (white dots) corresponding to vortex and antivortex solutions in which all phase space trajectories converge towards, underlining the deterministic nature of the system.Here we set p = R = 1 without loss of generality.(B) Reproduced results from Fig. 4A using the variational Gross-Pitaevskii model for a single ( 6) and a pair of coupled condensates (10).The blue stars denote the onset of vorticity for increasing power.
model evidences that the only stable condensate solutions in a single cell are vortices which is in agreement with experiment.
These vortex solutions bifurcate from the uncondensed normal phase N ± = 0 at condensation threshold of the single trap written, The above integral quantifies the amount of overlap between the condensate and the pumped region.If this overlap increases, then the threshold is lowered as expected.
Vorticity onset for two coupled condensates.The equations describing two coupled trapped condensates read, where n = 1, 2 denote the left and right condensate.
Here, J a,b ∈ C are the non-Hermitian tunneling rates between co-rotating and counter-rotating vortices (i.e.OAM conserving and non-conserving coupling strengths, respectively), , where Here, Q tot (r) describes the two pumped traps and the separation between them can be written r ′ −r = dx.The coordinates of the two pumps are related through r ′ = |r ′ | = r 2 + d 2 − 2rd cos (θ) and sin (θ ′ ) = r sin (θ)/r ′ .Physically, Re(J a,b ) and Im(J a,b ) correspond to a Josephson (particle conserving) and dissipative (particle non-conserving), respectively, coupling between the condensates.Calculating the integral (S7) numerically using the single-cell wavefunction found from 2DGPE simulations we find that I a,b < 0 which means Im(J a,b ) < 0 and therefore the dissipative coupling favours anti-phase synchronisation between neighbouring condensates, see Supplementary Note 7.This is in agreement with observations in experiment Fig. 2E-H and 2DGPE simulations Fig. 4E-G.The negative value of the integral intuitively makes sense because anti-phase displaced vortices constructively interfere between the two cells.Thus anti-phase condensates overlap more strongly with the pumped region between the traps and are populated more efficiently.
The lowest threshold solution of the coupled system (10)   48 .Such anisotropic low-threshold solution intuitively makes sense since polaritons flow stronger out of the trap parallel to the condensate dipole axis, thus creating optimal coupling (overlap) conditions between two spatially interacting condensates (e.g., similar to σbonding between atoms).
At higher pump powers these dipoles in each condensate loose stability and bifurcate each into either l = ±1 vortices.When |I b | > |I a | the system favours AFM ordering, like reported in Fig. 2E-H Note that I a,b < 0 which means that if the coupling anisotropy I b /I a is large then the critical power of the AFM vortex formation reduces.This means that by tailoring the pump P tot (r) one can adjust I a,b to shift the position of the transition point.
We obtain good results shown in Fig. 6B using I a,b as fitting parameters and solving the variational Gross-Pitaevskii equation for a single cell (6) and two cells (10)  as a function of increasing power.Our results are in excellent agreement with 2DGPE simulations previously shown in Fig. 4A, and thus explain the essential observations of our experiment.The blue stars indicate the onset of vorticity in the system.Going further, we have also extended Eq. (10)    realizations) we classified the state as dipole-dipole-vortex (bottom panels in Fig. S6).We believe this state can be explained by some tiny imperfections of the pump or sample disorder given the nonlinear nature of exciton-polariton system.The rest of the families are in a good agreement in our theoretical prediction.The numbers, in yellow, shown in Fig. 3B-D of the manuscript are renormalized such as to show statistical occurrence out of 93 single-shot realizations.

Supplementary Note 6: Variational Gross-Pitaevskii model for coupled polariton vortices
Extending Eq. ( 10) in the main text to arbitrary geometries of coupled nearest-neighbour polariton vortices gives, Here, ψ n,± is the phase and amplitude of the nth condensate component with OAM l = ±1, J a,b ∈ C are the tunneling rates between co-rotating and counter-rotating vortices, corresponds to the repulsive polariton-polariton interactions, R represents a gain saturation mechanism in the adiabatic exciton-reservoir limit, and p is the combined non-resonant optical pumping rate and cavity losses.The sum runs over nearest neighbours and Θ n,m is the angle of the link between two condensates in separate traps [notice the double winding in the exponent of Eq. (S1)] [47].In the special case of triangular geometry we have Θ n,m ∈ {0, 2π/3, 4π/3}.

3-cell system: AFM order in a single triangle
Here we numerically investigate the presence of AFM order in a single triangle when the power parameter p is scanned.It is convenient to characterise the behaviour of the three-cell system by defining a three-dimensional Bloch vector (or pseudospin) for each condensate similar to what is done in optics with light beams carrying OAM, where σx,y,z are the three Pauli matrices.Physically, projection on the S n x,y components means that the nth condensate has some dipolar structure.Projection on the S n z component means that the nth condensate has vorticity.Specifically, S n z > 0 corresponds to a counterclockwise vortex and S n z < 0 to a clockwise vortex.We numerically solve Eq. (S1) for the triangle of traps, averaging over 1000 random initial conditions.We show in Fig. S10A the normalised time-average and cell-average projection of the condensates on the dipole states ⟨S 2 x ⟩ + ⟨S 2 y ⟩ and vortex states ⟨S 2 z ⟩ where, ⟨S x,y,z ⟩ = From an optics perspective, these quantities are also similar to a light source's degree of linear polarization (DLP) and degree of circular polarization (DCP) except now for OAM.We also show in Fig. S10B the "amount" of AFM order using the following order parameter, Namely, if the order parameter ⟨M ⟩ < 0 then the system is preferentially AFM aligned.This quantity is useful when the condensate dynamics are nonstationary, and time-averaging over a long time window T is more meaningful.
The results in Fig. S10 show two distinct regions of interest.At low powers a single fixed point solution is dominant corresponding to dipole condensates arranged 120 • with respect to each other (see insets a-i and a-ii of example solution projected into the spatial domain for clarity), reminiscent of an XY ground state.In this regime there is no vorticity and therefore ⟨M ⟩ = 0.At higher powers we pass the vorticity threshold and vorticity starts growing monotonically.Here, a new family of attractors in which (on-average) two parallel vortices and one antiparallel appear (see insets a-iii and a-iv for example solution).For the given parameters |J a | < |J b | the vorticity threshold is associated with clear AFM order as can be seen in Fig. S10B.These results are similar to the observations on coupled nanodisk lasers (see Fig. 2 in Ref. [40]) Notice also that this solution lacks discrete rotational-or mirror-symmetry just like we observe in experiment (see e.g.Fig. 3D in main text).main manuscript, and Θ n,m is the angle of the link between two condensates in the plane.Notice that in the case of a triangular lattice the angles belong the set Θ n,m ∈ {0, 2π/3, 4π/3} and are responsible for the geometric frustration.
As we have discussed in the Materials and Methods section of the main manuscript around Eq. ( 11), the coupling rate between condensates J α,β n,m is complex which means that J is complex.Since the system selects an OAM configuration to optimize the gain then we are interested in maximizing Im(J ).From here on we are only concerned with the complex part of J and will drop the "Im(.)"notation.
The first inner product term in Eq. (S8) favours co-rotating vortices (FM arrangement).However, this term turns out to be weaker than the second term (see Fig. S11) which instead favours AFM arrangement.For simplicity, we will investigate the case of Θ n,m = 0 corresponding to a non-frustrated linear chain of condensates with uniform coupling strengths J +,− n,m = J b and J +,+ n,m = J a = 0. We then have, When pumped strongly enough, the presence of dipole states diminishes (i.e., above blue stars in Fig. 4 in the main manuscript) and we can parametrize the state vector as follows, s n = (cos (θ n ) cos (ϕ n ), cos (θ n ) sin (ϕ n ), sin (θ n ) cos (ϕ n ), sin (θ n ) sin (ϕ n )) T (S10) Here, θ n ∈ [0, π/2] determines the direction and amount of vorticity and ϕ n ∈ [0, 2π) the overall phase at site n.That is, θ n = 0 corresponds to OAM= +1 and θ n = π/2 to OAM= −1.Our gain functional then becomes, The cosine term is the XY-Hamiltonian previously pointed out in the work (37 ).The sine term is novel and importantly only depends on the angle θ n which determines the vorticity in each component.
If we project the vortex angles onto their nearest extremes θ n → θ n ∈ {0, π/2} we can rewrite our gain functional using a binary variable σ n ∈ {±} denoting the projected normalized OAM at each site, We remind that J b < 0 and optimizing (S12), i.e. optimizing the system gain, is the same as minimizing Eq. ( 1) in the main manuscript because they differ by a sign factor.Notice that the binarized term 1 − σ n σ m > 0 in the sum is always positive.Therefore, the optimal value is obtained when the "Ising" term is as positive as possible and the cosine term is as negative as possible.This happens for anti-phase θ n − θ m = π and AFM ordered vortices σ n σ m = −1, in agreement with our experimental and numerical results.

Figure 1 .
Figure 1.Driven trapped polariton condensates with OAM l = ±1 resembling classical spins.(A) Example of frustrated triangle for the Ising spins with AFM coupling.(B) Schematic of three polariton condensates in vortex states (yellow-violet surfaces), localized in the potential minima of the pump-induced energy landscape (black-grey surface).(C) Normalized experimentally measured timeaveraged real-space polariton condensate photoluminescence for a single cell.(D),(E) Examples of corresponding phase maps, revealing vortex and antivortex states, respectively.Dashed lines with circles schematically denote the excitation pattern.Numbers, in yellow, indicate statistical occurrence of states over one hundred single-shot realizations.

Figure 2 .
Figure 2. Vorticity onset and formation of vortex-antivortex pair in 2-cell structure.(A)-(D) Normalized experimentally measured time-averaged real-space polariton PL as a function of pump power.(A) Below and (B) at threshold the polariton PL is mostly coming from the ridge of the pump-induced traps with no clear vortex formation.(C) At higher power the condensate collapses inside the traps forming first dipole-dipole state, and then (D) a vortex-antivortex state.(E) Measured real-space polariton PL at P = 1.2P thr and (F)-(H) corresponding examples of phase maps.Dashed lines with circles schematically denote the excitation pattern.Numbers, in yellow, in (F)-(H) indicate statistical occurrence of states over one hundred single-shot realizations.

Figure 3 .
Figure 3. Vortex states in 3-cell structure.(A) Normalized experimentally measured time-averaged real-space polariton PL. (B)-(D) Corresponding extracted phase maps with numbers, in yellow, indicating statistical occurrence of states over one hundred single-shot realizations.Dashed lines with circles schematically denote the excitation pattern.

Figure 4 .
Figure 4. Vorticity dependence on pump power.(A) Calculated dependence of the stationary amplitude of the polariton wavefunction emerging in one, two, and three cells with varying pump power.The threshold (start of solid curves) lowers when more pump spots are introduced.The blue stars denote the onset of vorticity with increasing power.Below these stars the polaritons dominantly localize on top of the pump spots with no apparent vorticity.For the 3-cell triangle (green curve) an additional asymmetric frustrated state appears (green circles) at higher powers.(B)-(J) Example solutions from 2DGPE simulation of the condensate density and phase in building blocks of the lattice showing agreement with experiment.

Figure 5 .
Figure 5. Vorticity dependence on pump power.(A) Experimental time-integrated real-space polariton PL intensity (above condensation threshold) and (B) measured singleshot realization of the condensate phase illustrating the formation of a large-scale vortex lattice.(C) Calculated normalized mean number of AFM, FM, and free spins in the 22-spin AFM triangular Ising lattice as a function of energy.The shaded areas correspond to standard deviation (SD).The three types of spins are schematically depicted on the right.The average vortex populations obtained from experiment is plotted with the empty markers with error bars denoting SD.The obtained average populations from a discretized Gross-Pitaevskii simulations are shown with filled black, green and red markers and, remarkably, correspond to the minimum energy E = −42J, see Supplementary Note 6 for details on model and calculations.Overlaid honeycomb lattice in (A),(B),(D),(E) schematically shown with dashed lines and semi-transparent circles indicates the pump positions.Blue and red arrows show schematically the phase winding direction in each cell.(D) Time-averaged realization of the condensate density above threshold reveals donut-shaped PL. (E) Instantaneous condensate phase obtained from mean field modelling confirms vortices formation.

Figure 6 .
Figure 6.Vorticity onset in the single and pair of trapped condensate cells.(A) The colorscale depicts the Lyapunov potential (8) of the system with two minima (white dots) corresponding to vortex and antivortex solutions in which all phase space trajectories converge towards, underlining the deterministic nature of the system.Here we set p = R = 1 without loss of generality.(B) Reproduced results from Fig. 4A using the variational Gross-Pitaevskii model for a single (6) and a pair of coupled condensates(10).The blue stars denote the onset of vorticity for increasing power.
is then written ψ = (ψ 1,+ , ψ 1,− , ψ 2,+ , ψ 2,− ) T = √ N (1, 1, −1, −1) T e −iωt with occupation N = [p + |Im(J a + J b )|]/3 R and frequency ω = 3N +|Re(J a + J b )|.This is a antiphase dipole-dipole orientated head-totails (∞) ↔ (∞) with a condensation threshold of pthr = −|Im(J a + J b )|, in agreement with previous predictions . If |I b | < |I a | the converse is true and FM ordering is dominant (which we show experimental evidence of in the Supplementary Material).This vortex transition point P 2cells 0,vort can be analytically derived for the AFM case |I b | > |I a | (The FM case follows similar treatment) by locating the onset of linear instability for the dipole-dipole solution, to the case of arbitrarily coupled cells and analyzed the onset of vorticity and the preference towards AFM order when |I b | > |I a | in the 3-cell and 22-cell configuration (see Supplementary Note 6).

Figure S4 .Figure S5 .
Figure S4.Polarization-resolved polariton photoluminescence in a single cell.(A) Experimentally measured map of the Stokes parameter S3 for the trapped condensate excited with CW circularly polarized laser emission.Corresponding degree of circular polarization (DCP) map in B confirms that the condensate is mostly circularly polarized.Dashed lines with circles schematically denote the excitation pattern.

Figure 7 VFigure S6 .
FigureS6shows different examples of experimentally measured real-space phase maps corresponding to predicted families of the states observed in 3-cell structure (shown in Fig.3of the manuscript).The numbers here correspond to statistical occurrence of the states over one hundred individual realizations.We note that only in 7% of cases(7

Figure S9 .
Figure S9.Simulated single-shot realizations of the condensate phase map for the 22-cell structure.Dashed lines with semi-transparent white circles schematically denote pump pattern.

Figure
FigureS9shows simulated instantaneous phase maps, corresponding to Fig.5Dof the manuscript.Similar to experimental observations in Fig.S8the vortices stochastically flip their topological charges from realization-torealization, yet maintaining the dominant AFM order across 22-cell structure.In simple words, the polaritons OAM in the given cell strongly depends on the neighbouring cells.

BFigure S10 .
Figure S10.Comparison of the average condensate population in dipole states (Sx,y) and vortex states (Sz).We numerically solve Eq. (S1) over 1000 random initial conditions for different values of the effective power parameter p. a Shows the time-average and cell-average magnitude of projection onto the Sx,y and Sz Bloch vector components.Left and right insets above (A) show examples states from simulation for low and high power respectively.(B) Corresponding average vortex "magnetism" order parameter.Positive values indicate vortex "ferromagnetism" and negative values vortex "antiferromagnetism".

Figure S11 .
Figure S11.Example superposition of counter-rotating spatially displaced vortex wavefunctions with in-phase (left upper panels) and anti-phase (right upper panels) configuration.Bottom panel shows the coupling integral corresponding to Eq. (11) in the main manuscript.