Polaritonic network as a paradigm for dynamics of coupled oscillators

TL;DR

This paper demonstrates that polariton condensate networks can emulate various coupled oscillator models, enabling the study of complex synchronization phenomena and spin Hamiltonian minimization in a controllable photonic system.

Contribution

It establishes that polariton networks can represent multiple oscillator models like Kuramoto and Stuart-Landau, offering a versatile platform for simulating complex dynamics.

Findings

Polariton networks can emulate Kuramoto, Sakaguchi-Kuramoto, Stuart-Landau, Lang-Kobayashi oscillators.

They can realize regimes like synchronization, spin glasses, and secondary oscillation synchronization.

Control over system parameters allows interaction of different oscillatory sublattices.

Abstract

Photonic and polaritonic lattices have been recently theoretically proposed and experimentally realised as many-body simulators due to the rich behaviors exhibited by such systems at the macroscale. We show that the networks of polariton condensates encapsulate a large variety of behaviours of systems of coupled oscillators. By eliminating spatial degrees of freedom in nonresonantly pumped polariton network, we establish that depending on the values of experimentally tunable parameters the networks of polariton condensates may represent Kuramoto, Sakaguchi-Kuramoto, Stuart-Landau, Lang-Kobayashi oscillators and beyond. The networks of polariton condensates are therefore capable of implementing various regimes acting as analogue spin Hamiltonian minimizers, producing complete and cluster synchronization, exotic spin glasses and large scale secondary synchronization of oscillations. We suggest that the recently implemented control of the system parameters for individual sites in polariton lattices allows addressing for the first time the interaction of sublattices that belong to different oscillatory classes.

Equations

$$i\frac{\partial \psi }{\partial t}$$

$$\frac{\partial n_R}{\partial t}$$

$$\dot{{\Psi }_i}$$

$$\dot{R_i}$$

$$\dot{{\Psi }_i}=i(g\xi {\gamma }_i-1)|{\Psi }_i{|}^2{\Psi }_i-\xi {\gamma }_i|{\Psi }_i{|}^2{\Psi }_i+(1-ig)[({\gamma }_i-1){\Psi }_i+\sum _{j\ne i}{J}_{ij}{\Psi }_j].$$

$$\frac{1}{2}{\dot{\rho }}_i(t)$$

$${\dot{\theta }}_i(t)$$

$${\gamma }_i^{\prime }(t)=ϵ[{\rho }_{\mathrm{th}}-{\rho }_i(t)],$$

$${\dot{\theta }}_i(t)=(g\xi {\gamma }_i-1){\rho }_{\mathrm{th}}-g({\gamma }_i-1)-\sum _{j;j\ne i}\widetilde{J_{ij}}\sin ({\theta }_{ij}+\alpha ).$$

Full text

For a long time, two pervasive topics of modern science – dynamics of coupled oscillators and simulations of many-body solid state systems – have barely crossed each other’s paths. Complex dynamic behaviour of networks of coupled oscillators arises in various scientific disciplines ranging from biology, physics, and chemistry to social and neural networks as well as in established and emerging technological applications. Such networks served as paradigmatic models for understanding the mechanism of various collective phenomena. For instance, the Kuramoto oscillators kuramoto ; kuramoto2 have been successfully used to represent, study or even predict a wide variety of pattern formation in spatiotemporal systems, such as biochemical systems, neural networks, convecting fluids, and laser arrays. The reason for such power of networks of coupled oscillators in describing vastly different systems lies in the underlying symmetries of the system: these are described by similar universal order parameter equations that share similar characteristics halperin . Such symmetries make it possible to divide systems into various universality classes that differ only by the nature of the dynamics kadanoff and allow one not only to draw similarities between very different physical systems but also predict the behaviour of the new systems that fall into previously known universality class KBuniversal .

Traditionally, at the other end of the spectrum of nonlinear dynamical studies lie the complex many-body solid-state systems that are often considered as powerful platforms for simulating various elaborate Hamiltonians. A number of systems were realised implementing lattices of various physical origin: neutral atoms, ions, electrons in semiconductors, polar molecules, superconducting circuits, nuclear spins etc. QSreview . These are typically equilibrium systems that realise ground or excited states of their structure Hamiltonians. Recently, photonic and polaritonic lattices have emerged as promising platforms for many-body quantum and classical simulations photonlattices2007 ; bloch2008 . These systems are typically of a gain-dissipative nature, capable of symmetry breaking and spontaneous pattern forming, and have constant nonzero particle fluxes even at the steady state. Furthermore, as we argue in this article when the lattice elements have photonic component and the gain-dissipative nature the wavefunction packets evolve, interact and synchronise in a close resemblance to the coupled oscillators that are governed by the universal order parameter equations. As a result, on the one hand, many classical phenomena found in such lattices can be explained or predicted by the behaviour of the corresponding system of coupled oscillator networks from the same universality class; on the other hand, strong nonlinearities, spin polarization, sensitivity to magnetic fields, and individual site control greatly enrich possible states and dynamical regimes that can be generated in such lattices. In this article, we propose and theoretically justify the use of networks of exciton-polaritons (polaritonic networks) as a flexible universal platform to realise a vast array of known and extensively studied systems of coupled oscillators, to probe new exotic dynamical regimes and to create novel states of matter that may result from hybridisation of several different networks on one platform. In particular, we find the states that suggest possible practical implementations towards optical transistors and clarify the requirements for such a network to minimise classical spin Hamiltonians.

In the last decade, it emerged that strong light-matter interactions in semiconductor microcavities offer a versatile platform to realise nontrivial states, dynamics, localised structures. They consist of exciton-polaritons (polaritons) that are bosonic quasi-particles with a tiny effective mass which is typical $10^{-4}-10^{-5}$ of the bare exciton mass. Their energy-momentum dispersion curves can be controlled by appropriate detuning, and their properties and dynamics can be readily accessed by angular-resolved photo- or electroluminescence spectroscopy. A wealth of experimental results have been demonstrated with these systems, including Bose-Einstein condensates (BEC) bec , polariton lasers lasers , polariton parametric amplifiers ppa , and cavity quantum electrodynamics cavity . The polariton BEC or lasing has been demonstrated in various materials such as CdTe bec , GaAs gaas ; balili2007 , GaN gan , organic polymers organic and using optical pumping or electrically pumped exciton-polariton emitters yamamotoElPump .

We are interested in networks of $N$ polariton condensates created at lattice sites – the vertices of a two-dimensional graph at positions ${\bf r}_{i},i=1,\cdot\cdot\cdot,N.$ Many techniques are available to engineer a variety of the potential landscapes of polaritons chneiderRev ; amo . Polariton can be confined by strain-induced traps snoke , surface acoustic waves Cerda_PRB2012 , direct fabrication with the gold deposition technique Kim_NatPhys2011 , by using hybrid air gap microcavities Daffer , or by coupled mesas etched during the growth of the microcavity Winkler_NewJPhys2015 , by micropillars micropillars and in various geometries: square Kim_NatPhys2011 , triangular Kim_NewJPhys2013 , hexagonal Kusudo_PRB2013 , fully etched honeycomb Jacqmin_PRL2014 , Kagome Masumoto_NewJPhys2012 or even in quasi-periodic potentials amoFibb . However, the potential traps in a gain-dissipative system lead to complicated dynamics as the flow dynamics present in gain-dissipative systems even in the steady state is highly nontrivial in this geometry. To avoid such complications, polarion lattices can be optically engineered by exploiting the interactions between polaritons and reservoir excitons that can be injected in specific areas of the sample. Excitons barely move from the point where they are excited as they are orders of magnitude heavier than polaritons. Experimentally, the lattice is achieved by using a spatial light modulator that creates polariton condensates at the vertices of any prescribed graph wertz ; manni ; pendulum ; baumbergGeometrical ; BerloffNatMat2017 . This technique also allows controlling the intensities, $p_{i},$ of individual sites indexed by $i$ depending on the density of the polariton condensate at this site, if needed. In what follows we derive and discuss the network behaviour bearing in mind this technique of the lattice formation, however, other ways to create polariton lattices affect the parameters, but not the universality of the derived coupled oscillators equations.

The mean-field behaviour of polariton condensates is described by the generalised complex Ginzburg-Landau equation (cGLE) (often also referred to as a driven-dissipative Gross-Pitaevskii equation) coupled to the reservoir dynamics goveq ; carusotto ; reviewCarusotto . Although the process of Bose-Einstein condensation includes quantum effects, when condensate is formed it is accurately described by the mean-field equations as was shown in numerous experimental works meanfield ; meanfield2 ; meanfield3 ; meanfield4 ; meanfield5 ; meanfield6 ; meanfield7 ; meanfield8 ; meanfield9 . The equation on the wavefunction $\psi({\bf r},t)$ of the condensed system is coupled to the rate equation on the density of the hot reservoir $n_{R}({\bf r},t)$ so that $i\psi_{t}=-\frac{1}{2m}(1-i\hat{\eta}n_{R})\nabla^{2}\psi+U_{0}|\psi|^{2}\psi+g_{R}n_{R}\psi+\frac{i}{2}[R_{R}n_{R}-\gamma_{C}]\psi,$ and $n_{Rt}=-\left(\gamma_{R}+R_{R}|\psi|^{2}\right)n_{R}+P({\bf r},t),$ where we set $\hbar=1$, $U_{0}$ and $g_{R}$ are the polariton-polariton and polariton-exciton interaction strengths respectively, $\hat{\eta}$ is the energy relaxation wouters12 ; KBuniversal , $R_{R}$ is the rate of scattering from the hot reservoir into the condensates. The condensate ($\gamma_{C}$) and the reservoir ($\gamma_{R}$) relaxation rates describe photon losses in the cavity and hot exciton losses other than scattering into condensates. The incoherent pump source is described by the pumping intensity $P({\bf r},t)$. We nondimensionalise these equations by $\psi\rightarrow\sqrt{\gamma_{C}/2U_{0}}\psi$, $t\rightarrow 2t/\gamma_{C}$, ${\bf r}\rightarrow\sqrt{1/m\gamma_{C}}{\bf r}$, $n_{R}\rightarrow\gamma_{C}n_{R}/R_{R},P\rightarrow P\gamma_{C}^{2}/2R_{R}$ and introduce the dimensionless parameters $g=2g_{R}/R_{R}$, $b_{0}=2\gamma_{R}/\gamma_{C}$, $b_{1}=R_{R}/U_{0}$, $\eta=\hat{\eta}\gamma_{C}/R_{R}$. The resulting model yields

[TABLE]

A unique property of the exciton-polariton system is the flexibility with which the parameters $g,b_{0},b_{1},\eta$ can be controlled and changed to allow the system span various regimes bridging lasers or other nonequilibrium systems with equilibrium condensates and entering novel physical regimes. The lifetime of polaritons $\gamma_{C}$ is controlled by the accuracy of the cavity DBRs and spans two orders of magnitude bec ; nelsen . The detuning between the cavity photon energy and the exciton resonance determines the proportion of photon and exciton in the polariton and, therefore, the strength of the polariton-polariton and polariton-exciton interactions and effective mass reviewCarusotto . The repulsive interactions between excitons and polaritons $g_{R}$ can be further tuned by the pumping geometry, for instance, considering trapped condensates separated from the pumps baumbergTrapped .

The building block of our network is a single stationary condensate described by a wavefunction $\psi=\phi(r)$, created by a spatially localised radially symmetric incoherent pumping source $P=p(r)$. For instance, a Gaussian pump $p(r)=A\exp[-wr^{2}]$ where $w$ determines the inverse width has been widely used in experiments baumbergGeometrical ; BerloffNatMat2017 . In what follows we assume that the pumping intensity $A$ is chosen so that $\phi$ is normalised, so that $\int_{Q}|\phi|^{2}\,d{\bf r}=1$, where $Q$ is the entire plane of the cavity. We define the corresponding stationary reservoir profile as the steady state of Eq. (2), so that $n_{R}(r)=n(r)=p/(b_{0}+b_{1}|\phi|^{2})$. The networks of $N$ polariton condensates are created at lattice sites $i$ using a time and space varying pumping profile $P({\bf r},t)=\sum_{i=1}^{N}f_{i}(t)p(|{\bf r-r}_{i}|)$. The total wavefunction $\psi$ and the reservoir density $n_{R}$ can be approximated by $\psi({\bf r},t)\approx\sum_{i=1}^{N}a_{i}(t)\phi(|{\bf r-r}_{i}|)$ and $n_{R}({\bf r},t)\approx\sum_{i=1}^{N}k_{i}(t)n(|{\bf r-r}_{i}|)$ respectively. Such approximation is valid if the distance between the lattice sites exceeds the width of the condensate and reservoir GD-polariton . We use the shorth-hand notation $p_{i}\equiv p(|{\bf r-r}_{i}|)$ and similarly for $\phi$ and $n$. We eliminate the spatial degrees of freedom by multiplying Eq. (1) by $\phi_{i}^{*}$, Eq. (2) by $|\phi_{i}|^{2}$ and integrating both equations over the plane of cavity $Q.$ Previously, we have used this approach on a single Ginzburg-Landau equation that is relevant to the description of polariton condensates under several stringent assumptions: negligible blue-shift due to interactions of polaritons with the reservoir, short lifetime sample, fast reservoir relaxation, near threshold pumping intensity GD-polariton . Now we shall drop such restrictions and show how systems of different universality classes become relevant. We use the smallness of the overlap integrals for the wavefunctions of the different lattice sites polaritonGraph so that $l_{ij}\equiv\int_{Q}n_{i}\phi_{j}\phi_{i}^{*}\,d{\bf r}\gg\int_{Q}\phi_{j}\phi_{i}^{*}\,d{\bf r},$ $l=l_{ii}\gg\int_{Q}n_{i}|\phi_{j}|^{2}\,d{\bf r},$ and $H=\int_{Q}n\nabla^{2}\phi\phi^{*}\,d{\bf r}\gg\int_{Q}n_{i}\nabla^{2}\phi_{j}\phi_{i}^{*}\,d{\bf r}$ if $i\neq j.$ We also assume that for sufficiently smooth condensate profiles $\int_{Q}n_{i}\nabla^{2}\phi_{j}\phi_{i}^{*}\,d{\bf r}\ll l_{ij}.$ The dynamical equations on $\Psi_{i}(t)=a_{i}(t)\exp[-idt]$ and $R_{i}(t)=lk_{i}$ become

[TABLE]

where we used the notation $d=\int_{Q}\phi^{*}\nabla^{2}\phi\,d{\bf r},$ $h=\eta H/l,$ $\gamma_{i}=f_{i}\int_{Q}p|\phi|^{2}\,d{\bf r}/b_{0},\xi=b_{1}\int_{Q}n|\phi|^{4}\,d{\bf r}/lb_{0},$ $J_{ij}=(R_{i}l_{ij}+R_{j}l_{ji}^{*})/l.$ The energy relaxation parameter $\eta\ll 1$, therefore, $|H|<l$, so the term $|h|R_{i}\Psi_{i}$ will be neglected in comparison with $R_{i}\Psi_{i}$, whereas the imaginary part of $h$ will be assumed be absorbed by $g$. The coupling strength $J_{ij}$ is generally a complex number, so we write $J_{ij}\equiv{\cal J}_{ij}\exp[iv_{ij}]$ for real ${\cal J}_{ij}$ and $v_{ij}.$ In deriving Eqs. (3-4) we neglected higher order nonlinearities in $\Psi_{i}$ in the view of their smallness close to the condensation threshold. We consider several special cases of Eqs. (3-4).

Fast reservoir relaxation limit $b_{0}\gg 1$. In this limit, the reservoir dynamics can be replaced with its steady state, so $R_{i}=\gamma_{i}/(1+\xi|\Psi_{i}|^{2})\approx\gamma_{i}-\xi\gamma_{i}|\Psi_{i}|^{2}$ reducing the system of Eqs. (3-4) to the single equation

[TABLE]

For uniform pumping $\gamma_{i}=\gamma$ this is a celebrated Stuart-Landau system of coupled oscillators slref . This model can approximate a wide range of different oscillatory systems as it represents the normal form of an Andronov-Hopf-bifurcation. Operating close to an instability threshold lasers represent an example of the system close to such a bifurcation. We substitute $\Psi_{i}(t)=\sqrt{\rho_{i}(t)}\exp[{\rm i}\theta_{i}(t)]$ into Eq. (5) and separate real and imaginary parts to get

[TABLE]

where $\theta_{ij}=\theta_{i}-\theta_{j},$ $\tan\alpha=g$ and $\tilde{J_{ij}}={\cal J}_{ij}/\cos\alpha$. Note, that for the Gaussian pumping profile and wide reservoir, $|v_{ij}|\ll|{\cal J}_{ij}|$, so term $v_{ij}$ can be neglected since $l_{ij}\approx l_{ji}^{*}$ GD-polariton . For other network geometries such assumption may not be valid in which case we can absorb $v_{ij}$ into $\alpha_{ij}=\alpha-v_{ij}$.

Experimentally, the feedback can be applied to bring all the sites to the same density $\rho_{i}(t)=|\Psi_{i}(t)|^{2}=\rho_{\rm th}$ by combining Eq. (5) with an equation on the pumping adjustments

[TABLE]

where the parameter $\epsilon$ characterises the rate of such adjustment or its discrete version applied at discrete times $t_{n}$ (more appropriate for the current experimental control techniques), so that $\gamma_{i}(t_{n}<t\leq t_{n+1})=\gamma_{i}(t_{n})+\epsilon(t_{n+1}-t_{n})(\rho_{\rm th}-\rho_{i}(t_{n})).$ Under this control, close to the threshold $\rho_{i}\approx\rho_{\rm th}$ and Eqs. (6-7) reduce to a single equation

[TABLE]

This is the Sakaguchi-Kuramoto model of coupled oscillators sakaguchi with $\alpha$ representing a phase lag. Synchronisation and desynchronisation in this system has been extensively studied in the contexts as vastly different as a network of Wien-bridge oscillators in an experimental regime for which they can be approximated as phase oscillators wien , power grids consisting of many oscillating generators powergrid , and earthquake sequencing studies earth . The phase lag appears as a result of synaptic organisations in neuroscience systems, time delays in sensor networks, or transfer conductances in power networks. The Sakaguchi-Kuramoto model is a special case of the Winfree model with delta-function pulse shape $W_{1}(\theta)$ and a sinusoidal response curve $W_{2}(\theta)$, so that $\dot{\theta_{i}}=\omega_{i}+W_{2}(\theta_{i})\sum_{i=1}^{N}W_{1}(\theta_{j}).$ If the coupling is sufficiently weak and the oscillators are nearly identical, the phase can be replaced by its average over an entire period of oscillations leading to the Sakaguchi-Kuramoto model.

If $g=0$ ($\alpha=0$), then Eq. (9) reduces to the paradigmatic Kuramoto model kuramoto ; kuramoto2 that was a first tractable mathematical model for describing how coherent behaviour emerges in complex systems. This model exhibits a phase transition at a critical coupling, beyond which a collective behaviour is achieved. In our case, all natural frequencies are identical (and equal to $\rho_{\rm th}$) and the equation describes the negative gradient flow $\dot{\theta}=-\partial U(\theta)/\partial\theta$ for the smooth function $U(\theta)=-\sum_{i,j}\tilde{J_{ij}}\cos\theta_{ij}.$ Therefore, by LaSalle Invariance Principle (e.g. in khalil ) every trajectory converges to a minimum of the XY Hamiltonian $H_{XY}=-\sum_{i=1}^{N}\sum_{j=1}^{N}\tilde{J_{ij}}\cos(\theta_{i}-\theta_{j})$. In case of $\mathbb{S}^{1}$-synchronizing graphs all critical points are hyperbolic, so the synchronised state is the global minimum of $U(\theta)$, and all other critical points are local maxima or saddle points s1graphs . For arbitrary graphs the global minimum can be achieved by implementing the lowest pumping intensity that leads to threshold $\rho_{\rm th}$ GD-polariton . Polaritons graphs as global minimisers of the XY Hamiltonian has been theoretically justified and experimentally realised in our previous work BerloffNatMat2017 ; polaritonGraph . Using a resonant pumping in addition to nonresonant one (adding the terms proportional to $\Psi_{i}^{*}$ to the right-hand side of Eq. (1)) allows one to minimise the Ising or Potts Hamiltonians GD-resonant . Several other driving-dissipative platforms exploited this idea for minimisation of spin Hamiltonians: injection-locked laser systems yamamoto11 , the networks of optical parametric oscillators, yamamoto14 ; yamamoto16a ; yamamoto16b ; takeda18 , coupled lasers coupledlaser , and photon condensates KlaersNatPhotonics2017 . As the analogy with the coupled oscillators we presented above suggests, the minimisation of a spin Hamiltonian is realised for the type of couplings that allow for the flow to be represented by the negative gradient flow with real function $U$. This condition is satisfied only if the coupling matrix $J_{ij}$ in Eq. (5) is self-adjoint ${\cal J}_{ij}={\cal J}_{ji}$ and $v_{ij}=-v_{ji}$. For instance, if the couplings in Eq. (5) are of a pure Josephson type (e.g. $\dot{\Psi_{i}}=\cdot\cdot\cdot+{\rm i}\sum_{j}K_{ij}\Psi_{j}$ with real couplings $K_{ij}$) or have a non-negligible $g$, such a network will not necessarily minimise spin or any other Hamiltonian.

In addition, the parameter $g$ has a destabilising effect on the fixed points of Eq. (9). Different $\gamma_{i}$ that have to be maintained to allow all densities to reach the same value $\rho_{\rm th},$ provide each lattice element with its own “natural frequency,” $\omega_{i}=(g\xi\gamma_{i}-1)\rho_{\rm th}-g(\gamma_{i}-1),$ and, therefore, favour desynchronisation. In the network described by Eq. (9) synchronisation occurs when the coupling dominates the dissimilarity introduced by natural frequencies and the phase lag. The smaller is $g$ the more likely the global synchronisation is achieved. Concise results for complex networks are known for specific topologies such as, for instance, complete graphs, highly symmetric ring or linear graphs, acyclic graphs, and complete bipartite graphs with uniform weights. Also, Sakaguchi phase lag parameter $\alpha$ contributes to desynchronisation as it provides attraction and repulsion between the oscillator phases similar to the coupling time delay. The dependence of synchronisation and desynchronisation in polariton condensates on $g$ has been noted experimentally, but the reasons have not been previously identified wertz ; manni ; pendulum ; ostrovskaya2018 . Such a behaviour, however, is easily explained from the point of view of the dynamics of coupled oscillators.

Slow reservoir relaxation limit. Eq. (5) describes the direct coupling scheme: pumping at the mean field, calculated algebraically from the states of all oscillators, enters the coupling. The coupling scheme of Eqs. (3-4) is more complex: the mean field acts on the reservoir densities that obey its own nonlinear differential equations, and the acting force is a function of the reservoir state. This is similar to the famous example of synchrony on London’s Millennium Bridge where equations for the swinging mode of the bridge are coupled to the equations on individual pedestrians bridge , or to electronic or electrochemical oscillators that are coupled through the common macroscopic current or voltage, which obeys macroscopic equations describing the coupling circuit electric . By tuning the photonic component of polaritons one can change the polariton-polariton interactions up to 4 orders of magnitude ostrovskaya2018 which allows one to neglect the term $|\Psi_{i}|^{2}\Psi_{i}$, so that Eqs. (3-4) become similar to the Lang-Kobayashi equations (with $\Psi_{i}$ replaced by the electric field and $R_{i}$ by the population inversion of the $i-$th laser) obtained using Lamb’s semiclassical laser theory and capable of describing the dynamical behavior of coupled lasers lang80 ; acebron . We summarise all the regimes and models described above schematically in Figure 1.

The flexibility to tune the system parameters, the shape and geometry of the polariton lattice BerloffNatMat2017 , existence and tunability of the nonlocal couplings beyond the next neighbour interactions exotic , strong self-interactions of polariton condensates allow one to not only recreate the intriguing patterns, states and structures that fascinated the nonlinear dynamics community in the last couple of decades but also to enter novel regimes.

To illustrate this we consider two polaritonic networks configurations where the dynamic behaviour is the result of the interplay between two or more subsystems that belong to different types of coupled oscillators.

Triangular lattice in fast reservoir relaxation regime, low $g$, short lifetime, constant pumping. Any finite two-dimensional polariton lattice with the same pumping intensity across all the sites gives rise to an example of an interplay between Kuramoto (for $g=0$) or Sakaguchi-Kuramoto (for $g\neq 0$) oscillators away from the network boundaries with the Stuart-Landau oscillators at the sites close to the boundary. When all the lattice sites are equally pumped, the internal sites have the same number of neighbours and, therefore, receive equal particle fluxes from them. On the boundary, however, the number of neighbours is diminished, so fewer fluxes bring fewer particles decreasing the density of the oscillators. The difference in the “natural frequency” (given by the first term of the right-hand side of Eq. (7)) between the boundary and bulk condensates tends to desynchronise the lattice, but for lower pumping intensities this effect is overwhelmed by the synchronisation effect of the couplings, so the lattice is frequency synchronised (see Fig. 2(b)). As the pumping intensity is increased further above the threshold, the frequencies may still be synchronised, but the phases corresponding to a spin glass emerge as shown in Fig. 2(c,d). The appearance of the spin glass is due to nonlocal coupling and an increasing role of Dzyaloshinskii-Moriya interactions (DMIs), since the density misbalance between the sites makes the interaction directional and, therefore, of the DMIs type ourDM ; SkyrmionsNatureNano2013 . For even higher pumping intensities the desynchronisation takes place (Fig. 2(e-g)) resulting in a peculiar quasi-one-dimensional form of large-scale collective density oscillations in the direction of the larger axis of symmetry of the lattice. The collective density oscillates between left and right parts of the lattice with the appearance (and disappearance) of $\pi-$ phase difference between two clusters, indicating temporal formation of the domain wall. Figure 2 illustrates these regimes using the full spatially resolved simulations of the polaritonic networks (Eqs. (1,2)) in the fast reservoir relaxation limit. These states are robust, sustain an external noise, exist for a wide range of parameters, and reached from random initial conditions. The collective density oscillations also appear when larger lattices are considered.

We note that the system behaviour as the pumping intensity increases is reminiscent of the Mott metal-insulator transition. The frequency of the collective oscillations quench towards the transition point (see the inset in Fig. 2(a)) after which the oscillation frequency increases as a linear function, repeating the linear trend of the original Mott transition. We can induce reversible changes between these two states of the system by increasing or decreasing pumping intensity. Such oscillations suggest an interesting application towards implementation of polariton transistors polaritonTransistor that operate in the gigahertz range.

Triangular lattice, in both fast and slow reservoir relaxation regimes, low $g$, short and long lifetime, constant pumping. Our second example concerns cluster synchronization – a particular synchronization phenomenon that requires that synchronization occurs in each group, but there is no synchronization among the different groups cluster . The importance of cluster synchronization has been found in various applications including biological sciences and communication engineering, and various control schemes were designed to drive the network to cluster synchronization cluster2 . In polaritonic networks, cluster synchronization presents itself even in the case of equally pumped lattices. To illustrate cluster formation, we consider regular triangular lattices and numerically integrate Eqs. (1,2) in two opposite limits: fast (Fig. 3(a)) and slow (Fig. 3(b)) reservoir relaxation times. During a fraction of a nanosecond the system is frequency locked into several clusters, some of them are shown on Fig. 3(b,d) with dashed circles of the same color. The cluster state on Fig. 3(a,b) is a long-lived transient state that after a few nanoseconds evolves into chaotic oscillations. However, both reported states are stable to addition of the random noise and to intrinsic roughness of the sample modelled by adding noise potential to the right hand side of Eq. (1) and, therefore, may be detected experimentally.

The established links between polaritonic networks and other coupled oscillators systems suggest that these two examples are infinitely far from describing all possible dynamical regimes and hint that other intriguing states can be obtained. A particular interest is in finding Chimera states. A Chimera state is a spatio-temporal pattern in a network of identical coupled oscillators in which some of the oscillators synchronise (become frequency locked) while others remain incoherent (desynchronised) chimera1 ; stationarypattern . The simplest model that supports such states was reported as a pair of Sakaguchi-Kuramoto oscillator populations in which each oscillator has different coupling to others in the same group then to those in another group chimeraStrogatz . The type of the pattern depends strongly on the parameter $g$, and, therefore, the phase lag $\alpha$. The spiral pattern of chimeras appears when $\alpha$ is close to zero whereas the spot chimeras only appear when $\alpha$ is close to $\pi/2$ chimerasReview . The ‘turbulent chimeras’ were observed in some special cases of nonlocally coupled Stuart-Landau oscillators in which regions of local synchronisation appeared and vanished randomly over time turbulentChimeras . The Stuart-Landau system Eq. (5) has also been shown to produce a wealth of various stationary and dynamical behaviours, e.g. amplitude death, Hopf oscillations, large oscillations, quasi-periodicity, and chaos matthew91 . Polaritonic networks is a promising platform to study the wealth of these phenomena.

To summarize, we propose polaritonic networks as a paradigm for the dynamics of disparate systems of coupled oscillators. The dynamics of coupled oscillators often appear in the context of lasers and other driven-dissipative systems. It is therefore not surprising that systems of coupled oscillators are closely related to polaritonic systems. However, we show that depending on the system parameters and experimental controls used such networks can not only reproduce the behaviour of various known coupled oscillators and allow one to address intermediate regimes between different types of systems, but also to dramatically increase the physics of the system by combining different types of oscillators together in one interacting platform. This opens intriguing possibilities for entering novel hybrid regimes that have never been implemented before. Polaritonic network is a flexible and robust platform for achieving this.