Time-Delay Polaritonics

TL;DR

This paper explores a new regime of polariton condensates acting as non-linear oscillators with time-delayed interactions, enabling control over their dynamics and potential for extended configurations.

Contribution

It demonstrates the use of microcavity exciton-polariton condensates as macroscopic, optically controllable non-linear oscillators with time-delayed coupling, modeled by equations similar to the Lang-Kobayashi equation.

Findings

Demonstrated macroscopic polariton condensates as non-linear oscillators

Achieved deterministic tuning between fixed point and limit cycle regimes

Modeled the system with time-delayed coupled equations similar to Lang-Kobayashi

Abstract

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.

Equations

$${\Psi }_1(t)$$

$${\Psi }_2(t)$$

$${\mathcal{V}}_{11}=\left|\cos (\Delta \bar{\tau }/2)\right|,$$

$${\mathcal{V}}_{12}=\left|\sin (\Delta \bar{\tau }/2)\right|.$$

$$E\Psi (x)=\left(-\frac{{\hslash }^2{\partial }_x^2}{2m}+V(x)-i\frac{\hslash {\gamma }_{\mathrm{c}}}{2}\right)\Psi (x),$$

$$\Psi (x)=\{\begin{array}{ll}A{e}^{-ikx},& x\le -d/2\\ B{e}^{ikx}+C{e}^{-ikx},& |x|<d/2\\ D{e}^{ikx},& x\ge d/2\end{array}$$

$${\left[\frac{V_0}{\frac{i{\hslash }^{2}k}{m}-V_0}{e}^{ikd}\right]}^2=1.$$

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

$$E_{n,\pm }=\frac{{\hslash }^2{k}_{n,\pm }^2}{2m}-\frac{i\hslash {\gamma }_{\mathrm{c}}}{2},$$

$$\Psi (x,t)=c_1(t){\psi }_1(x)+c_2(t){\psi }_2(x).$$

$$i\hslash {\dot{c}}_i=\left[\frac{{\hslash }^2{k}^2}{2m}-i\frac{\hslash {\gamma }_{\mathrm{c}}}{2}+\kappa \left(V_0-\frac{i{\hslash }^{2}k}{m}\right)\right]c_i+V_0\kappa e^{ikd}{c}_j.$$

$$e^{ikd}{c}_j(t)=e^{-\kappa d}{c}_j(t-\tau ).$$

$$i{\dot{c}}_i$$

$${\dot{n}}_i$$

$$J(d)=J_0\left|H_0^{(1)}(k_{0}d)\right|,$$

$$i{\dot{c}}_i=\left(\Omega +({\alpha }_{\text{eff}}-i\sigma )|c_i{|}^2\right)c_i+J{e}^{i\beta }{c}_j(t-d/v).$$

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

$$\frac{\partial n_{\mathrm{A}}}{\partial t}$$

$$\frac{\partial n_{\mathrm{I}}}{\partial t}$$

$$i\hslash \frac{d\Psi (x)}{dt}=\left(-\frac{{\hslash }^2{\partial }_x^2}{2m}+V(x)-i\frac{\hslash {\gamma }_{\mathrm{c}}}{2}\right)\Psi (x),$$

$$i\hslash {\dot{c}}_i+i\hslash {\mathcal{N}}_{12}{\dot{c}}_j=F(c_i,c_j),$$

$${\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),$$

$$F(c_i,c_j)$$

$$+\left(\left[V_0-\frac{i{\hslash }^{2}k}{m}\right]c_j+V_0{e}^{ikd}{c}_i\right)\kappa e^{-i{k}^{\ast }d}.$$

$$\hat{M}=\left(\begin{array}{cc}1& {\mathcal{N}}_{12}\\ {\mathcal{N}}_{12}& 1\end{array}\right),$$

$$i\hslash \left(\begin{array}{c}{\dot{c}}_1\\ {\dot{c}}_2\end{array}\right)={\hat{M}}^{-1}\left(\begin{array}{c}F(c_1,c_2)\\ F(c_2,c_1)\end{array}\right).$$

$${\left[\frac{V_0}{\frac{i{\hslash }^{2}k}{m}-V_0}{e}^{ikd}\right]}^2=1,$$

$$i\hslash {\dot{c}}_i=\left[\frac{{\hslash }^2{k}^2}{2m}-i\frac{\hslash {\gamma }_{\mathrm{c}}}{2}+\kappa \left(V_0-\frac{i{\hslash }^{2}k}{m}\right)\right]c_i+V_0\kappa e^{ikd}{c}_j.$$

$$i\hslash {\dot{c}}_i=\left[\frac{{\hslash }^2{k}^2}{2m}-i\frac{\hslash {\gamma }_{\mathrm{c}}}{2}+\kappa \left(V_0-\frac{i{\hslash }^{2}k}{m}\right)\right]c_i+\left(V_0\cos (k_{\mathrm{c}}d)-\frac{{\hslash }^{2}k}{m}\sin (k_{\mathrm{c}}d)\right)e^{-\kappa d}\kappa c_j.$$

$$i{\dot{c}}_i$$

$${\dot{n}}_i$$

Full text

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: J.D.Toepfer@soton.ac.uk, Pavlos.Lagoudakis@soton.ac.uk

(March 15, 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

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

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 .

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

[TABLE]

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

[TABLE]

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.

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:

[TABLE]

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,

[TABLE]

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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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,

[TABLE]

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),

[TABLE]

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,

[TABLE]

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 ,

[TABLE]

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 [math]-order Hankel function of the first kind describes the cylindrically symmetric radial outflow of particles in the linear regime Wouters_2008PRB , therefore we choose

[TABLE]

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).

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,

[TABLE]

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

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 .

[TABLE]

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],

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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

[TABLE]

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:

[TABLE]

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.