# Time-Delay Polaritonics

J. D. Töpfer1,2∗ H. Sigurdsson1,2 L. Pickup2 P. G. Lagoudakis1,2∗ 1Skolkovo
Institute of Science and Technology Novaya St., 100, Skolkovo 143025, Russian
Federation  2Department of Physics and Astronomy, University of Southampton,
Southampton, SO17 1BJ, United Kingdom  ∗email<EMAIL_ADDRESS><EMAIL_ADDRESS>
(August 27, 2024)

Non-linearity and finite signal propagation speeds are omnipresent in nature,
technologies, and real-world problems, where efficient ways of describing and
predicting the effects of these elements are in high demand. Advances in
engineering condensed matter systems, such as lattices of trapped condensates,
have enabled studies on non-linear effects in many-body systems where exchange
of particles between lattice nodes is effectively instantaneous. Here, we
demonstrate a regime of macroscopic matter-wave systems, in which
ballistically expanding condensates of microcavity exciton-polaritons act as
picosecond, microscale non-linear oscillators subject to time-delayed
interaction. The ease of optical control and readout of polariton condensates
enables us to explore the phase space of two interacting condensates up to
macroscopic distances highlighting its potential in extended configurations.
We demonstrate deterministic tuning of the coupled-condensate system between
fixed point and limit cycle regimes, which is fully reproduced by time-delayed
coupled equations of motion similar to the Lang-Kobayashi equation.

## Introduction

Time-delay is widespread in nature and occurs when the constituents of a given
system interact via signals with a finite propagation time
mackey1977oscillation . If the characteristic timescale of the system, such as
the period of a simple pendulum, is much longer than the signal propagation
time, then one arrives at familiar examples such as Huygens clock
synchronisation described by instantaneous interactions oliveira2015huygens .
When propagation times are appreciably long, the role of the system’s history
is enhanced and the interactions are said to be time-delayed. Such systems
dictate the neurological function of our brains, affect traffic flow,
influence economic activities, define population dynamics of biological
species, regulate physiological systems, determine the stability of lasers,
and have application in control engineering erneux2009applied ;
atay2010complex . Besides their ubiquity in nature and science, coupled
systems with continuous time-delayed interactions exhibit interesting
mathematical properties such as an infinite dimensional state space, i.e. for
a fixed time-delay $\tau$ there are infinitely many initial conditions of the
system in the time interval $-\tau\leq t\leq 0$ needed to predict the dynamics
for $t>0$ schuster1989mutual . Time-delayed interaction or self-feedback is
known to greatly increase the dynamical complexity of a system, giving rise to
chaotic motion, chimera states Larger_PRL2013 , as well as being able to both
stabilise and destabilise fixed point and periodic orbit solutions
atay2010complex .

Technological applications of delay-coupled systems appear in diverse areas
such as control engineering scholl2008handbook , high speed random-bit
generation uchida2008fast , secure chaos communication argyris2005chaos , but
have also recently emerged in machine learning, where the demand for neuro-
inspired computing units has led to various hardware realisations of
artificial neural networks van2017advances . The intrinsic high-dimensional
state space of delay-coupled systems enables an efficient platform for tasks
such as pattern recognition, speech recognition, and time series prediction as
demonstrated in electronic Appeltant2011 , photonic larger2012photonic ;
Brunner_NatCom2013 ; Vandoorne2014 , and optoelectronic
paquot2012optoelectronic ; Martinenghi_PRL2012 ; Larger_PRX2017  systems.

In this work, we demonstrate the prospect of an ultrafast microscale platform
for engineering time-delayed coupled oscillator networks based on condensates
of microcavity exciton-polaritons (from here on polaritons)
kavokin_microcavities_2007 ; li_tunable_AQT2016 . Polaritons are bosonic
quasi-particles that can undergo a power-driven quantum phase transition to a
macroscopically occupied state with long-range phase coherence
Kasprzak_Nature2006 . This phase transition is associated with the formation
of a polariton condensate, a non-linear matter-wave quantum fluid, which
differs from classical cold atom Bose-Einstein condensates, and liquid light
droplets wu_cubic-quintic_2013 ; zhang_particlelike_2019  due to its non-
equilibrium nature. One of the greatest advantages of polariton condensates
for optoelectronic applications is the easy implementation of arbitrary
geometries, or graphs, using adaptive optical elements for the excitation
laser beam Berloff2017a . We show controllable tuning between steady state
(single-colour) and limit-cycle (two-colour) regimes of the two coupled
condensate system and explain our observations through time-delayed equations
of motion.

## Results

### Time-delayed coupled matter-wave condensates

Figure 1: Macroscopically coupled matter-wave condensates. (a) Schematic
showing the cavity plane with two focusing laser beams (red cones) which are
spatially displaced by a distance $d\approx 114\;\mathrm{\upmu m}$ and the
experimentally observed photoluminescence of two synchronised, ballistically
expanding and interfering condensate centres. Finite particle transfer time
$\tau$ results in time-delayed interaction between condensates $\Psi_{1}$ and
$\Psi_{2}$. Unlike in the case of a conventional bosonic Josephson junction
(b) the reported time-delayed coupling mechanism between two tightly-pumped
polariton condensates (c) is not mediated by a tunneling current $J$ (blue
dashed line) but by a radiative transfer of particles (blue wavy line).

Conventionally, networks and lattices of condensates consisting of cold atoms
bloch2005ultracold , Cooper pairs cataliotti2001josephson , photons
dung2017variable  or polaritons kim2011dynamical ; Abbarchi_NatPhy2013  are
studied in trap geometries, weakly coupled via tunneling currents, that allow
for the study of solid state physics phenomena such as superfluid-Mott phase
transitions greiner2002quantum , magnetic frustration struck2011quantum ,
$PT$-symmetric non-linear optics Zhang_PRL2016 , and Josephson physics
albiez2005direct ; Abbarchi_NatPhy2013 . Here we investigate the inverse case,
wherein polariton condensates are freely expanding from small (point-like)
sources experiencing dynamics reminiscent to macroscopic systems such as time-
delayed coupled semiconductor lasers Soriano_RevModPhy2013 . Ballistic
expansion extends over two orders of magnitude beyond the pump beam waist, and
occurs due to the repulsive potential formed by the uncondensed exciton
reservoir injected by the non-resonant pump Wouters_2008PRB . In the case of
spatially separated polariton condensates, propagation from one condensate
centre to another results in a substantial phase accumulation, interpreted as
a retardation of information flow between the condensates. To date, coupled
polariton condensates were restricted to distances $d\lesssim
40\;\mathrm{\upmu m}$ christmann2014oscillatory ; Ohadi2016 ; Berloff2017a ,
wherein coupling was assumed to be instantaneous although propagation time was
comparable or longer than the characteristic timescale of the systems. To
unravel the role of time-delayed interactions, we demonstrate the
synchronisation between two tightly-pumped condensates separated by up to
$d\approx 114\;\mathrm{\upmu m}$, as shown in Fig. 1(a). Figures 1(b) and (c)
compare schematically the conventional regime of coupled condensates separated
by a potential barrier and described by a tunneling current $J$, to the
macroscopically coupled driven-dissipative matter-wave condensates interacting
via radiative particle transfer subject to finite propagation time $\tau$.

### Dynamics of two interacting condensates

We utilise a strain compensated semiconductor microcavity to enable
uninhibited ballistic expansion of polaritons over macroscopic distances
cilibrizzi2014polariton . We inject a polariton dyad with identical Gaussian
spatial profiles, with full-width-at-half-maximum (FWHM) of $\approx
2\;\mathrm{\upmu m}$, using non-resonant optical excitation at the first Bragg
minima of the reflectivity stop-band, and employ spatial light modulation to
continuously vary the separation distance from $6\;\mathrm{\upmu m}$ to
$94\;\mathrm{\upmu m}$, while keeping constant the excitation density of each
pump beam at $P_{1,2}\approx 1.5\times P_{\mathrm{thr}}^{(1)}$, where
$P_{\mathrm{thr}}^{(1)}$ is the condensation threshold power of an isolated
condensate. For each separation distance $d$ (more than 400 positions), we
record simultaneously the spatial profile of the photoluminescence (real-
space), the dispersion (energy vs in-plane momentum in the direction of
propagation), and the photoluminescence in reciprocal (Fourier) space (see
Supplementary Movie 1). We observe that opposite to the two-fold hybridisation
of two evanescently coupled condensates Abbarchi_NatPhy2013 , the
macroscopically coupled system is characterised by a multitude of accessible
modes of even and odd parity (i.e., $0$ and $\pi$ phase difference between
condensate centres) that alternate continuously between opposite parity states
with increasing dyad separation distance. For a range of separation distances
only one resonant mode is present in the gain region of the dyad, wherein the
polariton dyad is stationary, occupying a single energy level. Between the
separation distances, wherein only one mode is present, we observe the
coexistence of two resonant modes of opposite parity resulting in non-
stationary periodic states.

In Fig. 2 we present an example of the two different regimes, stationary and
non-stationary states. Figures 2(a,b) show real-space photoluminescence, Fig.
2(c,d) show Fourier-space photoluminescence and Fig. 2(e,f) show the polariton
dispersion along the axis of the dyad with separation distances of (a,c,e)
$d=12.7\;\mathrm{\upmu m}$ and (b,d,f) $d=37.3\;\mathrm{\upmu m}$. Fig. 2(g)
depicts the integrated spectra of Fig. 2(e,f) using black dots and red
squares, respectively. Absolute values of the energy levels are given as a
blueshift with respect to the ground-state of the lower polariton branch.
Although from the clear fringe pattern in Fig. 2(a) one could infer that only
one mode of even parity is present, it is through the simultaneous recording
of either the Fourier-space or the dispersion that the coexistence of an odd
parity mode becomes apparent. Each of the two modes has a well defined but
opposite parity to the other. We note that as long as we maintain the mirror
symmetry of the system by pumping both of the two condensate centres with the
same power $P_{1}=P_{2}$, we do not observe formation of non-trivial phase
configurations $\phi_{1}-\phi_{2}\neq 0,\pi$.

Figure 2: Emission characteristics of two non-resonantly pumped polariton
condensates. Measured (a,b) real-space photoluminescence $|\Psi(x,y)|^{2}$,
(c,d) momentum-space photoluminescence $|\Psi(k_{x},k_{y})|^{2}$ and (e,f)
spectrally-resolved momentum-space photoluminescence along $k_{y}=0$ for two
condensates with separation distances (a,c,e) $d=12.7\;\mathrm{\upmu m}$ and
(b,d,f) $d=37.3\;\mathrm{\upmu m}$, respectively. (g) Normalised spectra,
which are obtained by integrating (e) and (f), are illustrated with black dots
and red hollow squares. For better visibility of low-intensity features all
images are illustrated in logarithmic grey-scale saturated below 0.001 as
shown in (f). Scale bars in (a,c,e) correspond to ($10\;\mathrm{\upmu
m}$,$1\;\mathrm{\upmu m}^{-1}$, $1\;\mathrm{\upmu m}^{-1}$).

To unravel the dynamics of the coupling on the separation distance between the
two condensates, we obtain the spectral position, parity, and spectral weight
of both energy-levels for pump spot separation distances from
$d=5\;\mathrm{\upmu m}$ to $d=66\;\mathrm{\upmu m}$. We have recorded
configurations with more than two occupied energy levels for several distances
$d$, but with the relative spectral weight of the third peak less than a few
percent. In the following, we focus our analysis to the two brightest energy
levels for each configuration. The measured normalised spectral weights and
spectral positions of the two states versus pump spot separation distance $d$
are illustrated in Fig. 3 (a) and (b) using black dots for even parity states
and red hollow squares for odd parity states, respectively. Figure 3(b) shows
that the system follows an oscillatory behaviour in the spectral weights of
the two parity states giving rise to continuous transitions between even and
odd parity states interleaved by configurations exhibiting both parity states.
Interestingly, each period of these oscillations in relative intensity
(starting and ending with a vanishing spectral weight) displays an ‘energy
branch’ featuring a notable reduction in energy with increasing pump spot
separation distance $d$ as is shown in Fig. 3(b). The energy of an isolated
condensate is measured in the same sample space and its value $\approx
2.22\;\mathrm{meV}$, above the ground state energy of the lower polariton
dispersion, is illustrated with a blue dashed horizontal line. We observe that
the phase-coupling of two spatially separated condensate centres is dominated
by a spectral red-shift with respect to the energy of an isolated condensate.
The spectral size of each ‘energy branch’, i.e., the measurable red-shift, is
decaying from branch-to-branch with increasing pump spot separation distance
$d$. The two dyad configurations exhibiting single-colour and two-colour
states shown in Fig. 2 are indicated with gray vertical dashed lines in Fig.
3. Interestingly, it has been shown that the observation of phase-flip
transitions, accompanied with changes in oscillation frequency to another mode
are a universal characteristic of time-delayed coupled non-linear systems
schuster1989mutual ; prasad2008universal . Such dynamics can be associated
with neuronal systems adhikari2011time ; barardi2014phase ;
dotson2016experimental , and coupled semiconductor lasers
Soriano_RevModPhy2013 ; clerkin2014multistabilities , but have also been
demonstrated experimentally for other types of time-delayed coupled systems
such as non-linear electronic circuits reddy2000experimental ;
prasad2008universal , living organisms takamatsu2000time , chemical
oscillators cruz2010phase  and candle-flame oscillators manoj2018experimental
.

Figure 3:  Spectra of ballistically coupled polariton condensates. (a)
Spectral weight and (b) measured blueshift of even (black dots) and odd (red
hollow squares) parity states formed by the coupling of two spatially
separated polariton condensates with distance $d$. The energy level of a
single condensate pumped with the same non-resonant excitation power density,
i.e. $1.5$ times the condensation threshold of a single condensate, is
illustrated with a horizontal blue dashed line. The vertical dashed lines at
$d=12.7\;\mathrm{\upmu m}$ and $d=37.3\;\mathrm{\upmu m}$ indicate a two-
colour and single-colour state, respectively.

Time averaged measurements over $\sim 1000$ realisations of the system as
presented in Fig. 2(a,c,e) cannot reveal whether the system of two coupled
polariton condensates exhibits periodic dynamics in the form of two-colour
states or whether it stochastically picks one of the two parity states.
However, it is well known that the power spectrum of a time-dependent signal
is related to its auto-correlation function by means of Fourier-transform
(Wiener-Khinchin theorem). A two-colour solution, as shown in Fig. 2(g),
results in the periodic disappearance and revival of the first order temporal
coherence function $g^{(1)}(\bar{\tau})$. To characterise the temporal
evolution of $g^{(1)}(\bar{\tau})$, we perform interferometric measurements of
the two coupled condensates and extract the fringe visibility $\mathcal{V}$
which - up to a normalisation factor - is proportional to the magnitude of the
coherence function $\left|g^{(1)}\right|$. Assuming two coexisting states of
opposite parity but equal spectral weight, one may simplify the complex
amplitudes of the two coupled condensates $\Psi_{1}$ and $\Psi_{2}$ as

$\displaystyle\Psi_{1}(t)$
$\displaystyle=\psi_{0}\left(e^{-i\mu_{\mathrm{e}}t/\hbar}+e^{-i\mu_{\mathrm{o}}t/\hbar}e^{i\phi)}\right),$
(1) $\displaystyle\Psi_{2}(t)$
$\displaystyle=\psi_{0}\left(e^{-i\mu_{\mathrm{e}}t/\hbar}-e^{-i\mu_{\mathrm{o}}t/\hbar}e^{i\phi)}\right),$
(2)

with time-independent complex amplitude $\psi_{0}$ and relative phase $\phi$
between the two modes of opposite parity and energy
$\mu_{\mathrm{e}},\mu_{\mathrm{o}}$ for the even and odd mode respectively.
The mixture of two modes, with energy splitting
$\hbar\Delta=\left|\mu_{\mathrm{e}}-\mu_{\mathrm{o}}\right|$, causes an
antiphase temporal beating in the intensities $\left|\Psi_{1,2}\right|^{2}$
similar to oscillatory population transfer in coupled bosonic Josephson
junctions. Therefore, while the occupation amplitude of both condensate
centres oscillates with a period $T=2\pi/\Delta$, they feature a relative
phase shift of $\pi$ due to mixture of even and odd parity states. We use a
Michelson interferometer, comprising of a retroreflector mounted on a
translational stage and measure both the time averaged interference of the
photoluminescence of the same condensate
$I_{11}(\bar{\tau})=\left<\left|\Psi_{1}(t)+\Psi_{1}(t+\bar{\tau})\right|^{2}\right>$
as well as the time averaged interference of the photoluminescence of opposite
condensates
$I_{12}(\bar{\tau})=\left<\left|\Psi_{1}(t)+\Psi_{2}(t+\bar{\tau})\right|^{2}\right>$,
where $\bar{\tau}$ is the relative time-delay controlled by the adjustable
position of the retroreflector. Examples of the interferometric images are
illustrated in Fig. 4(a) for two condensates with a separation distance of
$d=10.3\;\mathrm{\upmu m}$, which demonstrates condensation in both an even
and an odd parity state as shown in Fig. 4(b). Assuming the spectral
composition as noted in Eqs. (1) and (2) the corresponding interferometric
visibilities $\mathcal{V}$ can be written as

$\displaystyle\mathcal{V}_{11}=\left|\cos(\Delta\bar{\tau}/2)\right|,$ (3)
$\displaystyle\mathcal{V}_{12}=\left|\sin(\Delta\bar{\tau}/2)\right|.$ (4)

The experimentally extracted visibilities $\mathcal{V}_{11}(\bar{\tau})$ and
$\mathcal{V}_{12}(\bar{\tau})$ versus relative time-delay $\bar{\tau}$ are
illustrated in Fig. 4(c) with blue circles and red squares, respectively. In
agreement with the predicted correlations [Eqs. (3) and (4)] we observe high
fringe visibility $\mathcal{V}_{11}(0)>0.8$ for the interference of one
condensate with itself and almost vanishing fringe visibility
$\mathcal{V}_{12}(0)<0.1$ for the interference of opposite condensates at zero
time-delay. Furthermore, we find periodic disappearance and revival of fringe
visibility with period $T\approx 15.3\;\mathrm{ps}$, which is consistent with
the observed energy-splitting of the two modes $\hbar\Delta=270\;\mathrm{\upmu
eV}$ [see Fig. 4(b)]. The relative phase shift between the two visibilities
$\mathcal{V}_{11}(\bar{\tau})$ and $\mathcal{V}_{12}(\bar{\tau})$ is a direct
result of the periodic population transfer between the two condensates. For
comparison the inset in Fig. 4(c) shows the expected visibilities [Eqs. (3)
and (4)] multiplied with an exponential decay accounting for the finite
coherence time of the system.

In Fig. 4(d) we show the decay of temporal coherence for two coupled
condensates in a single-colour state ($d=20\;\mathrm{\upmu m}$), a two-colour
state ($d=20.5\;\mathrm{\upmu m}$) and for an isolated condensate excited with
the same pump power density. We note that the coherent exchange of particles
in both regimes of the coupled condensate system results in an enhanced
coherence time. For the single-colour state ($d=20\;\mathrm{\upmu m}$) and the
isolated condensate the coherence time $\bar{\tau}_{\mathrm{c}}$ can be
extracted from exponential fits yielding $25.5\;\mathrm{ps}$ and
$10.2\;\mathrm{ps}$, respectively.

Figure 4: Periodic population transfer between two ballistically coupled
polariton condensates. The emission of a polariton dyad with condensation
centres $\Psi_{1,2}$ pumped with $P_{1,2}=1.5P^{(1)}_{\mathrm{thr}}$ at
locations $(x=\pm d/2,y=0)$ with $d\approx 10.3\;\upmu\mathrm{m}$ is
interferometrically and spectrally investigated in (a-c). (a) Time integrated
real space photoluminescence of both condensates and the interference patterns
of one condensate centre $\Psi_{1}(t)$ interfering with a delayed version of
itself $\Psi_{1}(t+\bar{\tau})$ or with a delayed version of the spatially
displaced condensate $\Psi_{2}(t+\bar{\tau})$. (b) Corresponding spectrally
resolved real-space photoluminescence along the axis of the dyad $y=0$. The
colour-scale of the normalised counts is saturated above $0.5$ for better
visibility. (c) Extracted fringe visibilities for the interference of same
condensate centres (blue hollow circles) and opposite condensate centres (red
hollow squares) versus relative time-delay $\bar{\tau}$ of the two
interferometer arms. The inset in (c) illustrates the expected visibilities
[Eq. (3) and (4)] multiplied with an exponentially decaying envelope. (d)
Comparison of the temporal decay of coherence extracted from the interference
of the same condensate
$\langle\left|\Psi_{1}(t)+\Psi_{1}(t+\bar{\tau})\right|^{2}\rangle$ for two
coupled condensates with separation distances $d=20\;\upmu\mathrm{m}$ and
$d=20.5\;\upmu\mathrm{m}$, as well as a single isolated condensate. In all
three cases each condensate is pumped equally with $P\approx
1.7P^{(1)}_{\mathrm{thr}}$. The error bars are calculated as the standard
deviation of the extracted visibility within the full width at half maximum
($\approx 2\;\mathrm{\upmu m}$) of each condensate. Lines represent
exponential fits for the decay of coherence for the single-mode polariton dyad
and the single condensate, respectively.

### Numerical analysis

The dynamics of polariton condensates can be modelled via the mean field
theory approach where the condensate order parameter $\Psi(\mathbf{r},t)$ is
described by a 2D semiclassical wave equation often referred as the
generalised Gross-Pitaevskii equation coupled with an excitonic reservoir
which feeds non-condensed particles to the condensate Wouters_PRL2007 . In
Supplementary Note 1 we give numerical results of the spatiotemporal dynamics
of the dyad reproducing semi-quantitatively the experimentally observed
results. We also provide a supplemental animation showing the calculated
evolution of the two-colour condensate (Supplementary Movie 2).

In the following, we show that our experimental observations are described as
a system of time-delayed coupled non-linear oscillators. For simplicity we
consider the 1D Schrödinger equation corresponding to the problem of $s$-wave
scattering of the condensate wavefunctions with $\mathrm{\delta}$-shaped
complex-valued, pump induced, potentials. We start by characterising the
energies of the time-independent non-hermitian single particle problem:

$E\Psi(x)=\left(-\frac{\hbar^{2}\partial_{x}^{2}}{2m}+V(x)-i\frac{\hbar\gamma_{\mathrm{c}}}{2}\right)\Psi(x),$
(5)

where $V(x)$ describes complex-valued $\delta$-shaped potentials separated by
a distance $d$, $V(x)=V_{0}\left(\delta(x+d/2)+\delta(x-d/2)\right)$, and
$V_{0}$ lies in the first quadrant of the complex plane, i.e. repulsive
interactions and gain. The eigenfunctions of Eq. (5) describing normalisable
solutions of outwards propagating waves from the potential (non-resonant pump)
centres are written as,

$\Psi(x)=\begin{cases}Ae^{-ikx},&x\leq-d/2\\\ Be^{ikx}+Ce^{-ikx},&|x|<d/2\\\
De^{ikx},&x\geq d/2\end{cases}$ (6)

where $k$ also belongs to the first quadrant of the complex plane. This
problem is well known for the case of lossless attractors ($\Re{(V_{0})}<0$,
$\Im{(V_{0})}=0$) describing electron states in a 1D diatomic Hydrogen
molecule ion. The resonance condition of the system is written,

$\left[\frac{V_{0}}{\frac{i\hbar^{2}k}{m}-V_{0}}e^{ikd}\right]^{2}=1.$ (7)

The solutions of Eq. (S10) can explicitly be written as,

Figure 5:  Resonances of two complex-valued $\delta$-shaped potentials.
Calculated (a) imaginary part and (b) real part of the eigenvalues $E$ in Eq.
(9) for $\Re{(V_{0})}=1\;\mathrm{meV}\;\mathrm{\upmu m}$,
$\Im{(V_{0})}=2\mathrm{meV}\;\mathrm{\upmu m}$,
$m=0.28\;\mathrm{meV}\;\mathrm{ps}^{2}\;\mathrm{\upmu m}^{-2}$,
$\gamma_{\mathrm{c}}=1/5.5\;\mathrm{ps}^{-1}$. The two solutions with largest
gain (imaginary part of $E$) are illustrated as lines, other solutions are
shown as dotted branches.

$k_{n,\pm}=-i\tilde{V}+\frac{i}{d}W_{n}(\mp d\tilde{V}e^{d\tilde{V}}),\quad
n\in\mathbb{Z},$ (8)

with $\tilde{V}=mV_{0}/\hbar^{2}$. Equation (8) describes infinitely many
solutions of the system of even $(+)$ and odd $(-)$ parity, where $W_{n}$ are
the branches of the Lambert $W$ function. The corresponding complex-valued
eigenvalues,

$E_{n,\pm}=\frac{\hbar^{2}k_{n,\pm}^{2}}{2m}-\frac{i\hbar\gamma_{\mathrm{c}}}{2},$
(9)

are illustrated in Fig. 5 and qualitatively reproduce the experimental
findings of the multiple energy branches shown in Fig. 3. We interpret the
experimental occurrence of predominantly two lasing modes with the behaviour
of the imaginary values of $E_{n}$ for the simplistic 1D-toy model. While
there are periodically alternating regions of even and odd-parity solutions
dominating the gain, we expect distances at which two modes (of opposite
parity) have equal gain and thus can operate with equal intensity. It is worth
noting that the Lambert $W$ function naturally arises for problems involving
delay differential equations Asl2003 .

The time-dependent problem can be formulated as a superposition of two
displaced, normalisable, and ballistically propagating waves $\psi_{1,2}(x)$,
each emerging from one of the condensate centres,

$\Psi(x,t)=c_{1}(t)\psi_{1}(x)+c_{2}(t)\psi_{2}(x).$ (10)

Here, in analogy with Eq. (6), the normalised ansatz is written as
$\psi_{1,2}(x)=\sqrt{\kappa}e^{ik|x\pm d/2|}$, with a complex-valued
wavevector $k=k_{\mathrm{c}}+i\kappa$. When the coupling between condensates
is weak, i.e. small $\xi=\exp(-\kappa d)$, one can omit all terms of order
$\mathcal{O}(\xi^{2})$ and higher. Then plugging Eq. (10) into the time-
dependent form of Eq. (5) and integrating out the spatial degrees of freedom,
assuming that Eq. (S10) is satisfied, one gets (see Supplementary Note 2),

$i\hbar\dot{c}_{i}=\left[\frac{\hbar^{2}k^{2}}{2m}-i\frac{\hbar\gamma_{\mathrm{c}}}{2}+\kappa\left(V_{0}-\frac{i\hbar^{2}k}{m}\right)\right]c_{i}+V_{0}\kappa
e^{ikd}c_{j}.$ (11)

where $j=3-i$ and $i=1,2$ are the condensate indices. When setting $c_{i}=\pm
c_{j}$, and solving for stationary states, the above equation recovers the
exact resonant solutions dictated by Eq. (8). Equation (11) then shows that
inter-condensate interaction is in the form of a coherent influx of particles
from condensate $j$ onto the condensate centre of condensate $i$ (and vice
versa), with a phase retardation of $k_{\mathrm{c}}d$. When $c_{i}$ and
$c_{j}$ oscillate at a fixed frequency $\omega$ we can transform the phase-
shifting term $\exp(ik_{\mathrm{c}}d)$ into an effective time-delay,

$e^{ikd}c_{j}(t)=e^{-\kappa d}c_{j}(t-\tau).$ (12)

The time-delay $\tau$ corresponds to an interaction lag between condensate
centres $i$ and $j$ caused by their spatial separation and is given by
$\tau=k_{\mathrm{c}}d/\omega$. In the case of weak coupling, i.e. small
changes in oscillation frequency $\omega$ compared to the frequency
$\omega_{0}$ of a single unperturbed condensate, the time-delay is
approximately proportional to the dyad separation distance $d$, i.e.
$\tau\approx k_{\mathrm{c},0}d/\omega_{0}$ where we also use the notation
$k_{0}=k_{\mathrm{c},0}+i\kappa_{0}$ for the wavevector of the single
unperturbed condensate. The exponentially decaying term on the right-hand side
of Eq. (12) accounts for the one-dimensional spatial decay of particles
propagating in between the two condensate centres. Introducing local non-
linear interactions, reservoir gain and blueshift one can write the full non-
linear equation of motion for the two coupled condensates as Wouters_PRL2007 ,

$\displaystyle i\dot{c}_{i}$
$\displaystyle=\left[\Omega+\left(g+i\frac{R}{2}\right)n_{i}+\alpha|c_{i}|^{2}\right]c_{i}+Je^{i\beta}c_{j}(t-d/v),$
(13) $\displaystyle\dot{n}_{i}$
$\displaystyle=-(\Gamma_{\mathrm{A}}+R|c_{i}|^{2})n_{i}+P.$ (14)

Here, $n_{i}$ correspond to the pump induced exciton reservoirs providing
blushift and gain into their respective condensates, $g$ is the polariton-
reservoir interaction strength, $R$ is the rate of stimulated scattering of
polaritons into the condensate from the active reservoir, $\alpha$ is the
interaction strength of two polaritons in the condensate,
$v=\omega_{0}/k_{\mathrm{c},0}$ is the phase velocity of the polariton
wavefunction, and $\Gamma_{\mathrm{A}}$ is the radiative decay rate of the
bottleneck reservoir excitons. The parameter $\Omega=\Omega_{0}-i\Gamma$
captures the self energy of each condensate which will, in general, have a
contribution from a background of optically inactive dark excitons generated
at the pump spot, and $\Gamma$ denotes the effective linewidth of the
polaritons expanding away from the pump spot. The complex valued coupling is
written $J\exp{(i\beta)}=V_{0}\kappa\exp{(-\kappa d)}$ for brevity. Equation
(13) is then in the form of a discretised complex Gross-Pitaevskii equation
but with time-delayed interaction between the bosonic particle ensembles;
which greatly increases the dimensionality of phase-space and complexity of
the coupled system. Our system then has strong similarity with the famous
Lang-Kobayashi equation Kobayashi_JQE1980 ; kozyreff2000global  where in our
case each condensate acts as a radiating antenna of symmetrically expanding
waves. The two spatially separated antennas interfere and maximise their gain
by adjusting both their common frequency and their relative phase difference
$\phi=\arg{(c_{i}^{*}c_{j})}$. Similarities of the dynamics of coupled
polariton condensates to equations of motion in the form of the Lang-Kobayashi
equation was discussed in theoretical works recently for instantaneously
coupled condensates with complex-valued couplings kalinin2019polaritonic .
Unlike trapped ground state bosonic systems, such as cold atoms, the
ballistically expanding polariton matter-wave condensates necessarily
experience time-delayed coupling since $k_{\mathrm{c},0}d\gg 1$ similar to
inter-cavity coupling of semiconductor lasers Soriano_RevModPhy2013 .

In order to accurately reproduce experimentally observed spectra using Eqs.
(13) and (14) we fit the distance dependence of the coupling amplitude $J(d)$
to the spatial envelope of a single condensate. It is known that in a 2D
system the $0$-order Hankel function of the first kind describes the
cylindrically symmetric radial outflow of particles in the linear regime
Wouters_2008PRB , therefore we choose

$J(d)=J_{0}\left|H^{(1)}_{0}(k_{0}d)\right|,$ (15)

where $J_{0}$ is a parameter describing the coupling strength (see
Supplementary Note 3). Numerical integration of Eqs. (13) and (14) is computed
for separation distances $d>10\;\mathrm{\upmu m}$ using Gaussian white noise
as initial conditions which, in close analogy to the experimental findings,
reproduces periodic parity-flip transitions accompanied by cyclic solutions in
the transition region. Fig. 6(a) depicts the experimentally measured and
normalised emission spectra of the polariton dyad versus pump spot separation
distance $d$ on a grey-scale colour-map, for which the extracted two most
dominant spectral peaks are depicted in Fig. 3. The red dots represent the
numerically calculated spectral peak from Eqs. (13) when in single mode
operation, or the two most dominant spectral peaks when multiple spectral
components exist; showing excellent agreement with experiment. As an example,
we illustrate a continuous transition of the system from an anti-phase to in-
phase state described by the time-delayed coupled model in Fig. 6(b) and (c)
showing the corresponding spectral decomposition and phase-space diagrams of
the system for an increasing set of distances $d$ from $20\;\mathrm{\upmu m}$
to $21\;\mathrm{\upmu m}$. Similar stationary and oscillatory behaviour, for
separation distances $d=20\;\mathrm{\upmu m}$ and $d=20.5\;\mathrm{\upmu m}$
respectively, is shown in Fig. 4(d). The numerically simulated phase-space
diagrams depict periodic orbits in the transition region, involving periodic
oscillations of the phase difference $\phi=\arg(c_{1}^{*}c_{2})$ and
population imbalance $z=(|c_{1}|^{2}-|c_{2}|^{2})/(|c_{1}|^{2}+|c_{2}|^{2})$,
which is also confirmed by 2D-simulations of the generalised Gross-Pitaevskii
equation (see Supplementary Note 1). We have verified through numerics that
time-delay physics are indeed an accurate representation of the coupled
condensate dynamics (see Supplementary Note 2).

Figure 6: Periodic parity-transitions described by time-delayed coupled
oscillators. (a) Comparison of experimentally measured spectra of two coupled
polariton condensates coupled with separation distance $d$ and numerically
simulated spectral peaks using Eqs. (13), (14) and (15). Experimentally
measured spectra are normalised for each distance $d$ and the grey-scale
colour-map is saturated above $0.5$ for better visibility. (b) and (c)
illustrate the spectral decomposition and phase-diagram for the anti-phase to
in-phase transition of the system by increasing the distance $d$ from
$20\;\mathrm{\upmu m}$ to $21\;\mathrm{\upmu m}$. Simulation parameters:
$\hbar\Omega=(1.22-i0.5)\;\mathrm{meV}$, $\hbar\alpha=0.1\;\mathrm{\upmu eV}$,
$\hbar R=0.5\;\mathrm{\upmu eV}$, $\hbar g=0.5\;\mathrm{\upmu eV}$,
$v=1.9\;\mathrm{\upmu m}\;\mathrm{ps}^{-1}$, $P=100\;\mathrm{ps}^{-1}$,
$\Gamma_{\mathrm{A}}=0.05\;\mathrm{ps}^{-1}$, $\hbar J_{0}=1.1\;\mathrm{meV}$,
$k_{0}=(1.7+i0.014)\mathrm{\upmu m}^{-1}$ and $\beta=-1$.

We note that in the limit of fast active reservoir relaxation,
$\Gamma_{\mathrm{A}}^{-1}\ll\Gamma^{-1}$, one can adiabatically eliminate Eq.
(14) and introduce an effective non-linear term to Eq. (13),
$(\alpha_{\text{eff}}-i\sigma)|c_{i}|^{2}$, accounting for polariton-polariton
and polariton-reservoir interactions, as well as condensate gain saturation.
The dynamical equations of two coupled condensates are then described by time-
delayed coupled Stuart-Landau oscillators,

$i\dot{c}_{i}=\left(\Omega+(\alpha_{\text{eff}}-i\sigma)|c_{i}|^{2}\right)c_{i}+Je^{i\beta}c_{j}(t-d/v).$
(16)

Numerical analysis of Eq. (16) is given in Supplementary Note 4 showing
qualitative agreement with experiment.

## Discussion

We present an extensive experimental and theoretical study of a system of two
ballistically expanding (untrapped) interacting polariton condensates, the
fundamental building block of polariton graphs with higher connectivity. We
demonstrate a regime for coupled matter-wave condensates, wherein the coupling
is not instantaneous but mediated by a particle flow inherently connected with
time retardation effects. We observe deterministic selection of steady state
(single-colour) or dynamical (two-colour) modes of the system by controlling
the separation distance between the condensates. Time-delay polaritonics
potentially offers an ultrafast platform for simulating the dynamics of real-
world systems of time-delayed coupled non-linear oscillators that appear in
photonics, electronics and neural circuits. Given the high non-linearities of
polaritons and the ease of optically imprinting multiple condensates of
arbitrary geometries on planar microcavities, the system offers promising
applications for neuromorphic devices based on lattices of history dependent,
non-trapped, strongly interacting, polariton condensates with a wide range of
coupling strengths, fast optical operations (input), dynamics (processing),
and readout.

## Methods

### Microcavity sample and experimental methods

The microcavity sample used is a strain compensated $2\lambda$ GaAs
microcavity with three pairs of 6 nm InGaAs quantum wells embedded at the
anti-nodes of the electric field cilibrizzi2014polariton . The intracavity
layer contains a wedge and the position on the sample is chosen such that the
cavity photonic mode is red-shifted from the excitonic mode at zero inplane
wavevector ($|\mathbf{k}|=0$) by $\approx-5.5$ meV. For all experiments the
microcavity is held in a cold finger cryostat (temperature $T\approx 6\;K$)
and is optically pumped with a circularly polarised continuous wave monomode
laser blue-detuned above the cavity stopband ($\lambda\approx
785\;\mathrm{nm}$). To prevent heating of the sample an acousto-optic
modulator is used to generate square wave packets at a frequency of 10 kHz and
duty cycle of $5\%$. A liquid crystal spatial light modulator (SLM) imprints a
phase pattern such that when the beam is focused through the 0.4 numerical
aperture microscope objective lens it excites the sample with the desired
spatial geometry. The phase patterns are carefully designed so that when
changing the pump separation distance $d$ of the two-spot excitation pattern,
the diffraction efficiency of the SLM remains constant and both pump spots
retain equal excitation power and width. The photoluminescence is collected in
reflection geometry and spectrally resolved using an 1800 grooves/mm grating
in a $750\;\mathrm{mm}$ spectrometer, which is equipped with a charge-coupled
device (CCD). Real-space, Fourier-space and dispersion images are acquired
using exposure times in the order of milliseconds.

### Interferometry

A modified Michelson-interferometer, where one mirror is replaced with a
retroreflector mounted on a translational stage, is used for measuring the
interference and temporal coherence of the emission. By tilting the angle of
the emission entering the interferometer we select to spatially overlap the
emission of either opposing condensates or the same condensates onto a CCD
camera. Interference fringe visibilities are extracted from the normalised
$1^{\mathrm{st}}$ diffraction order of the computed discrete Fourier transform
of each interference image.

### Image processing

Image displayed in Figure 1(a) has been digitally processed with a low-pass
filter to increase the visibility of interference fringes.

## Acknowledgements

The authors are grateful to N. G. Berloff for fruitful discussions and
acknowledge the support of the UK’s Engineering and Physical Sciences Research
Council (grant EP/M025330/1 on Hybrid Polaritonics).

## Author contributions

P.G.L. led the research project. P.G.L., J.D.T and L.P designed the
experiment. J.D.T. and L.P. carried out the experiments and analysed the data.
J.D.T and H.S. developed the theoretical modelling. H.S performed numerical
simulations. All authors contributed to the writing of the manuscript.

## Data availability

The data that support the findings of this study are openly available from the
University of Southampton repository (DOI: 10.5258/SOTON/D1149) dataset .

## Competing interests

The authors declare no competing interests.

Supplemental Information

Figure S1: (a-b) Calculated relative visibility and spectral position of even
(red) and odd (blue) parity states formed by the coupling of two spatially
separated polariton condensates for $P_{0}=1.1P_{\text{thr}}^{(1)}$. The
measured energy level of a single condensate pumped with the same non-resonant
excitation power is illustrated with a black horizontal dashed line. (c-f)
Spatial density of the synchronised dyad $|\Psi(\mathbf{r},t)|^{2}$ with size
$d=12\;\mathrm{\upmu m}$ calculated at $P_{0}=1.5P_{\text{thr}}^{(1)}$.
Colour-scale is linear and saturated at 30% maximum intensity. (g) Spectrally
resolved momentum space image of the dyad wavefunction density. Colour-scale
is logarithmic. (h) Trajectory of the condensates density imbalance $z$ and
relative phase $\phi_{1}-\phi_{2}$ forming a closed loop indicating persistent
oscillations with a period $T\approx 10\;\mathrm{ps}$. Vertical dashed line in
(a,b) indicates the spectral decomposition of the dyad with
$d=12\;\mathrm{\upmu m}$.  Figure S2: Spectra of two coupled condensates.
Comparison between normalised spectral intensity from (a) experiment, (b) Eqs.
(13) and (14) [main text], and (c) Eqs. (S13) and (S14). Red squares and blue
circles denote the only stable energies of in-phase (blue) and antiphase (red)
solutions determined from Bogoliubov-de-Gennes stability analysis. Simulation
parameters: $\hbar\Omega=(1.22-i0.5)\;\mathrm{meV}$,
$\hbar\alpha=0.1\;\mathrm{\upmu eV}$, $\hbar R=0.5\;\mathrm{\upmu eV}$, $\hbar
g=0.5\;\mathrm{\upmu eV}$, $v=1.9\;\mathrm{\upmu m}\;\mathrm{ps}^{-1}$,
$P=100\;\mathrm{ps}^{-1}$, $\Gamma_{\mathrm{A}}=0.05\;\mathrm{ps}^{-1}$,
$\hbar J_{0}=1.1\;\mathrm{meV}$, $k_{0}=(1.7+i0.014)\mathrm{\upmu m}^{-1}$ and
$\beta=-1$. Figure S3: Spatial emission characteristics of a single polariton
condensate pumped at $1.5$ times the condensation threshold. (a) Radially-
resolved photoluminescence in momentum-space visualises a radial outflow from
the condensate centre with well-defined wavevector
$k_{\mathrm{c}}=1.7\;\mathrm{\upmu m}^{-1}$. (b) Measured radial expansion in
real-space shows good agreement with a 0-order Hankel-function
$H_{0}^{(1)}(kr)$ for distances $r\geq 5\;\mathrm{\upmu m}$ from the
condensate centre. Figure S4: Time-delayed coupled Stuart-Landau oscillators.
Comparison of experimentally measured spectra (grey-scale image) of two
coupled polariton condensates coupled with separation distance $d$ and
numerically calculated fixed-point solutions using a system of two time-
delayed coupled Stuart-Landau oscillators (red dots) [see Eq. (16) in main
text]. Parameters: $\hbar\Omega=2.22\;\mathrm{meV}$,
$\hbar\alpha_{\text{eff}}=-0.0229\;\mathrm{meV}$,
$\hbar\sigma=0.0115\;\mathrm{meV}$,
$v=1.9\;\upmu\mathrm{m}\;\mathrm{ps}^{-1}$, $\hbar J_{0}=1.1\;\mathrm{meV}$,
$k_{0}=(1.7+i0.014)\upmu\mathrm{m}^{-1}$ and $\beta=-1$.

## Supplementary Note 1: Numerical spatiotemporal simulations

The dynamics of polariton condensates can be modelled via the mean field
theory approach where the condensate order parameter $\Psi(\mathbf{r},t)$ is
described by a 2D semiclassical wave equation often referred as the
generalised Gross-Pitaevskii equation coupled with an excitonic reservoir
which feeds non-condensed particles to the condensate Wouters_PRL2007 . The
reservoir is divided into two parts: An active reservoir
$n_{\mathrm{A}}(\mathbf{r},t)$ belonging to excitons which experience bosonic
stimulated scattering into the condensate, and an inactive reservoir
$n_{\mathrm{I}}(\mathbf{r},t)$ which sustains the active reservoir
Lagoudakis2010a ; Lagoudakis_PRL2011 .

$\displaystyle i\frac{\partial\Psi}{\partial t}$
$\displaystyle=\left[-\frac{\hbar\nabla^{2}}{2m}+\frac{g}{2}(n_{\mathrm{A}}+n_{\mathrm{I}})+\frac{\alpha}{2}|\Psi|^{2}+\frac{i}{2}\left(Rn_{\mathrm{A}}-\gamma\right)\right]\Psi,$
(S1) $\displaystyle\frac{\partial n_{\mathrm{A}}}{\partial t}$
$\displaystyle=-\left(\Gamma_{\mathrm{A}}+R|\Psi|^{2}\right)n_{\mathrm{A}}+Wn_{\mathrm{I}},$
(S2) $\displaystyle\frac{\partial n_{\mathrm{I}}}{\partial t}$
$\displaystyle=-\left(\Gamma_{\mathrm{I}}+W\right)n_{\mathrm{I}}+P(\mathbf{r}).$
(S3)

Here, $m$ is the effective mass of a polariton in the lower dispersion branch,
$\alpha$ is the interaction strength of two polaritons in the condensate, $g$
is the polariton-reservoir interaction strength, $R$ is the rate of stimulated
scattering of polaritons into the condensate from the active reservoir,
$\gamma$ is the polariton decay rate, $\Gamma_{{\mathrm{A,I}}}$ is the decay
rate of active and inactive reservoir excitons respectively, $W$ is the
conversion rate between inactive and active reservoir excitons, and
$P(\mathbf{r})$ is the non-resonant CW pump profile.

We perform numerical integration of Eqs. (S1), (S2) and (S3) in time using a
linear multistep method in time and spectral methods in space. The polariton
mass and lifetime are based on the sample properties: $m=0.28$ meV ps2
$\upmu$m-2 and $\gamma=1/5.5$ ps-1. We choose values of interaction strengths
typical of InGaAs based systems: $\hbar\alpha=3.3$ $\upmu$eV $\upmu$m2,
$g=20\alpha$. The choice of $g>\alpha$ stems from negatively detuned cavity
where polaritons in the condensate are more photonic and thus interact weakly
together than those directly with the incoherent reservoir. The non-radiative
recombination rate of inactive reservoir excitons is much smaller than the
condensate decay rate $\Gamma_{\mathrm{I}}^{-1}=500$ ps whereas the active
reservoir is more ambiguous. It has been argued that the active reservoir
(also referred as bottleneck polaritons) decay rate should be larger than the
condensate decay rate due to fast thermalisation to the exciton background
Wouters_2008PRB  and partly because they can decay radiatively. The choice of
this value often depends on the type of experiment and ranges over several
orders of magnitude Lagoudakis_PRL2011 ; Lagoudakis2010a . In the current
study we choose an intermediate value of $\Gamma_{\mathrm{A}}^{-1}=20$ ps and
find it produces results in good agreement with experiment. It should be noted
that the results are not highly sensitive to the exact value of
$\Gamma_{\mathrm{A}}$. The final two parameters are then found by fitting
experimental results of a single condensate, i.e., the onset of a sharp ring
in $k$-space around $k_{\mathrm{c}}\approx 1.7$ $\upmu$m-1 and a blueshift of
$\sim 500$ $\upmu$eV when raising the pump power from threshold to twice
threshold power $P_{\text{thr}}^{(1)}$ for a single isolated condensate. The
above two measured features of a single condensate are reproduced using the
values $\hbar R=129$ $\upmu$eV $\upmu$m-2, and $W=0.1$ ps-1. The pump is
written $P(\mathbf{r})=P_{0}e^{-r^{2}/2w^{2}}$ where $P_{0}$ denotes the pump
power and $w$ corresponds to a $2\;\mathrm{\upmu m}$ full width at half
maximum.

A single numerical scan (i.e., no averaging over many stochastic initial
conditions) of the dyad distance dependence is shown in S1(a-b) and reproduces
semi-quantitatively the results of Fig. 3 [main text]. It is worth noting that
the current model overestimates the stability of the single energy state in
the dyad seen from the plateaus appearing at 0 and 1 in the calculated
relative visibility. This can possible stem from the lack of energy relaxation
mechanisms, parameter dependence on the polariton Hopfield fraction, and/or
overestimation of the active reservoir depletion which can be artificially
tuned Lagoudakis2010a ; Wouters_2010PRB . Nevertheless, regimes of one-colour
and two-colour operation are clearly visible in simulations. All simulations
are performed with some natural weak disorder present in the system to account
for realistic non-ideal cavity conditions.

In S1(c-h) we show the condensate wavefunction calculated at
$d=12.7\;\mathrm{\upmu m}$ at $P_{0}=1.5P_{\text{thr}}^{(1)}$. Real space
dynamics shown in panels (c)-(f) are characterised by periodic beatings in the
spatial intensity of the polariton dyad, and the spectrally resolved momentum
space image of the wavefunction (g) shows occupancy of two dominant energy
levels. In S1(h) we plot the normalised density imbalance
$z=(\rho_{1}-\rho_{2})/(\rho_{1}+\rho_{2})$ against the condensate relative
phase. Here $\rho_{1,2}=|\Psi(\mathbf{r}-\mathbf{r}_{1,2})|^{2}$ and
$\mathbf{r}_{1,2}$ are the location coordinates of the left and the right pump
respectively and $\phi_{1,2}=\arg{[\Psi(\mathbf{r}-\mathbf{r}_{1,2})]}$. The
results show a closed trajectory corresponding to persistent oscillations in
both phase and particle density with a period $T\approx 10$ ps.

## Supplementary Note 2: Derivation of time delay equation

Plugging Eq. (10) [main text] into the time-dependent form of Eq. (5) [main
text],

$i\hbar\frac{d\Psi(x)}{dt}=\left(-\frac{\hbar^{2}\partial_{x}^{2}}{2m}+V(x)-i\frac{\hbar\gamma_{\mathrm{c}}}{2}\right)\Psi(x),$
(S4)

and integrating out the spatial degrees of freedom gives the following
algebraic differential equation,

$i\hbar\dot{c}_{i}+i\hbar\mathcal{N}_{12}\dot{c}_{j}=F(c_{i},c_{j}),$ (S5)

where $j=3-i$ ($i=1,2$). The overlap between the two wavefunctions is written,

$\mathcal{N}_{12}=e^{-\kappa
d}\left(\cos(k_{\mathrm{c}}d)+\frac{\kappa}{k_{\mathrm{c}}}\sin(k_{\mathrm{c}}d)\right)\simeq
e^{-\kappa d}\cos(k_{\mathrm{c}}d),$ (S6)

where $\kappa/k_{\mathrm{c}}\ll 1$ since the chemical potential of the
condensate is much larger than its linewidth. The right-hand side of Eq. (S5)
is written,

$\displaystyle F(c_{i},c_{j})$
$\displaystyle=\left(\frac{\hbar^{2}k^{2}}{2m}-i\frac{\hbar\gamma_{\mathrm{c}}}{2}\right)\left(c_{i}+\mathcal{N}_{12}c_{j}\right)+\left(\left[V_{0}-\frac{i\hbar^{2}k}{m}\right]c_{i}+V_{0}e^{ikd}c_{j}\right)\kappa$
$\displaystyle+\left(\left[V_{0}-\frac{i\hbar^{2}k}{m}\right]c_{j}+V_{0}e^{ikd}c_{i}\right)\kappa
e^{-ik^{*}d}.$ (S7)

We can decouple the above algebraic differential equation by taking into
account that the mass matrix,

$\hat{M}=\begin{pmatrix}1&\mathcal{N}_{12}\\\
\mathcal{N}_{12}&1\end{pmatrix},$ (S8)

has a well defined inverse since $\mathcal{N}_{12}<1$ for $d>0$. This allows
us to write,

$i\hbar\begin{pmatrix}\dot{c}_{1}\\\
\dot{c}_{2}\end{pmatrix}=\hat{M}^{-1}\begin{pmatrix}F(c_{1},c_{2})\\\
F(c_{2},c_{1})\end{pmatrix}.$ (S9)

Assuming that that the resonance condition is satisfied [Eq.(7) in main text],

$\left[\frac{V_{0}}{\frac{i\hbar^{2}k}{m}-V_{0}}e^{ikd}\right]^{2}=1,$ (S10)

and that coupling between condensates is weak, i.e., small $\xi=\exp(-\kappa
d)$, one can omit all terms of order $\mathcal{O}(\xi^{2})$ and higher. This
then gives,

$i\hbar\dot{c}_{i}=\left[\frac{\hbar^{2}k^{2}}{2m}-i\frac{\hbar\gamma_{\mathrm{c}}}{2}+\kappa\left(V_{0}-\frac{i\hbar^{2}k}{m}\right)\right]c_{i}+V_{0}\kappa
e^{ikd}c_{j}.$ (S11)

It is worth noting that without the resonance condition [Eq. (S10)], instead
of arriving at Eq. (S11) one would get,

$i\hbar\dot{c}_{i}=\left[\frac{\hbar^{2}k^{2}}{2m}-i\frac{\hbar\gamma_{\mathrm{c}}}{2}+\kappa\left(V_{0}-\frac{i\hbar^{2}k}{m}\right)\right]c_{i}+\left(V_{0}\cos{(k_{\mathrm{c}}d)}-\frac{\hbar^{2}k}{m}\sin{(k_{\mathrm{c}}d)}\right)e^{-\kappa
d}\kappa c_{j}.$ (S12)

Where the inter-condensate coupling term in Eq. (S12) is similar to previously
established results Ohadi2016 ; Berloff2017a , but with one critical addition,
the complex wavevector $k$ is innately related to the potential $V_{0}$.

To illustrate the difference between equations of motion with and without
time-delay we replot a portion of the experimental and numerical spectra from
Fig.6 [main text] in S2(a,b). In S2(c) we then plot the spectra coming from
two coupled equations similar to Eqs. 13 and 14 [main text] but but now
written:

$\displaystyle i\dot{c}_{i}$
$\displaystyle=\left[\Omega+\left(g+i\frac{R}{2}\right)n_{i}+\alpha|c_{i}|^{2}\right]c_{i}+Je^{i\beta}c_{j},$
(S13) $\displaystyle\dot{n}_{i}$
$\displaystyle=-(\Gamma_{\mathrm{A}}+R|c_{i}|^{2})n_{i}+P.$ (S14)

Here $J=J_{0}H_{0}^{(1)}(k_{0}d)$ is the zeroth order Hankel function of the
first kind describing the interference experienced between the two condensates
as a function of distance $d$. Such a model was recently theoretically
investigated in the context of polariton condensates kalinin2019polaritonic
and its relatability to other models in different universality classes,
including time-delay models. A major notable feature missing from the no-time-
delay model [Eqs. (S13) and (S14) ] is the overlap between different energy
branches indicating cyclical/non-stationary states and the cusp forming at the
ends of the branches. Indeed, for small non-linearities ($|c_{i}|^{2}\simeq
0$) Eq. (S13) is just a two level model with in-phase and anti-phase
eigensolutions $c_{\mathrm{L,U}}=(c_{1}\pm c_{2})/\sqrt{2}$ whose stability in
the non-linear regime ($|c_{i}|^{2}\neq 0$) can be studied through the
Bogoliubov-de-Gennes method. Equation (S13) also supports non-stationary
states at larger $P$ but the spectra qualitatively different from the
experimental observations. We therefore believe that time-delay physics as
given by Eq. 13 [main text] are an accurate representation of the coupling
between ballistically expanding condensates. It describes the effects of
interference and how the renormalised self-energy modifies the outflow
polariton momentum, giving rise to stable two-colour operation over wide
distance intervals.

## Supplementary Note 3: Single condensate

We characterise the luminescence of a single condensate pumped in the same
sample area and with the same pump power, i.e. $1.5$ times the condensation
threshold of a single condensate, as the two condensates forming the polariton
dyad presented in Fig. 2 and Fig. 3 [main text]. The measured and radially-
resolved intensity of far-field and near-field emission are illustrated in
S3(a) and (b), respectively. From (a) we can identify a well-defined outflow-
wavevector $k_{\mathrm{c}}=1.7\;\mathrm{\upmu m}^{-1}$. Further, by fitting
the radial intensity in (b) for $r>5\;\mathrm{\upmu m}$ with the squared
magnitude of the 0-order Hankel function of the first kind
$|H_{0}^{(1)}(kr)|^{2}$ we obtain the complex-valued wavevector
$k=(1.7+i0.014)\;\mathrm{\upmu m}^{-1}$.

## Supplementary Note 4: Time-delayed coupled Stuart-Landau oscillators

S4 illustrates the numerically calculated brightest fixed point solutions (red
dots) of the time-delayed coupled Stuart-Landau oscillators [Eq. (16) in main
text] versus condensate separation distance $d$ matching the experimentally
measured and normalised spectra, illustrated in grey-scale. The observed red-
shift of the coupled condensate system with respect to the isolated condensate
energy level at $2.22\;\mathrm{meV}$ (horizontal dashed line) is the result of
an effectively attractive non-linearity, i.e. $\alpha_{\text{eff}}<0$
estrecho_single-shot_2018 ; PhysRevB.89.235310 . Its physical origin is the
reduction in reservoir population at each condensate centre due to the inter-
condensate stimulated scattering leading to an effectively reduced exciton-
polariton interaction energy.

## References

  * (1) Mackey, M. C. & Glass, L.  Oscillation and chaos in physiological control systems.  _Science_ 197, 287–289 (1977).
  * (2) Oliveira, H. M. & Melo, L. V.  Huygens synchronization of two clocks.  _Scientific reports_ 5, 11548 (2015).
  * (3) Erneux, T.  _Applied Delay Differential Equations_ , vol. 3 (Springer-Verlag, New York, 2009).
  * (4) Atay, F. M.  _Complex Time-Delay Systems: Theory and Applications_ (Springer-Verlag, Berlin Heidelberg, 2010).
  * (5) Schuster, H. G. & Wagner, P.  Mutual Entrainment of Two Limit Cycle Oscillators with Time Delayed Coupling.  _Progress of Theoretical Physics_ 81, 939–945 (1989).
  * (6) Larger, L., Penkovsky, B. & Maistrenko, Y.  Virtual chimera states for delayed-feedback systems.  _Phys. Rev. Lett._ 111, 054103 (2013).
  * (7) Schöll, E. & Schuster, H. G.  _Handbook of Chaos Control_ , vol. 2 (Wiley-VCH, Weinheim, 2008).
  * (8) Uchida, A. _et al._ Fast physical random bit generation with chaotic semiconductor lasers.  _Nature Photonics_ 2, 728 (2008).
  * (9) Argyris, A. _et al._ Chaos-based communications at high bit rates using commercial fibre-optic links.  _Nature_ 438, 343 (2005).
  * (10) Van der Sande, G., Brunner, D. & Soriano, M. C.  Advances in photonic reservoir computing.  _Nanophotonics_ 6, 561–576 (2017).
  * (11) Appeltant, L. _et al._ Information processing using a single dynamical node as complex system.  _Nature Communications_ 2, 468 (2011).
  * (12) Larger, L. _et al._ Photonic information processing beyond turing: an optoelectronic implementation of reservoir computing.  _Optics express_ 20, 3241–3249 (2012).
  * (13) Brunner, D., Soriano, M. C., Mirasso, C. R. & Fischer, I.  Parallel photonic information processing at gigabyte per second data rates using transient states.  _Nature Communications_ 4, 1364 (2013).
  * (14) Vandoorne, K. _et al._ Experimental demonstration of reservoir computing on a silicon photonics chip.  _Nature Communications_ 5, 3541 (2014).
  * (15) Paquot, Y. _et al._ Optoelectronic reservoir computing.  _Scientific reports_ 2, 287 (2012).
  * (16) Martinenghi, R., Rybalko, S., Jacquot, M., Chembo, Y. K. & Larger, L.  Photonic nonlinear transient computing with multiple-delay wavelength dynamics.  _Phys. Rev. Lett._ 108, 244101 (2012).
  * (17) Larger, L. _et al._ High-speed photonic reservoir computing using a time-delay-based architecture: Million words per second classification.  _Phys. Rev. X_ 7, 011015 (2017).
  * (18) Kavokin, A., Baumberg, J. J., Malpuech, G. & Laussy, F. P.  _Microcavities_ (Oxford University Press, Oxford, 2007).
  * (19) Li, F. _et al._ Tunable Open-Access Microcavities for Solid-State Quantum Photonics and Polaritonics.  _Advanced Quantum Technologies_ 0, 1900060 (2019).
  * (20) Kasprzak, J. _et al._ Bose-einstein condensation of exciton polaritons.  _Nature_ 443, 409–414 (2006).
  * (21) Wu, Z. _et al._ Cubic-quintic condensate solitons in four-wave mixing.  _Physical Review A_ 88, 063828 (2013).
  * (22) Zhang, Z. _et al._ Particlelike Behavior of Topological Defects in Linear Wave Packets in Photonic Graphene.  _Physical Review Letters_ 122, 233905 (2019).
  * (23) Berloff, N. G. _et al._ Realizing the classical xy hamiltonian in polariton simulators.  _Nature Materials_ 16, 1120–1126 (2017).
  * (24) Bloch, I.  Ultracold quantum gases in optical lattices.  _Nature Physics_ 1, 23 (2005).
  * (25) Cataliotti, F. _et al._ Josephson junction arrays with bose-einstein condensates.  _Science_ 293, 843–846 (2001).
  * (26) Dung, D. _et al._ Variable potentials for thermalized light and coupled condensates.  _Nature Photonics_ 11, 565 (2017).
  * (27) Kim, N. Y. _et al._ Dynamical d-wave condensation of exciton–polaritons in a two-dimensional square-lattice potential.  _Nature Physics_ 7, 681 (2011).
  * (28) Abbarchi, M. _et al._ Macroscopic quantum self-trapping and josephson oscillations of exciton polaritons.  _Nature Physics_ 9, 275–279 (2013).
  * (29) Greiner, M., Mandel, O., Esslinger, T., Hänsch, T. W. & Bloch, I.  Quantum phase transition from a superfluid to a mott insulator in a gas of ultracold atoms.  _Nature_ 415, 39 (2002).
  * (30) Struck, J. _et al._ Quantum simulation of frustrated classical magnetism in triangular optical lattices.  _Science_ 333, 996–999 (2011).
  * (31) Zhang, Z. _et al._ Observation of parity-time symmetry in optically induced atomic lattices.  _Phys. Rev. Lett._ 117, 123601 (2016).
  * (32) Albiez, M. _et al._ Direct observation of tunneling and nonlinear self-trapping in a single bosonic josephson junction.  _Physical review letters_ 95, 010402 (2005).
  * (33) Soriano, M. C., García-Ojalvo, J., Mirasso, C. R. & Fischer, I.  Complex photonics: Dynamics and applications of delay-coupled semiconductors lasers.  _Rev. Mod. Phys._ 85, 421–470 (2013).
  * (34) Wouters, M., Carusotto, I. & Ciuti, C.  Spatial and spectral shape of inhomogeneous nonequilibrium exciton-polariton condensates.  _Phys. Rev. B_ 77, 115340 (2008).
  * (35) Christmann, G. _et al._ Oscillatory solitons and time-resolved phase locking of two polariton condensates.  _New Journal of Physics_ 16, 103039 (2014).
  * (36) Ohadi, H. _et al._ Nontrivial phase coupling in polariton multiplets.  _Phys. Rev. X_ 6, 031032 (2016).
  * (37) Cilibrizzi, P. _et al._ Polariton condensation in a strain-compensated planar microcavity with ingaas quantum wells.  _Applied Physics Letters_ 105, 191118 (2014).
  * (38) Prasad, A. _et al._ Universal occurrence of the phase-flip bifurcation in time-delay coupled systems.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 18, 023111 (2008).
  * (39) Adhikari, B. M., Prasad, A. & Dhamala, M.  Time-delay-induced phase-transition to synchrony in coupled bursting neurons.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 21, 023116 (2011).
  * (40) Barardi, A., Sancristóbal, B. & Garcia-Ojalvo, J.  Phase-coherence transitions and communication in the gamma range between delay-coupled neuronal populations.  _PLoS computational biology_ 10, e1003723 (2014).
  * (41) Dotson, N. M. & Gray, C. M.  Experimental observation of phase-flip transitions in the brain.  _Phys. Rev. E_ 94, 042420 (2016).
  * (42) Clerkin, E., O’Brien, S. & Amann, A.  Multistabilities and symmetry-broken one-color and two-color states in closely coupled single-mode lasers.  _Phys. Rev. E_ 89, 032919 (2014).
  * (43) Ramana Reddy, D. V., Sen, A. & Johnston, G. L.  Experimental evidence of time-delay-induced death in coupled limit-cycle oscillators.  _Phys. Rev. Lett._ 85, 3381–3384 (2000).
  * (44) Takamatsu, A., Fujii, T. & Endo, I.  Time delay effect in a living coupled oscillator system with the plasmodium of physarum polycephalum.  _Phys. Rev. Lett._ 85, 2026–2029 (2000).
  * (45) Cruz, J. M. _et al._ Phase-flip transition in coupled electrochemical cells.  _Phys. Rev. E_ 81, 046213 (2010).
  * (46) Manoj, K., Pawar, S. A. & Sujith, R. I.  Experimental evidence of amplitude death and phase-flip bifurcation between in-phase and anti-phase synchronization.  _Scientific Reports_ 8, 11626 (2018).
  * (47) Wouters, M. & Carusotto, I.  Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons.  _Phys. Rev. Lett._ 99, 140402 (2007).
  * (48) Asl, F. M. & Ulsoy, A. G.  Analysis of a system of linear delay differential equations.  _Journal of Dynamic Systems, Measurement, and Control_ 125, 215–223 (2003).
  * (49) Lang, R. & Kobayashi, K.  External optical feedback effects on semiconductor injection laser properties.  _IEEE Journal of Quantum Electronics_ 16, 347–355 (1980).
  * (50) Kozyreff, G., Vladimirov, A. & Mandel, P.  Global coupling with time delay in an array of semiconductor lasers.  _Physical Review Letters_ 85, 3809 (2000).
  * (51) Kalinin, K. P. & Berloff, N. G.  Polaritonic network as a paradigm for dynamics of coupled oscillators.  _arXiv_ arXiv:1902.09142 (2019).
  * (52) Töpfer, J. D., Sigurdsson, H., Pickup, L. & Lagoudakis, P. G.  _Data for Time-Delay Polaritonics_ (University of Southampton repository, 2019).  URL http://dx.doi.org/10.5258/SOTON/D1149.
  * (53) Lagoudakis, K. G., Pietka, B., Wouters, M., André, R. & Deveaud-Plédran, B.  Coherent oscillations in an exciton-polariton josephson junction.  _Phys. Rev. Lett._ 105, 120403 (2010).
  * (54) Lagoudakis, K. G. _et al._ Probing the dynamics of spontaneous quantum vortices in polariton superfluids.  _Phys. Rev. Lett._ 106, 115301 (2011).
  * (55) Wouters, M., Liew, T. C. H. & Savona, V.  Energy relaxation in one-dimensional polariton condensates.  _Phys. Rev. B_ 82, 245315 (2010).
  * (56) Estrecho, E. _et al._ Single-shot condensation of exciton polaritons and the hole burning effect.  _Nature Communications_ 9, 2944 (2018).
  * (57) Smirnov, L. A., Smirnova, D. A., Ostrovskaya, E. A. & Kivshar, Y. S.  Dynamics and stability of dark solitons in exciton-polariton condensates.  _Phys. Rev. B_ 89, 235310 (2014).