A driven-dissipative spin chain model based on exciton-polariton condensates

TL;DR

This paper presents an analytical and numerical study of a driven-dissipative spin chain model based on exciton-polariton condensates, revealing controllable magnetic orderings and novel paired-spin states through bifurcation and coupling manipulation.

Contribution

It introduces a new driven-dissipative spin chain model based on exciton-polariton condensates, analyzing magnetic order formation and phase-controlled spin states.

Findings

Spontaneous spin bifurcation leads to magnetized states.

Multiple magnetic orders can be controlled via adiabatic pumping.

A novel paired-spin ordered state is demonstrated.

Abstract

An infinite chain of driven-dissipative condensate spins with uniform nearest-neighbor coherent coupling is solved analytically and investigated numerically. Above a critical occupation threshold the condensates undergo spontaneous spin bifurcation (becoming magnetized) forming a binary chain of spin-up or spin-down states. Minimization of the bifurcation threshold determines the magnetic order as a function of the coupling strength. This allows control of multiple magnetic orders via adiabatic (slow ramping of) pumping. In addition to ferromagnetic and anti-ferromagnetic ordered states we show the formation of a paired-spin ordered state $\left|\dots \uparrow \uparrow \downarrow \downarrow \dots \right. \rangle$ as a consequence of the phase degree of freedom between condensates.

Figures

Equations

$$i{\dot{\Psi }}_n$$

$$+\frac{1}{2}(\bar{\alpha }{S}_n+\alpha S_{zn}{\hat{\sigma }}_z){\Psi }_n-\frac{J}{2}\sum _{⟨nm⟩}{\Psi }_m,$$

$$S_n\equiv \frac{|{\Psi }_{n+}{|}^2+|{\Psi }_{n-}{|}^2}{2},$$

$$S_{zn}\equiv \frac{|{\Psi }_{n+}{|}^2-|{\Psi }_{n-}{|}^2}{2},$$

$$W_{\text{bif}}=W_{\text{lin}}+\eta \frac{{ϵ}^2+{\gamma }^2}{\alpha ϵ}.$$

$${\Psi }_{n+1}=e^{i{\phi }_{n+1}}{\Psi }_n,\text{(FM bonds)},$$

$${\Psi }_{n+1}=e^{i{\phi }_{n+1}}{\hat{\sigma }}_x{\Psi }_n,\text{(AFM bonds)},$$

$$i{\dot{\Psi }}_n$$

$$+\frac{1}{2}(\bar{\alpha }{S}_n+\alpha S_{zn}{\hat{\sigma }}_z){\Psi }_n.$$

$${ϵ}_J=ϵ+J$$

$${\omega }_J=-J($$

$${ϵ}_J^{\text{AFM}}=ϵ+2J\cos (2\pi m/N),{\omega }_J^{\text{AFM}}=0,$$

$${ϵ}_J^{\text{FM}}=0,{\omega }_J^{\text{FM}}=-2J\cos (2\pi m/N).$$

$$\delta {\psi }_{\pm }^{(m)}=u_{\pm }^{(m)}{e}^{ikma+\lambda t}+v_{\pm }^{(m)\ast }{e}^{-ikma+{\lambda }^{\ast }t}.$$

$$\mathbf{M}=\left(\begin{array}{cccc}{\mathcal{M}}^{(1)}& {\mathcal{M}}_J& 0_{4,4}& {\mathcal{M}}_J\\ {\mathcal{M}}_J& {\mathcal{M}}^{(2)}& {\mathcal{M}}_J& 0_{4,4}\\ 0_{4,4}& {\mathcal{M}}_J& {\mathcal{M}}^{(3)}& {\mathcal{M}}_J\\ {\mathcal{M}}_J& 0_{4,4}& {\mathcal{M}}_J& {\mathcal{M}}^{(4)}\end{array}\right).$$

$${\mathcal{M}}^{(n)}={\mathcal{M}}_E^{(n)}+{\mathcal{M}}_W^{(n)}+{\mathcal{M}}_{\gamma }+M_{ϵ}+M_{{\alpha }_1}^{(n)}+M_{{\alpha }_2}^{(n)}.$$

$${\mathcal{M}}_E^{(n)}=\left(\begin{array}{cccc}-E^{(n)}& 0& 0& 0\\ 0& E^{(n)}& 0& 0\\ 0& 0& -E^{(n)}& 0\\ 0& 0& 0& E^{(n)}\end{array}\right).$$

$${\mathcal{M}}_W^{(n)}=-\frac{i\eta }{4}\left(\begin{array}{cccc}2|{\psi }_{+}^{(n)}{|}^2+|{\psi }_{-}^{(n)}{|}^2& ({\psi }_{+}^{(n)}{)}^2& {\psi }_{+}^{(n)}{\psi }_{-}^{(n)\ast }& {\psi }_{+}^{(n)}{\psi }_{-}^{(n)}\\ ({\psi }_{+}^{(n)\ast }{)}^2& 2|{\psi }_{+}^{(n)}{|}^2+|{\psi }_{-}^{(n)}{|}^2& {\psi }_{+}^{(n)\ast }{\psi }_{-}^{(n)\ast }& {\psi }_{+}^{(n)\ast }{\psi }_{-}^{(n)}\\ {\psi }_{+}^{(n)\ast }{\psi }_{-}^{(n)}& {\psi }_{+}^{(n)}{\psi }_{-}^{(n)}& 2|{\psi }_{-}^{(n)}{|}^2+|{\psi }_{+}^{(n)}{|}^2& ({\psi }_{-}^{(n)}{)}^2\\ {\psi }_{+}^{(n)\ast }{\psi }_{-}^{(n)\ast }& {\psi }_{+}^{(n)}{\psi }_{-}^{(n)\ast }& ({\psi }_{-}^{(n)\ast }{)}^2& 2|{\psi }_{-}^{(n)}{|}^2+|{\psi }_{+}^{(n)}{|}^2\end{array}\right)-\frac{i}{2}G\mathbb{I},$$

$${\mathcal{M}}_{\gamma }=\frac{i\gamma }{2}\left(\begin{array}{cccc}0& 0& -1& 0\\ 0& 0& 0& -1\\ -1& 0& 0& 0\\ 0& -1& 0& 0\end{array}\right),$$

$${\mathcal{M}}_{{\alpha }_1}^{(n)}$$

$${\mathcal{M}}_{{\alpha }_2}^{(n)}$$

$${\mathcal{M}}_J=\frac{\cos (ka)}{2}\left(\begin{array}{cccc}-J& 0& 0& 0\\ 0& J^{\ast }& 0& 0\\ 0& 0& -J& 0\\ 0& 0& 0& J^{\ast }\end{array}\right).$$

$${\omega }_J=-J(e^{i{\phi }_i}+e^{i{\phi }_j}),$$

$${ϵ}_J=ϵ+J(e^{i{\phi }_i}+e^{i{\phi }_j}),$$

$$e^{i{\phi }_1}+e^{-i{\phi }_4}=e^{i{\phi }_2}+e^{-i{\phi }_1}=e^{i{\phi }_3}+e^{-i{\phi }_2}=e^{i{\phi }_4}+e^{-i{\phi }_3}.$$

$$a_1=a_3,$$

$$a_2=a_4,$$

$$b_1=b_2=b_3=b_4.$$

$$\cos ({\phi }_1)=\mp \cos ({\phi }_2).$$

$$i\dot{\Psi }=\frac{1}{2}\left[-i\left(g(S)+\gamma {\sigma }_x\right)+(1-i\Lambda )\left(\bar{\alpha }S+\alpha S_z{\sigma }_z-ϵ{\sigma }_x-\frac{\hslash {\nabla }^2}{m^{\ast }}+g_{P}P(\mathbf{r})\right)\right]\Psi .$$

Full text

Many-body spin systems, both classical and quantum, have found applications in a number of fields of rising complexity. Their Hamiltonians (Ising, $XY$, Heisenberg, Sherrington-Kirkpatrick, etc.) have been used to study collective behaviors such as familiarity recognition in neural networks Hopfield and Tank (1986), hysteresis in DNA interactions Vtyurina et al. (2016), combinatorial optimization problems in logistics, patterning, and economics Stein and Newman (2013); Bouchaud (2013). Besides their wide application, controllable spin lattices also offer insight into physical problems such as frustration Kim et al. (2010); Diep (2013), spin-ice Nisoli et al. (2013); Drisko et al. (2017), spin-wave dynamics Katine et al. (2000); Glazov and Kavokin (2015), domain wall motion Stenger et al. (1998); Liew et al. (2008); Adrados et al. (2011), and spin-glass formation Mezard et al. (1986); Stein and Newman (2013). A driven-dissipative spin lattice, where both phase and spin of the vertices are free, has yet to be addressed. Here, in contrast to entropy and minimum energy principles (as in the Ising model), the stationary physics of the system is governed by the balance of gain and decay with remarkably different solutions Le Boité et al. (2013). Currently, only limited investigation has been devoted to driven-dissipative lattice systems where recent works have proposed “simulators” based on interacting exciton-polariton condensates Berloff et al. (2016) and Ising machines with degenerate optical parametric oscillators Inagaki et al. (2016).

Nonresonantly excited spinor exciton-polariton (or simply polariton) condensates Fraser et al. (2016); Sanvitto and Kena-Cohen (2016a) have developed into a popular platform for cutting edge opto-electronic and opto-spintronic technologies Liew et al. (2011); Sanvitto and Kena-Cohen (2016b). The driven-dissipative condensates are realized by matching the gain and the decay of polaritons through continuous external driving of either optical or electrical nature. These macroscopic coherent states possess a spin and a phase degree of freedom, strong nonlinearities and a small effective mass, allowing them to interact and synchronize with other spatially separate condensates over long distances (hundreds of microns) Nelsen et al. (2013), making them interesting candidates for driven-dissipative spin lattices. Recently it was reported that a spinor polariton condensate bifurcates at a critical pump intensity into either of two highly circularly polarized states using a continuous linearly polarized nonresonant excitation Ohadi et al. (2015). The emission polarization is explicitly related to the polariton condensate pseudospin orientation (from here on spin) Kavokin et al. (2004). The system has since then been extended to polariton condensate spin pairs Ohadi et al. (2016) which can controllably display alignment of antiferromagnetic (AFM) and ferromagnetic (FM) nature. This spin degree of freedom offers a unique way to study ordering amongst coupled spin vertices in various lattices.

In this paper, we extend such polariton condensates to an infinite chain model and present methods of controllably producing different spin-ordered chains. We solve exactly and numerically analyze the stationary states of the infinite chain of spin-bifurcated condensates with nearest-neighbor same-spin coupling in the tight binding approach. The stationary solutions correspond to ferromagnetic, antiferromagnetic, and paired-spin order states of two-up and two-down spins (P) (see Fig. 1). States characterized by FM bonds with zero phase-slip and AFM bonds with $\pi$ phase-slip are shown to have a minimum bifurcation threshold, and are stable against long-wavelength fluctuations. Monte-Carlo trials with adiabatic ramping of the pump intensity on a cyclic system of 4 condensates give a phase diagram in full agreement with the predicted minimum threshold winners as a function of coupling strength. This clear hierarchy for the probability of formation is an important prerequisite for a spin-lattice simulator. Non-adiabatic trials on the other hand result in a complex phase diagram, as a result of the initial condition progressing to its nearest phase space attractor. In addition to spatially uniform stationary states, we find that frustrated or defect states, with oscillating spinors can appear in this system.

Condensation of the bosonic spin $\pm 1$ quasiparticles known as exciton-polaritons Kasprzak et al. (2006); Balili et al. (2007); Lai et al. (2007); Carusotto and Ciuti (2013); Byrnes et al. (2014) is regarded as the solid-state analog of cold-atom Bose-Einstein condensates Bloch et al. (2008). Its spin structure and strong interactions allow one to realize spinor condensates where macroscopic coherence and superfluid character give birth to many intriguing phenomena such as the nonlinear optical spin Hall effect Kammann et al. (2012); Antón et al. (2015), the formation of polariton half-vortices Lagoudakis et al. (2009), and spontaneous symmetry breaking Ohadi et al. (2012). They offer a new path towards spin manipulation Žutić et al. (2004) with already promising results on spin-switches Amo et al. (2010); Dreismann et al. (2016), and transistors Ballarini et al. (2013).

Driven-dissipative polariton condensates can be accurately modeled using a coherent macroscopic spinor order parameter $\Psi=(\Psi_{+},\Psi_{-})^{\text{T}}$ where $\Psi_{\pm}$ are the spin-up and spin-down components respectively. Similarly to the pair of polariton condensates, where the transport of polaritons from one condensate to another can be regarded as a form of coherent coupling in the tight-binding approximation Ohadi et al. (2016), a system of many condensates labeled by index $n$ can be described by coupled dynamical equations

[TABLE]

[TABLE]

[TABLE]

where the sum is over nearest neighbors. Here we define $g(S_{n})\equiv-W+\Gamma+\eta S_{n}$ as the pumping-dissipation imbalance, $\Gamma$ is the (average) dissipation rate, $W$ is the incoherent in-scattering (or pump rate), and $\eta$ defines the gain-saturation nonlinearity Keeling and Berloff (2008). The birefringence of the system corresponds to the splitting of the $XY$-polarized states in both energy ($\epsilon$) and decay-rate ($\gamma$). The interaction parameters are written as $\bar{\alpha}=\alpha_{1}+\alpha_{2}$, $\alpha=\alpha_{1}-\alpha_{2}$, where $\alpha_{1}$ and $\alpha_{2}$ are the same-spin and opposite-spin polariton-polariton interaction constants, respectively. Finally, $\Re{(J)}>0$ depicts the strength of the coherent same-spin Josephson type coupling between the condensates, whereas $\Im{(J)}$ gives the dissipative coupling between the condensates Aleiner et al. (2012). In particular, $\Im{(J)}<0$ defines the inter-site damping to account for energy relaxation. Eq. (2) is the average condensate population and Eq. (3) is the circular polarization intensity (considered as a spin here).

The critical pump intensity for condensation in a single condensate is determined by the lowest decay rate mode and can be written, $W_{\text{lin}}=\Gamma-\gamma$, resulting in a linearly polarized emission. At higher pump intensities the order parameter bifurcates into either a spin-up or spin-down state due to instability in the linearly polarized modes due to their splitting ($\epsilon+i\gamma$) and polariton-polariton interactions. For a single condensate, the critical bifurcation threshold is Ohadi et al. (2015)

[TABLE]

In the following, we work above this critical pump threshold such that each condensate is either in a ‘spin-up’ or a ‘spin-down’ state.

Formally, for identical lattice sites, the symmetry-conserving stationary solutions can be found by using the following spinor ansatz:

[TABLE]

[TABLE]

where ($\varphi_{n+1}$) is the phase shift moving from the condensate in question to its nearest neighbor $n+1$, and is to be determined. Eq. 1 can now be written as:

[TABLE]

This corresponds to a single condensate with complex renormalized splitting $\epsilon_{J}$ and energy shift $\omega_{J}$, arising from AFM bonds and FM bonds, respectively. The strength of these parameters depends on the relative phases between nearest-neighbors in the system. The modified parameters of a chain with two nearest neighbors can be written as

[TABLE]

where $\delta_{i}=1,0$ for AFM or FM bonding respectively.

The problem of FM- and AFM bonded condensates has thus been reduced to a single condensate with a known solution Ohadi et al. (2015). The requirement for site-independent $\epsilon_{J}$ and $\omega_{J}$, and cyclic boundary condition of integer $2\pi$ for the phase accumulated around the closed chain, restricts the possible phases of the bonds. For chain systems of either FM or AFM ordering it can be shown that the coupling results in (see Sec. E):

[TABLE]

where $m=0,1,2,\dots$ and $N$ is the number of condensates in the chain. For AFM chains, $N$ must be an even number since the spin unit-cell is $\left|\uparrow\downarrow\right.\rangle$. In addition, there is a paired-spin (P) state, where each site has one FM and one AFM bond, and the spin unit-cell is $\left|\uparrow\uparrow\downarrow\downarrow\right.\rangle$. P solutions have $\omega_{J}^{\text{P}}=\pm J$ and $\epsilon_{J}^{\text{P}}=\epsilon\pm J$, where the signs are independent. Due to the periodicity of the stationary solutions, the essential physics of the spin-bifurcated condensate chain system can be captured within a chain of 4 condensates characterized by 10 distinct solutions (see Sec. A). Furthermore, we confirm the analogy between the solutions of the tight binding model (Eq. 1) to a 4 condensate chain accounting for the $(x,y)$ spatial degrees of freedom (see Sec. F).

To identify the stable stationary solutions, we perform a long wavelength stability analysis (see Sec. B) for the set of coupled equations describing linear fluctuations along the periodic 4 condensate chain. Three lowest bifurcation threshold solutions of FM, AFM, and P spin order are found to have negative real-part Lyapunov exponents $\lambda$ within the first Brillouin zone of the chain (see Fig. 2), and hence are completely stable against fluctuations travelling along the chain. These solutions are characterized by a 0 phase slip between FM bonded condensates and $\pi$ phase slip between AFM bonded condensates (see Sec. D for parameter values).

To calculate the phase diagram and verify the analysis, we perform Monte-Carlo simulations of the periodic 4 condensate chain as a function of pump intensity $W$ and coupling strength $J$. Fig. 3 shows the result of 10 Monte-Carlo (MC) trials at each site over a $100\times 100$ pixel map in parameter space. For each realization of the numerical experiment, the pump intensity is linearly increased from $W_{0}=0.7W_{\text{bif}}$ to $W$ at a rate $10^{-4}\times W_{\text{bif}}$ ps*-1*, similar to the ramp-times achieved in experiments Ohadi et al. (2015). Three distinct phases corresponding to the stable AFM, P, and FM stationary spin patterns are observed, as identified by the long-wavelength stability analysis.

To explain the regimes of each state (red areas in Fig. 3), we plot the spin bifurcation threshold power $W_{\text{bif}}$ against $J$ in Fig. 4. The thresholds are calculated using Eq. 4, with $\epsilon\rightarrow\Re{(\epsilon_{J})}$, $\gamma\rightarrow\gamma+\Im{(\epsilon_{J})}$, $\Gamma\rightarrow\Gamma+\Im{(\omega_{J})}$. As the pump power is slowly increased, the state that reaches the bifurcation threshold first wins, since it has time to stabilize before competing states can bifurcate. The calculated phase-boundaries of $J/\epsilon=0.42,0.91$ for the AFM-P and P-FM boundaries are in close agreement with the Monte-Carlo simulations of Fig. 3. The asymptotic behavior in Fig. 4 for the AFM and P states at $J/\epsilon=0.5$ and $J/\epsilon=1$ respectively, is a consequence of $\epsilon_{J}$ approaching zero, destabilizing the stationary solution. We note that the ramp-time of the pump can influence the phase diagram. Fast ramp times soften the competitive advantage of a low spin-bifurcation threshold, resulting in a blurring of the phase-boundaries (see Sec. C). The depression at low $W$ in Fig. 3(c) corresponds to an area of instability outlined by the black dashed line calculated using linear stability analysis for fluctuations at $k=0$ (see Sec. B).

We note that the MC iterations do not always result in a stationary AFM, P, or FM steady-state described above. Nonstationary symmetry-breaking solutions can arise close to stability boundaries due to the finite ramp rate of the pump. The analysis of highly nontrivial evolutions of the system in this case is beyond the scope of this work.

In addition to the stationary states, an oscillating limit cycle solution composed of three spins against one opposite spin is often observed for low $J$ and high $W$, as shown in Fig. 3(d) and Fig. 5(a). A time-trace of the $S_{z}$ spin component of this state is plotted in Fig. 5(b). Though the average spin on each site has converged, the spin precession indicates a superposition of states that are phase locked. Interestingly, the energy of the spinor-components of the limit cycle state correspond to that of P state, $\omega_{J}=-J$ and $\epsilon_{J}=\epsilon-J$, except for the minority spin population in the opposing condensate (e.g., $\Psi_{+4}$ polaritons from Fig. 5) which also populates a separate peak in energy. Thus the limit cycle solution can be characterized as a ‘frustrated P-state’, described by multiple energies $\omega_{J}$ and splittings $\epsilon_{J}$, and resulting in an oscillating spinor and frequency comb emission, similar to discussed in Ref. Rayanov et al. (2015). In larger chains, the limit cycle states can appear as a result of inhomogeneity in the chain couplings.

In conclusion, we have solved analytically and investigated numerically solutions in an infinite chain of coupled driven-dissipative spinor polariton condensates. A mixture of intra-spin coupling and nearest-neighbor inter-coupling allows not only controllable formation of antiferromagnetic states or ferromagnetic states, but also shows solutions with mixed antiferromagnetic and ferromagnetic bonding. We find that minimum bifurcation threshold determines the spin order in the chain. The one-to-one correspondence between the spin and the phase-slips of the lowest threshold states makes this system binary and opens the possibility of mapping it to binary models such as the 1D Ising Hamiltonian where the minimization of loss (bifurcation threshold) replaces minimization of energy. Our work is an important step towards understanding and controlling spin order in open-dissipative nonlinear spin lattices.

Acknowledgements.—This work was supported by the Research Fund of the University of Iceland, The Icelandic Research Fund, Grant No. 163082-051, grant EPSRC EP/L027151/1, the Mexican Conacyt Grant No. 251808, and the Singaporean MOE grants 2015-T2-1-055 and 2016-T1-1-084. I.A.S. acknowledges support from a mega-grant № 14.Y26.31.0015 and GOSZADANIE № 3.2614.2017/ПЧ of the Ministry of Education and Science of Russian Federation

Appendix A Complete solutions of 4 condensate chain system

It has been established that stationary chain solutions of either FM or AFM ordering results in a modified single-condensate dynamical equation (Eq. 7) with only shifted parameters according to Eqs. 10-11. Stationary chain solutions of mixed FM and AFM bonding (P solution) follow the same procedure but with only $\pi$ or zero phase slips possible between condensates. Focusing on a chain of 4 condensates (smallest cell to encompass all spin orderings periodically) we find 10 distinct solutions which are summarized in Fig. 6. As the number of condensates in the chain increases, more solutions of FM or AFM ordering become available but the number of P solutions remains fixed. It’s worth mentioning that panel (c) and (f) are special cases where the phases between neighboring condensates result in a cancellation such that $\omega_{J}=0$ and $\epsilon_{J}=\epsilon$ in Eqs. 8-9.

Appendix B Linear stability analysis

In this section we formulate the linear stability analysis for a periodic solution for an infinite chain of condensates. This solution is constructed by periodic repetition of a particular solution for the closed ring of four condensate, and it has, in general, the period $4a$, where $a$ is the nearest-neighbor distance. The perturbed solution can be written as $\psi_{\pm}^{(m)}+\delta\psi_{\pm}^{(m)}$, where $m$ is the number of the condensate. The unperturbed solution is periodic, $\psi_{\pm}^{(m)}=\psi_{\pm}^{(m+4)}$, and the perturbation is chosen in the form of plane wave ($\hbar=1$)

[TABLE]

Here the complex amplitudes are also set to be periodic, $u_{\pm}^{(m)}=u_{\pm}^{(m+4)}$ and $v_{\pm}^{(m)}=v_{\pm}^{(m+4)}$. There are 16 linearized equations for the amplitudes and 16 Lyapunov exponents $\lambda(k)$. The solution is stable when all of them satisfy $\Re\{\lambda(k)\}\leq 0$.

The linearized equations for the vector $\mathbf{U}=\{u_{+}^{(1)},v_{+}^{(1)},u_{-}^{(1)},v_{-}^{(1)},\dots,u_{-}^{(4)},v_{-}^{(4)}\}^{\mathrm{T}}$ can be written in matrix form as $i\lambda\mathbf{U}=\mathbf{M}\cdot\mathbf{U}$, where the $16\times 16$ matrix $\mathbf{M}$ can be presented in the $4\times 4$ blocks form:

[TABLE]

Here, $\mathcal{M}^{(n)}$ is the matrix describing the fluctuations in the $n$-th condensate within the elementary cell, $n=1,2,3,4$, $0_{4,4}$ is the $4\times 4$ zero matrix, and $\mathcal{M}_{J}$ is the same-spin coupling between nearest neighboring condensates.

For the matrices in (13) we have

[TABLE]

The first matrix is defined by the energy of $n$-th condensate

[TABLE]

The second matrix arises from the harvest and saturation rates of the condensate from the static reservoir

[TABLE]

where $G=\Gamma-W$ and $\mathbb{I}$ is the identity matrix. The third and the fourth matrices arise from coupling between up and down components

[TABLE]

The fifth and the sixth matrices arise from the interactions:

[TABLE]

The final matrix describes the same-spin coupling between the nearest-neighboring condensates

[TABLE]

Appendix C Non-adiabatic Monte-Carlo for undamped 4 condensate chain

Here we give results analogous to Fig. 3 but with damping absent ($\Im{(J)}=0$) and instantaneous switching of the pump intensity at its mark value. Unlike Fig. 3 where states with the lowest bifurcation threshold were dominant, we uncover a more complex probability map in Fig. 7 through 30 MC trials bounded by their $k=0$ stability regions (black dashed lines) predicted by Eq. 13.

Fig. 7 shows 77% of the data divided between 6 solutions of the 4 condensate chain. The remaining 4 solutions from Fig. 6 are not observed since they are unstable over the entire $J-W$ map. Another 11% of the data (not shown here) ended in oscillating limit cycle states discussed in Fig. 5. The remaining 12% were categorized as nonstationary/chaotic.

The complex features of the probability maps in Fig. 7 as opposed to the more simplistic ones in Fig. 3 highlight the important role of damping in the system and adiabatic switching of the pump intensity. Intuitively, the complicated features in Fig. 7 arise from the order parameter overshooting many possible stable minima in the phase space of the system. It then becomes a matter of the nearest and strongest attractor to stabilize the solution.

Appendix D Numerical methods

A QR algorithm is implemented to solve the eigenvalue problem of Eq. 13. Eq. 1 is solved using a variable-order Adams-Bashforth-Moulton predictor-corrector method. The parameters used for 0D simulations were: $\eta=0.02$ ps*-1*; $\Gamma=0.1$ ps*-1*; $\epsilon=0.04$ ps*-1*; $\gamma=0.2\epsilon$; $\alpha_{1}=0.01$ ps*-1*; $\alpha_{2}=-0.5\alpha_{1}$.

The parameters used for 2D simulations of Eq. 27 were: $m^{*}=5\times 10^{-5}m_{0}$; $\eta=0.01$ ps*-1*; $\Gamma=0.2$ ps*-1*; $\epsilon=0.015$ ps*-1*; $\gamma=0.2\epsilon$; $\alpha_{1}=0.003$ ps*-1*; $\alpha_{2}=-0.5\alpha_{1}$; $d=12$ $\mu$m; $\sigma=10.3$ $\mu$m. Where $m_{0}$ is the free electron rest mass.

Appendix E Derivation of the coupling contribution in chain systems

Consider the stationary condensate chain composed entirely of either FM- or AFM bonds. According to Eqs. 8-9, each condensate with two nearest neighbors is presented with a term

[TABLE]

for two FM bonds or

[TABLE]

for two AFM bonds. Here $\varphi_{i,j}$ is the phase difference of moving from the condensate in question to its neighbor. This contribution can in general be a complex number appearing equally in each condensate. In this section we show that this number must stay real for a chain system. The lattice unit cell of the chain system is one condensate and one bond. Assuming that the chain closes on itself, the number of free variables ($\varphi_{i}$) is then equal to the number of independent equations. The following can then be generalized to any number of condensates in a chain.

Let’s now imagine 4 condensates locked in a chain. We can classify the phase jumps going clockwise as $\{\varphi_{1},\varphi_{2},\varphi_{3},\varphi_{4}\}$ (see Fig. 8). It is obvious that $\omega_{J}$ and $\epsilon_{J}$ must be equal for all condensates in the chain in order to have a steady state. Thus the phase contribution $e^{i\varphi_{i}}+e^{i\varphi_{j}}$ must be the same for all condensates. This means that the 4 stationary condensates allow us to write,

[TABLE]

Writing $e^{i\varphi_{i}}=a_{i}+ib_{i}$ where $a_{i},b_{i}\in\mathbb{R}$ it’s then easy to show that,

[TABLE]

The real part of the contribution $e^{i\varphi_{i}}+e^{i\varphi_{j}}$ is thus equal for each condensate but the imaginary part gets canceled. Consequently, from $|e^{i\varphi_{i}}|^{2}=1$ we come to the solution $a_{1}=\mp a_{2}$, which can more clearly be written:

[TABLE]

The minus sign in Eq. 26 corresponds to a cancellation in the coupling with no shift in $\omega_{J}$ or $\epsilon_{J}$ whereas the plus sign mandates the opposite. Applying the constraint $e^{i(\varphi_{1}+\varphi_{2}+\varphi_{3}+\varphi_{4})}=1$ corresponding to a full cycle in our 4 condensate chain we come to the conclusion that the only possible values of coupling in the latter case are $\cos{(\varphi_{1})}=\cos{(2\pi m/N)}$ where $m=0,1,2,\dots$ and $N$ is the number of condensates in the chain. As the number of condensates increases in the chain, more solutions become available.

The same procedure can be applied to a state where each condensate has one FM- and one AFM bond (P solutions). Then only $e^{i\varphi_{i}}=\pm 1$ satisfies the chain.

Appendix F 4 condensate chain with spatial degrees of freedom

The tight binding model (Eq. 1) offers a simple solution to the stationary spin patterns in the condensate chain. We find that these exact solutions can also be produced with little difficulty accounting for the spatial degree of freedom. The complex Ginzburg-Landau equation can be written then Keeling and Berloff (2008); Wouters and Carusotto (2007); Wouters et al. (2010):

[TABLE]

Here, $m^{*}$ is the effective mass of the polaritons. The exciton reservoir is taken to be completely static and the induced repulsive potential is then given by an effective interaction constant $g_{P}$. The remainder of the parameters serve the same purpose here as in the tight binding model with $W=P(\mathbf{r})$. We note that modeling the system by coupling an exciton reservoir rate equation to the order parameter only requires rescaling of the parameters and does not critically affect the observed solutions in Eq. 27 when the decay rate of the reservoir is taken to be large compared to the polariton lifetime Smirnov et al. (2014).

The pump $P(\mathbf{r})$ is a $3\times 3$ square-arrangement of Gaussians separated by the lateral distance $d$ and with a FWHM $\sigma$ which then form four potential minimum in a $2\times 2$ arrangement where the polaritons condense. In Fig. 9 we show three solutions of different spin order in a closed 4 condensate chain. The AFM, P, and FM solutions possess phase slips corresponding to the lowest bifurcation threshold solutions from Fig. 6(b,d,g), and Fig. 3(a-c). Each solution can be achieved by either tuning the strength of polariton interaction with the pump $g_{P}$, or by increasing the strength of the center Gaussian pump spot causing an increased barrier between the condensates which effectively changes the coupling strength $J$.