Competing role of interactions in synchronization of exciton-polariton condensates

We present a theoretical study of synchronization dynamics in incoherently pumped exciton-polariton condensates in coupled polariton traps. Our analysis is based on a coupled-mode theory for the generalized Gross-Pitaevskii equation, which employs an expansion in non-Hermitian, pump-dependent modes appropriate for the pumped geometry. We find that polariton-polariton and reservoir-polariton interactions play competing roles and lead to qualitatively different synchronized phases of the coupled polariton modes as pumping power is increased. Crucially, these interactions can also act against each other to hinder synchronization. We map out a phase diagram and discuss the general characteristics of these phases using a generalized Adler equation.

Fundamentally, we study the synchronization of two coupled 'oscillators'; however, the unique platform of excitonpolariton condensates under incoherent pumping modifies this picture significantly. In these systems, an uncondensed fraction of polaritons (the "reservoir") is deposited by the pump, typically at high energies; interactions amongst them give rise to stimulated scattering towards lower energy states, which can lead to condensation at a threshold pump strength P L 1 . The condensate mode above this threshold power appears with a self-organized frequency and associated spatial pattern. Typically, above a second threshold P L 2 a second mode condenses with generally a different oscillation frequency [28,29]. While this threshold physics is similar to the photon laser, exciton-polaritons distinguish themselves with strong, pump-dependent nonlinear interactions, which come in two varieties. Quasi-particles belonging to the condensate interact, giving rise to purely energetic effects. On the other hand, the deposited reservoir polaritons provide both a repulsive potential and the source of gain that allows condensation in the first place. Crucially, these interactions strongly depend on quasi-particle density, and hence pump power, leading to an effective pump-dependent coupling between the oscillators that is different from Huygens' original clock oscillator model of synchronization [30,31]. The interactions mediate frequency modulation, which we find plays a crucial role in eventually locking the two oscillators to a single frequency at the synchronization threshold. In fact, over the same pumping range, synchronization is no longer observed if interactions are turned off. Strong interactions thus enable the synchronization threshold to be reached even for far detuned modes at low powers; this is different from the laser, where relatively weaker interactions fail to lower the threshold for synchronization below that of other instabilities (e.g. a third mode turning on), except under carefully constructed situations that maximise modal coupling and minimize detuning [32]).
This simple emergence of synchronization beyond a threshold power is not the full story, however; the different nature and origin of the two distinct nonlinear interactions, a feature quite unique to polariton condensates, reveals even richer dynamics. When one interaction is dominant relative to the other, synchronization is indeed aided by the interactionmediated frequency modulation. Crucially, the properties of the emergent synchronized state depend strongly on the nature of the dominant interaction. However, we also find that these strong interactions do not always help the two oscillators to synchronize. When both interactions have comparable effects -in a quantifiable way that we derive -the situation is reversed: the interactions play competing roles that actually prevent a synchronized state from being reached.
To efficiently capture the nonlinear, pump-dependent interactions and reservoir dynamics present in this system, we employ a temporal coupled mode theory (TCMT) of condensate dynamics introduced in our recent work [33]. This approach is built upon the standard description of polariton condensation via the generalized Gross-Pitaevskii equation (gGPE), but provides a more economical framework for both numerical and analytical studies. In particular, the modal theory enables us to develop an analytical model for the role of interactions based on a multivariable generalized Adler equation, and to numerically study the long-time limit of dynamics for a wide parameter regime. Our key result is a phase diagram which maps out the presence of distinct synchronized phases when either of the two interactions is dominant, and a desynchronized phase when both interactions are actively competing.
The rest of this paper is organized as follows: in Section I, we briefly review the non-Hermitian coupled-mode theory that is the foundation for the analysis to follow. In Section II, we introduce the model system for our study of synchronization: two coupled, detuned polariton traps under uniform incoherent pumping. The coupled-mode theory is specialized to this case, and an amplitude-phase-intensity basis is introduced to express the modal equations more efficiently. A generalized Adler equation for the phase dynamics is analyzed in Section III, which provides intuitive insight into the role played by interactions in the existence of a synchronized phase. Finally, in Section IV, the emergence of stable synchronized fixed points as a function of pump power is investigated, and a phase diagram is mapped out. Analytical results are compared both with full numerical simulations of the TCMT, as well as select simulations of the gGPE using a split-step integrator; excellent agreement is found amongst the various methods.

I. REVIEW OF NON-HERMITIAN COUPLED MODE THEORY
A. Generalized Gross-Pitaevskii Equation The dynamics of exciton-polariton condensation under incoherent pumping has been very effectively described via a generalized Gross-Pitaevskii equation (gGPE). This takes the form of a complex dynamical equation for the order parameter Ψ(r, t), coupled to a rate equation for the incoherent reservoir population n R deposited by the pump (we set = 1) [34]: Here, V(r) is an arbitrary confining potential; for our study of synchronization, V(r) describes two tunnel-coupled and detuned polariton traps. Such models are relevant to condensation in disordered semiconductor microcavities, as well as recently realized coupled micropillar systems [35,36]. The pump depositing the exciton reservoir has strength P and spatial profile f (r). We consider condensation under uniform incoherent pumping; the confining and pumping geometries are depicted in Fig. 1 (a). Scattering from the exciton reservoir provides gain (rate R) leading to the formation of a coherent condensate when polariton losses γ c are overcome. The reservoir relaxation rate γ R encapsulates all loss processes other than scattering into the condensate. Most important for the purpose of this work are the two distinct types of interparticle interactions present: between polaritons within the condensate (strength g), and between reservoir excitons and condensate polaritons (strength g R ). Within the mean-field description employed here, interactions provide repulsive potentials that lead to polaritons experiencing pump-dependent frequency blueshifts. Note that reservoir-mediated gain and blueshift arise from separate physical processes (scattering and repulsive interactions respectively), and are thus characterized by generally different parameters R and g R . Finally, an additional diffusion term for reservoir excitons may be included in Eq. (1b); however, this term is often neglected due to reservoir excitons having much heavier mass relative to that of condensate polaritons. We justify this approximation in Appendix G.
Unfortunately, the set of coupled nonlinear partial differential equations (PDEs) described by Eqs. (1a) are numerically expensive to solve and provide limited analytic insight except in some special cases [37]. However, from this starting point we are able to arrive at a modal description that proves much more efficient. This approach, originally developed in Ref. [33], projects Eqs. (1a), (1b) onto a spatial basis consisting of pump-dependent non-Hermitian modes, which we introduce next. The result of this projection is a set of ordinary differential equations for time-dependent basis expansion coefficients; these equations, expressed in Eqs. (7a), (7b), comprise our temporal coupled-mode theory (TCMT).

B. Non-Hermitian Pump modes
To define the spatial basis for the TCMT, we begin by considering a linearized version of the gGPE, in the absence of polariton-polariton interactions and pump depletion. We drop nonlinear terms ∝ |Ψ| 2 , and replace the reservoir density by its undepleted value, n R = P f (r)/γ R ; then Eqs (1a), (1b) reduce to the single linearized dynamical equation: Eq.
(2) provides an exact description of the physical situation before and just upto the formation of a condensate, but is clearly not equivalent to the full gGPE beyond the condensation threshold, where |Ψ| 2 = 0. For our purposes, however, Eq. (2) is used to define the linear non-Hermitian generator H L (P ), which incorporates the trapping potential V(r) and the pump-induced potential. Then, the pump power P is arbitrarily tunable, and serves to parametrize the pump dependence of this linear generator. The computational modes {ϕ n } we employ for our modal theory are then eigenmodes of this generator, defined as: H L (P )ϕ n (r; P ) = ν n (P )ϕ n (r; P ).
Since H L (P ) is non-Hermitian, its eigenmodes are not orthogonal relative to the standard inner product. However, it can be shown that by introducing a set of dual modes {φ n } which satisfy the dual non-Hermitian eigenproblem H * L (P )φ n = ν * n (P )φ n , a complete basis may be obtained, satisfying the biorthogonality relation [38,39]: where P is the minimal region beyond which the pump has vanishing strength. It is easy to see thatφ * n = ϕ n , so that Eq. (4) reduces to a self-orthogonality relation: P dr ϕ n (r; P )ϕ m (r; P ) = δ nm The linear dynamics described by Eq. (2) are encoded in the pump-dependent, complex eigenvalues ν n (P ) = ω n (P ) + iγ n (P ) of each eigenmode. The real part ω n (P ) represents the modal frequency, determined by the confining potential V(r) and the blueshift due to reservoir excitons (∝ g R ). On the other hand, the imaginary parts γ n (P ) of the eigenvalues characterize non-Hermitian pumping and dissipation: they describe the net gain experienced by the nth mode at a given pump power. For low enough pump powers, all modal eigenvalues have negative imaginary parts, indicative of a belowthreshold regime -all modes experience net loss. As the pump power increases, the modal eigenvalues flow in the complex plane as the gain from the pump increases. For a specific pump strength P = P L n , which we refer to as the linear threshold power, the nth eigenmode acquires a vanishing imaginary part; for P > P L n , this mode experiences gain [See Fig. 1 (c)]. A threshold pump power can be associated with each mode; it is determined by how effective the specific mode is at utilizing gain from the pump [33], in the presence of the confining and pump-induced potential. Modes that are spatially configured to utilize the pump more efficiently have lower threshold pump values, and vice versa.
Due to their explicit pump dependence and associated threshold physics, we refer to the non-Hermitian eigenmodes of the linear generator as pump modes from this point on. Furthermore, the linear threshold power value P L n associated with each pump mode {ϕ n } hints at a natural ordering principle for these non-Hermitian modes. This allows us to a define a suitable truncation scheme for an expansion in a basis of these pump modes, which we discuss next.

C. Coupled-Mode Equations
Having defined our computational basis, we are now equipped to tackle the full nonlinear problem posed by Eqs. (1a), (1b). We proceed to expand the condensate order parameter Ψ(r, t) in the set of pump-dependent non-Hermitian pump modes with time-dependent coefficients: Ψ(r, t) = N n=1 a n (t)ϕ n (r; P ) where the sum may extend over an arbitrary number N of pump modes. The choice of modes to include in this expansion may seem unclear. Note that a truncation scheme that retains only pump modes with the lowest frequencies would not be appropriate in all cases: the non-Hermitian dynamics can lead to condensation in a mode that is not the lowest frequency mode for a given pumping potential. Instead, a truncation based on the linear threshold powers of the pump modes is more appropriate. We choose a basis of size N to comprise of the modes with the N lowest linear threshold power values. In this way, the modes that best optimize gain from the pump are retained in the basis expansion [33].
With the expansion placed on firm footing, we may now project the fully nonlinear Eqs. (1a), (1b) onto this basis, and integrate out the spatial dependence using Eq. (5). We leave the details of this procedure for Appendix A; the resulting dynamical equations corresponding to Eq. (1a) and Eq. (1b) re-spectively are given by: i d dt a n = ν n (P )a n + g In the above, we have introduced dynamical variables describing the evolving reservoir density, referred to as the reservoir matrix elements N nm (t). These are given by: and hence describe the change in reservoir density relative to its unsaturated value. Also, we define the spatial overlap matrix elements: Eqs. (7a), (7b) constitute the TCMT with nonlinear interactions and pump depletion. In Eq. (7a), we clearly see the eigenvalues ν n (P ) = ω n (P ) + iγ n (P ) controlling linear dynamics. All remaining terms incorporate the nonlinear effects neglected in the linear theory; the term ∝ g describes polariton-polariton interactions, with the mode overlap integrals A nmkq describing self-interaction and cross-mode interaction contributions. The final term with reservoir matrix elements N nm describes the effect of reservoir depletion; note that these reservoir matrix elements are themselves dynamically determined by Eq. (7b).
The TCMT provides an efficient spectral description of the nonlinear, non-Hermitian dynamics of polariton condensation, for arbitrary pumping and trapping potentials. In Ref. [33], TCMT simulations were compared with direct splitstep integration of the gGPE, and the results found to agree very favourably in a variety of settings. In the next section, we employ this TCMT for an analytical and numerical study of synchronization in coupled polariton traps.

A. System and Parameters
For the study of synchronization, we consider a confining potential V(r) that describes two coupled polariton traps, and a uniform incoherent pump across both traps [See Fig. 1 (a)]. Generally, there is a (low) pumping range wherein it is accurate to truncate the TCMT to only the two most preferred pump modes, namely the modes with lowest linear threshold powers. We consider specifically a detuned trap regime here, where these two modes are predominantly confined to the left or right trap; hence we denote the modes by ϕ l and ϕ r respectively, as shown in Fig. 1 (b). We emphasize again that these modes are not simply eigenmodes of the individual traps; rather they account for tunneling between the traps, the unsaturated pump, and dissipation through cavity loss. Under pumping with a wide, homogeneous pump spot, it is the static trapping potential that predominantly determines the modes ϕ l , ϕ r ; the non-Hermitian pump-induced potential does not modify the modes directly in this case, but our analysis does not make any explicit simplifications using this fact. However the non-Hermitian potential controls the evolution of modal eigenvalues with pump power: the homogeneous pump spot ensures this evolution is similar for both eigenvalues [see Fig. 1 (c)]. Furthermore, both modes have an approximately equal linear threshold power P L ; this is evident in how the eigenvalues cross the real line concurrently as pump power is increased. In the absence of pump-induced repulsion, the modal frequencies are ω l0 and ω r0 , determined by the trapping potential alone; we choose ω l0 > ω r0 , and define the bare pump mode detuning ∆ω ≡ ω l0 − ω r0 ; for the potential landscape V(r) chosen here, ∆ω 0 = 0.12 meV, a typical value consistent with modal detunings in a variety of polariton trapping structures [21,36]. For concreteness, we set R = 0.1 µm 2 meV, m −1 = 0.59 µm 2 meV, γ c = 1 meV, and γ R = 10 meV.
Note that we have chosen to operate in the regime of fast reservoir relaxation, γ R γ c , which is relevant for condensation in disordered semiconductor quantum wells. More importantly, in the opposite case where γ R γ c , we find a tendency towards strongly multimode behaviour, in agreement with recent numerical studies in this non-adiabatic regime [33,40,41]. The complicated dynamics can in fact be captured quite effectively by the TCMT, provided we include more than two modes in our basis [See Appendix F]. Therefore, to study dynamics within a two-mode approximation, we restrict ourselves to the fast reservoir relaxation regime.

B. Amplitude-Phase-Intensity basis
The condensate wavefunction in the two-mode TCMT is given by Eq. (6) for n = l, r, with the time-dependence entirely included in the mode coefficients {a n (t)}. For the study of synchronization, it proves useful to access the amplitude and phase dynamics directly; to this end we write the coefficients in a phasor representation, a n (t) =ā n (t)e −iφn(t) , n = l, r, and then cast Eqs. (7a) ,(7b) in terms of the dynamical variables {φ, z, ρ}: where φ is the relative phase, z the normalized modal intensity imbalance, and ρ the total intensity. Note that both φ ∈ [−π, π], and z ∈ [−1, 1] are bounded variables, while ρ generally increases monotonically with the pump power P ; these observations prove very useful in our analysis later. Also, in the regime of fast reservoir relaxation, the reservoir matrix elements N nm can be solved for in terms of the coefficients {a n }, so that the {φ, z, ρ} variables are sufficient for a full description of the reservoir-condensate system. A large body of earlier work on synchronization [42] concerns itself with weakly coupled Hermitian systems: the total intensity ρ is a conserved quantity and oscillator amplitudes are deemed to be approximately stationary so thatż ≈ 0, thereby allowing a phase-only description of the dynamics. In the desynchronized regime, the modes typically evolve at two distinct frequencies, namely the blueshifted trap frequencies, but modified by nonlinear interactions. The phases φ l,r (t) then exhibit drift-like evolution. On the other hand, in the synchronized regime the two modes begin oscillating at the same unique frequency Ω 0 . When this occurs, both phases φ l,r (t) exhibit a linear drift in time with the same drift constant, given by the synchronized frequency, Ω 0 t, plus (in general) a constant offset. Crucially, then, it follows that the relative phase φ = φ l − φ r becomes stationary in time,φ = 0. This requirement is an analytic signature of frequency synchronization in such phase-only systems.
For the non-Hermitian, pump-dependent system under consideration here, we find it crucial to keep track of z and ρ dynamics in addition to the relative phase evolution. Oscillations in z and ρ are pump-dependent and may be large [See Appendix D], so that the aforementioned approximations are not always valid. More importantly, while the relative phase may become stationary momentarily (φ = 0), we find that this state can only persist ifż =ρ = 0 concurrently, and if the state is stable to fluctuations in any of these dynamical quantities. The requirement of such stable fixed points in {φ, z, ρ} space is then the analytic signature of synchronization; comparisons between analytic and numerical results in Section IV B confirm that the modal frequency detuning does in fact vanish at such fixed points. Although the present system is not described by a phaseonly model, the phase dynamics still remain quite informative: sinceφ = 0 is still a necessary (though not sufficient) condition for synchronization to occur, it can place strong constraints on the synchronized phase, which we will explore next.

III. GENERALIZED ADLER EQUATION
One of the earliest models of synchronization involves the standard single-variable Adler equation [43]: The Adler equation describes the dynamics of the relative phase φ, of two coupled oscillators for example, in the presence of a detuning term Ω that causes φ to drift linearly in time, and a coupling term F sin(φ) that encourages its pinning to a constant value. Clearly, the synchronizedφ = 0 solution is possible only if −F < Ω < F , namely when the detuning term is small enough compared to the coupling. This model very successfully describes a broad range of synchronization and phase locking dynamics, from injection locking of oscillators to an external drive [44,45], to the synchronization of coupled oscillators within a phase-only picture [46]. In the present case, a somewhat more complicated equation for the relative phase φ may also be obtained by transforming Eqs. (7a), (7b) to {φ, z, ρ} space as defined in Eq. (10); we refer to it as the generalized Adler equation: In analogy with the standard Adler equation, we refer to the φ-independent term Ω(z, ρ) as the detuning term (generally not just the same as the bare pump mode detuning ∆ω 0 ), and ρF (φ, z) as the coupling term, which is periodic in φ. The ρ dependence of both terms reflects the pump-dependent nature of the nonlinear modal interaction. Note that Ω in Eq. (11) represents the bare frequency detuning of the uncoupled oscillators, and not the actual emergent detuning in the presence of coupling. Similarly, the detuning term Ω(z, ρ) is not the emergent detuning between the two coupled polariton modes; the latter is determined via simulation of the modal equations [See Section IV B]. The condition forφ = 0 imposed by the generalized Adler equation becomes: where max φ /min φ {f } are the maximum/minimum values respectively taken by f (φ, z, ρ) as φ varies in [−π, π], for a given (z, ρ). Eq. 13) implies that synchronization is again possible only when the detuning is small enough compared to the coupling term; now, however, these terms are no longer constants like in Eq. (11), but are rather configuration dependent. To understand the implications of this constraint, we study the detuning and coupling terms separately, beginning with the former.

A. Nonlinear Detuning modification
The detuning term Ω(z, ρ) includes the bare modal frequency difference and its modification due to nonlinear inter-actions. Here, we will see that this detuning term provides intuitive insight into the effect of nonlinear interactions on synchronization dynamics. The full expression for the detuning term is given by Eq. (B5) in Appendix B; in the main text, we find it more instructive to present the form Ω(z, ρ) takes for specific regimes. In the absence of any repulsive interactions, g, g R = 0, the detuning term is a constant equal to the bare pump mode detuning, independent of the pump power. The situation becomes more interesting once either type of repulsive interaction is active. If only interactions of polaritons within the condensate are considered -we define this as the g-mediated regime (g R = 0)the detuning term becomes: We have momentarily reverted to amplitude variablesā n since each term is intuitively clearest here. The first bracketed term represents the blueshift of the the left trap mode due to repulsion from polaritons occupying that mode (∝ A llllā 2 l ), and from polaritons occupying the right trap mode (∝ A llrrā 2 r ); the latter arises since the modes have nonzero spatial overlap. The second bracketed term is the corresponding blueshift of the right trap mode. Here A llrr A llll ≈ A rrll Arrrr 0.1, so that the 'direct' blueshift is more important, and causes the mode with higher occupation to be more strongly blueshifted due to g. A reduction in the detuning term -which aids synchronization, as per Eq. (13) -thus requires the low frequency mode (r) to experience a stronger blueshift, which is possible when this mode has higher occupation (ā 2 r >ā 2 l ), namely z < 0, as depicted in Fig. 2 (a).
In the g R -mediated regime where only interactions between reservoir excitons and condensate polaritons are active (g = 0), we have instead: The exciton repulsion ∝ g R now linearly blueshifts the mode frequencies, lending them the P dependence shown in Fig. 1 (c); the homogeneous pump spot ensures an equal blueshift that leaves the pump mode detuning unchanged, ω l (P ) − ω r (P ) = ∆ω 0 . However, condensate formation depletes the exciton reservoir, since the latter serves as the source of gain. In particular, the higher the occupation of a mode, the more depleted is the exciton reservoir it sees. This reservoir depletion reduces the reservoir-mediated blueshift experienced by each mode; in the detuning term, this reduction is captured by the N nn elements, which measure the reservoir depletion seen by mode n, and are explicitly negative [c.f. Eq. (8)]. To reduce the detuning term in this case, the low frequency mode (r) again needs to be more strongly 2. (a) Schematic variation of detuning term Ω with z in the gmediated regime (middle) and gR mediated regime (bottom). In the g-mediated regime, the frequency blueshift comes from polaritons within the condensate. Hence the effective modal detuning term Ω is lowered when the low frequency mode has higher occupation (z < 0) and experiences a stronger blueshift. In the gR-mediated regime, the blueshift is now due to reservoir excitons, and Ω is instead reduced when the low frequency mode has lower occupation (z > 0) and sees a less depleted exciton reservoir. Right panel: Evolution of detuning (solid blue/red) and coupling terms (dashed black), scaled by ∆ω0, with ρ in the (b) g-mediated regime (gR = 0), and (c) the gR-mediated regime (g = 0). Synchronization is possible only in the shaded regions. These regions prefer z < 0 in the g-mediated regime and z > 0 in the gR-mediated regime, and grow with ρ.
blueshifted; however, due to the different blueshift mechanism, the lower frequency mode (r) must now have a lower occupation (ā 2 r <ā 2 l ), namely z > 0, since it then sees a less depleted exciton reservoir and can experience a stronger blueshift. This is exactly the opposite to the g-mediated case [see Fig. 2 (a)]. Hence, we see a configuration-dependence to the reduction of the detuning term in the generalized Adler equation, with the gand g R -mediated regimes preferring opposite configurations.
Finally, note the additional prefactor of γ c /γ R = 0.1 relative to the g-mediated case, which weakens the g R mediated blueshift. This can be explained as follows: for larger values of γ R , the pump threshold for condensation increases. The stronger pumping allows for a stronger scattering ∝ R|Ψ| 2 from the reservoir to the condensate before reservoir depletion becomes important, which in turn means a higher condensate occupation is possible (for fixed polariton loss rate γ c ). The saturated reservoir density, however, is unchanged, since the stronger scattering uses up the additional pump power. Therefore, reservoir-dependent terms are reduced by a factor of γ c /γ R relative to terms that depend on the condensate density alone.

B. Pump-dependent coupling
Moving on to the coupling term ρF (φ, z), our choice of variables immediately indicates its growth with pump power via the explicit scaling by ρ; the precise dependence is clarified later. More explicitly, one finds: a nonlinear function of z with multi-harmonic φ-dependence.
Note that the coupling term may also be divided into a gmediated term ∝ F g and a g R -mediated term ∝ F g R ; the latter is again weaker in this regime by the factor γ c /γ R . The explicit forms of these functions may be found in Appendix B.
Most importantly, we find the coupling term is approximately unchanged under z → −z, which is very different from the behaviour of the detuning term. Both detuning and coupling terms combine to determine the possibility of a synchronized phase via Eq. (13). This is best explored graphically; we first consider the g-mediated regime, where g R = 0. In Fig. 2 (b), we plot max φ /min φ {ρF (φ, z)} (dashed black) and Ω(z, ρ) (solid blue) for fixed ρ, over the entire range of z ∈ [−1, 1]. The value of ρ is proportional to the pump power: the left panel is for ρ = 10, while the right panel shows ρ = 40, corresponding to an increased pump power. Eq. (13) forφ = 0 is satisfied in the shaded regions; here the detuning term lies within the range of the coupling term as determined by the dashed curves. Hence, in the unshaded regions, synchronization is impossible; clearly, the z < 0 region is preferred for synchronization here, a result stemming from the reduced value of Ω for this configuration. With increasing pump power, the shaded region area grows, as the increasing coupling strength and detuning modification make synchronization easier. The analogous plot in the g R -mediated regime (g = 0) is shown in Fig. 2 (c), for ρ = 10 and ρ = 30. Here, the situation is effectively reversed: the different mechanism for frequency modification means that the detuning term (solid red) is reduced for z > 0 instead, whereas the coupling term (dashed black) is mostly unchanged. As a result, synchronization is preferred here for the z > 0 configuration.

A. Stability Maps
The intuitive description of the previous section places constraints on theφ = 0 state, but does not guarantee synchronization; a complete analysis requires studying the full system of equations given by: Here we have rewritten Eq. 12 to isolate the bare pump mode detuning, ∆ω 0 . G φ,z,ρ are functions of system parameters and the variables (φ, z) only. The explicit forms of these functions are provided in Appendix B; we will find that the expressions above can already yield useful insight. While Eqs. (18a)-(18c) cannot be analytically solved for the fixed points, progress can be made if both modal gains are taken to be equal, γ l (P ) ≈ γ r (P ). Recall that the γ n (P ) are pre-determined by solving the non-Hermitian problem for the pump modes; we find the modal gains are indeed numerically equal in the present case. Note that this is also the most interesting scenario, since neither mode is preferred over the other by the pump.
We consider fixing all system parameters other than the pump power. Theρ equation can then be solved to obtain the parametric dependence (on φ and z) of the fixed points of ρ, The dominant power dependence here comes from the evolution of γ n (P ), which is shown in Fig. 1 (c). The function G ρ (φ, z) in the denominator is independent of P , and hence for a given (φ, z) pair, the steady state value of ρ evolves linearly with pump power. In what follows, it is then justified to use ρ as a surrogate variable for the pump power P . Theż equation also simplifies toż = 0 = ρ G z (φ, z). This implies that the values of (φ, z) for whichż vanishes are unchanged with increasing ρ, and hence pump power.
In addition to the existence of fixed points in {φ, z, ρ} space, a persistent synchronized phase demands that such fixed points be dynamically stable. We find that a standard stability analysis of Eqs. (18a)-(18c) based on the eigenvalues λ of the associated Jacobian matrix J also benefits from the simplified ρ dependence. At the fixed points, we have: The matrix element ∝ ∆ω 0 /ρ 2 yields terms in the characteristic equation that are suppressed by 1/ρ relative to remaining contributions at the same order. For large enough ρ, this term may be dropped; then, the resulting Jacobian matrix has a characteristic equation χ(φ, z, λ/ρ) = 0, with a crucial implication: the ρ dependence of the characteristic equation serves only to scale its roots, the eigenvalues of J. This result is explicitly derived in Appendix C. Hence, the sign of any eigenvalue does not change as ρ, and therefore pump power, is changed. While this conclusion holds only approximately, and that too at the fixed points of Eqs. (18a)-(18c) at any pump power, we find that in practice the signs of eigenvalues of J are quite robust to changes in ρ, for values of ρ close to ρ FP at that pump power. Thus, for any fixed set of system parameters, stable regions where all eigenvalues of J have negative real parts, and unstable regions where at least one eigenvalue has positive real part, are approximately unchanged with pump power. Because φ-z space is bounded in [−π, π] × [−1, 1], the emergence and movement of fixed points with increasing pump power can be tracked on finite, fixed stability maps of stable and unstable regions to characterize the global behavior of the system (even though the entire space is not explored for any given initial condition). An example of a stability map is shown for the g R -mediated regime in Fig. 3 (a), with plain regions being stable and shaded regions unstable, and theż = 0 contours shown in solid red; all of these features are unchanged with pump power. Only theφ = 0 contours (dashed blue) change with P ; note that such a contour only exists at points where condition (13) is satisfied. Hence we can now see how the generalized Adler equation plays a defining role in the emergence of a synchronized phase with increasing pump power. In the g R -mediated regime, it is clear from Fig. 2 (c) that theφ = 0 contours prefer the z > 0 region. The trajectories of fixed points with P are shown in the lower panel of Fig. 3 (a). The fixed point in the z > 0 region flows with increasing P and enters the stable region; as this crossing occurs, a synchronized phase of the system becomes stable. Note that in the z < 0 region, theż = 0 contour exists in a stable region; however, a stable fixed point cannot emerge until theφ = 0 contour spreads in this region, which is clearly restricted based on our analysis of the generalized Adler equation [See Fig. 2 (c)]. In contrast, Fig. 3 (c) shows a stability map in the g-mediated regime. Here, theφ = 0 contours prefer instead the z < 0 region, again clear from Fig. 2 (b). As such, the fixed point that flows from an unstable to a stable region now has z < 0.
Most interestingly, an intermediate regime exists where the two interactions compete; the different scaling of detuning and coupling terms discussed earlier implies that for this regime, g ∼ g R (γ c /γ R ). A typical stability map here is shown in Fig. 3 (b). For the same range of pump strengths

B. Phase Diagram
The predictions of the previous sections manifest strikingly in a phase diagram obtained via full numerical integration of the TCMT. We simulate the dynamical mode coefficients a n (t) to long times until a steady state is reached, and then compute the modal detuning ∆ from the Fourier transform of mode coefficients, F{a n }. This procedure is repeated for varying pump strengths P and increasing values of the polariton-polariton interaction strength g, starting from g = 0; importantly, the reservoir-polariton interaction strength is kept fixed at g R = 0.2 µm 2 meV. A plot of the computed detuning ∆ in P -g space constitutes the phase diagram shown in Fig. 4. In the yellow shaded regions, a desynchronized solution persists (∆ = 0). Outside this region, the modal detuning van-ishes and the two modes synchronize. Three distinct dynamical regimes can be identified, which we will now describe. For weak enough values of g g R (γ c /γ R ), a synchronized state emerges beyond a threshold pump power. In this region, the typical pump dependence of the Fourier spectra F{a n } of mode coefficients is shown in Fig. 4 (a). The curves projected below the spectra depict the evolution of the mode frequencies as a function of pump power. Clearly, the modes experience a frequency blueshift due to interactions as the pump is increased (dotted lines indicate constant frequency). This continues until the threshold power is reached, beyond which a single frequency synchronized state emerges (indicated in solid black). In the synchronized case,φ =ż =ρ = 0 in the steady state; we thus superimpose the steady state population imbalance in this synchronized regime on Fig. 4. For weak g, the polariton configuration appears with z > 0 (shaded red), in what we label as the regime of g R -mediated synchronization.
With stronger g ∼ g R (γ c /γ R ), the predicted competition between the two types of frequency modulation effects does in fact arise, and the synchronized state disappears. Typical frequency spectra F{a n } in this regime, shown in Fig. 4 (b), still indicate a frequency blueshift due to interactions, but the modes remain detuned. In this interacting regime, then, dynamics are reminiscent of ac Josephson oscillations be- tween the coupled condensates [See Section IV C]. For values of g that are stronger still, a synchronized phase emerges once more, but now with z < 0 (shaded blue). In this gmediated regime where interactions are strongest, the frequency blueshift is most pronounced, as is clear from the Fourier spectra in Fig. 4 (c). Again, the blueshifted modes lock to a single frequency beyond a synchronization power threshold.
The threshold for synchronization in both the g R -mediated and g-mediated regimes is predicted very well by our analysis of fixed points withφ =ż =ρ = 0 moving on a fixed stability map. The stability map is computed at fixed P = 1.15P L , making use of our previous result that stable and unstable regions are approximately unchanged with pump power. The phase boundary computed via this analysis is shown in dashed purple in Fig. 4. The slight discrepancy near the phase boundary may be explained as follows. For pump powers above the dashed purple line, the synchronized solution becomes stable; however, in a narrow range of pump powers, the initial condition determines whether the system settles into the synchronized phase or remains desynchronized. We note that this bistability of synchronized and desynchronized solutions is similar to results found in another two-mode configuration by Borgh et. al. [47].
Interestingly, for the pump range studied here, a synchronized state does not emerge if both interactions are turned off. For this case where g = g R = 0, the Fourier spectra of mode coefficients as a function of pump power are shown in Fig. 4 (d). The modes experience no blueshift, as expected, and the modal detuning is unchanged with pump power. Therefore, while the presence of either a dominant polariton-polariton interaction or reservoir-polariton interaction is necessary for synchronization, the competition of both interactions actually hinders synchronization.
Finally, one may note that the value of g R used in Fig. 4 is strong relative to g; this is chosen to make clear the role played by reservoir-polariton interactions in synchronization. Phase diagrams for weaker g R values are included in Appendix E, Fig. 7. The g R -mediated synchronization region generally shrinks and can even disappear for weaker values of the reservoir-polariton interaction strength. However, the key element -the competing nature of the interactions -still remains: if for g R = 0, synchronization occurs at a given pump strength for a fixed value of g, then nonzero g R pushes up this pump strength for g-mediated synchronization, and for strong enough g R may even lead to g R -mediated synchronization instead.

C. Comparisons with SSI
Since the TCMT is derived from the generalized Gross-Pitaevksii equation, a direct comparison between TCMT simulations of the previous section and a symplectic split-step integration (SSI) of the gGPE [c.f. Eqs. (1a), (1b)] can be made. In the one-dimensional geometry under consideration, we find that the two-mode TCMT is between one to two orders of magnitude faster than the SSI at a given pump power, for an equivalent spatial resolution and final integration time. The efficiency is primarily due to the TCMT's avoidance of spatial integration at every time step via a modal expansion. This fact also saves the TCMT's computation times from an unfavourable scaling with spatial dimension. SSI computation times, on the other hand, scale exponentially with dimension d, as N d g log(N g ) for a spatial grid with (uniform) density N g [48], which should render the TCMT even more favourable in higher dimensions (not considered here).
Thus, instead of computing the full phase diagram using the SSI, we provide select comparisons of TCMT and SSI results; in particular this is done for distinct points along the horizontal dashed white line in Fig. 4, which stretches across the three distinct dynamical regimes. In the synchronized regimes, at positions labelled (a) and (c) on the phase diagram, the condensate density is stationary in the steady state as mentioned earlier. We thus compare |Ψ| 2 , shown in Figs. 5 (a), (c) respectively; TCMT results are in solid black, and SSI results in dashed red. The synchronized configurations in (a), (c) are distinguished from self-trapping in coupled condensates [35] within the φ-z plane: while frequency synchronization leads to a constant relative phase (and population imbalance) as discussed earlier, the self-trapped regime has a relative phase that drifts linearly in time, and population imbalance that is weakly oscillating [49].
In the desynchronized regime, the mode amplitudes contain multiple frequencies unlike the synchronized case [See Fig. 4 (b)]. The condensate density is then no longer station-ary in time. Here a plot of the polariton occupations in each trap, in Fig. 5 (b), indicates ac Josephson-like oscillations between the detuned traps, with a running phase [20,35,47,50,51]. The TCMT and SSI exhibit excellent agreement in all three cases, which helps validate the preceding results of the two-mode TCMT.
Lastly, the dynamics of the reservoir across the synchronized and desynchronized regimes may also be computed using both methods; these results are included in Appendix E for completeness, and are also found to agree well.

V. CONCLUSIONS AND OUTLOOK
In this paper we have investigated the synchronization of two detuned, coupled polariton modes as a function of pump power, with a particular focus on the role played by the different nonlinear interactions unique to pumped polaritonic systems. Our analysis is based on a temporal coupled-mode theory that allows both analytical insight and efficient numerical simulations. We find that the polariton configuration in the emergent synchronized phase is strongly dependent on the relative influence of polariton-polariton and reservoir-polariton interactions. Most interestingly, the two types of interactions can have competing effects that prevent the emergence of a synchronized phase altogether. This conclusion is verified against direct simulation of the generalized Gross-Pitaevksii equation and very good agreement is found.
A natural extension of this work is to larger coupled systems, such as polariton lattices [52]. In particular, in recently realized flat band condensation [36], polariton states are highly localized due to disorder, while still possessing a finite frequency dispersion. The original experiment was concerned with threshold physics under weak pumping, and thus the polariton-polariton nonlinearity played a negligible role. However under stronger pumping, the relative strengths of reservoir-polariton and polariton-polariton interactions could be crucial in determining whether or not a spatially extended, synchronized (single-frequency) condensate could be formed, and what spatial configuration such a condensate may take. For studies in modern polariton lattice structures, the work here can be extended to regimes where γ R γ c ; the regimes directly studied in this paper should still be relevant in disorder-generated polariton trap geometries [21,34]. More generally, our results uncover an additional role of reservoir-mediated interactions in condensate dynamics, adding to other effects studied recently including the dynamical reservoir regime, multimode dynamics, and instabilities [33,40,41].
follows, the overlap matrix elements A nmrs , B nmrs are as defined in Eq. (9), while the reservoir matrix elements are scaled to extract out the explicit dependence on ρ, N nm → ρ γc γ R N nm . We begin with the equation for the relative phase φ, The bare detuning ∆ω 0 = ω l (P ) − ω r (P ), as presented in the main text. The function G φ (φ, z) can be written as: Here, N φ (φ, z) contains the reservoir matrix elements: The function K φ (φ, z) describes polariton-polariton repulsion terms: To compare with the form of the generalized Adler equation of Section III, the full expression for the detuning term Ω(z, ρ) contains simply the φ-independent terms from Eq. (B1): Note that N nn (z) is the φ-independent part of the reservoir matrix element N nn . These parts will be made clear in due course. The coupling term of Eq. (17) contains all the φdependent terms from Eq. (B1): Here, N rl (φ, z) is the φ-dependent part of the matrix element N rl . The functions F g and F g R are defined as: Next, we move on to the dynamical equation for the modal intensity imbalance, z, defined in the main text as: Here, the function G z (φ, z) has the form: As before, N z (φ, z) contains the reservoir matrix elements: while the function K z (φ, z) describes polariton-polariton repulsion terms: Lastly, we move on to the dynamical equation forρ: In all of the above, note that the reservoir matrix elements N nm themselves are dynamical quantities that can be expressed in terms of the variables (φ, z). In the regime where γ R γ c , which is the case we consider in the main text, the reservoir dynamics can be adiabatically eliminated. This behaviour is inherited by the scaled reservoir matrix elements, which are then also bound to the condensate evolution such that the equations for {φ,ż,ρ} are sufficient to determine condensate dynamics. We write the reservoir matrix elements N nm in terms of a φ-dependent part N nm (φ, z) and a φindependent part N nm (z), such that: Then, for the diagonal reservoir matrix elements, we find: and: Finally, the 'off-diagonal' or 'coupling' reservoir matrix element N rl = N lr can be similarly expressed by defining: In all the above expressions, we have defined P = P fr P 0 , where P 0 = γ c γ R /R. Note that this 'off-diagonal' matrix element is proportional to the modal coupling at lowest order, while the 'diagonal' matrix elements do in fact contain a contribution that is independent of the modal coupling.   initial modal occupations. The transient dynamics therefore involve a growing condensate density due to pumping, before a steady state is eventually reached. In the desynchronized regime, steady state curves are traced out in φ-z space, with examples shown in Fig. 6 (a) in a typical case. Clearly, the large variation in z over a period indicates the complicated amplitude dynamics can not be assumed to be approximately static. Furthermore, with increasing pump power, larger amplitude oscillations for z are observed.
In the regime where a synchronized phase is possible, a stable fixed point exists as determined by the stability analysis described in the main text. Here, the system dynamically flows towards this fixed point in φ-z space. In the long time limit, the system localizes at this fixed point; this is shown for a typical case in Fig. 6 (b). We show results here by unwrapping the phase φ, for clarity. The exact details of the flow at a fixed pump power will depend on initial conditions, as expected.

Appendix E: Supplementary simulations and reservoir dynamics
In this section, we include additional simulation results to supplement those included in the main text. We begin by presenting phase diagrams in P -g space for values of g R weaker than the value used for Fig. 4 of the main text. In Fig. 7 (a), (b), and (c), phase diagrams are plotted for g R = (0.15, 0.1, 0) µm 2 meV respectively. Clearly, the g R mediated synchronized regime shrinks as g R becomes weaker. Note, however, that the reservoir-polariton interaction g R still competes with synchronization mediated by the polaritonpolariton interaction strength g; when g R is turned off, the gmediated synchronization region is larger, over the same range of values of g, than the situation when g R = 0.
In the main text we presented condensate dynamics; here, we briefly discuss the associated reservoir dynamics. In the synchronized regimes, corresponding to positions (a) and (c) in Fig. 4, the condensate has a single frequency. Both the condensate and reservoir densities become stationary in the long time limit; the latter is then obtained by solving for the steady state of Eq. (1b): Recall that in our modal description of condensate dynamics, the TCMT can be used to obtain |Ψ| 2 but generally not the full reservoir density; instead, the simulated reservoir variables are the matrix elements N nm (t) [c.f. Eq. (7b)]. However, in the single frequency case, n R (r) is entirely determined by the condensate density, as is clear from Eq. (E1). Thus the TCMT can be used to obtain n R (r) directly in this case. We obtain the reservoir density, scaled by P 0 L, where L is the extent of the pump, and P 0 = γ c γ R /R as introduced earlier, for positions (a) and (c) of Fig. 4. The results are shown in solid black in Fig. 8 (a) and (c) respectively, and corresponding SSI results are shown in dashed red; note the excellent agreement. For (a), the synchronized phase has z > 0, and so the higher frequency, left trap mode has higher occupation. This indicates the condensate density is higher in the left trap, and thus the reservoir must be more strongly depleted there. This is precisely what is observed. For (c), where z < 0 instead, the reservoir overlapping with the right trap is more depleted, which agrees well with the simulated results.
In the desynchronized regimes, the reservoir density has a time-dependent evolution. In this more general case, a comparison of reservoir dynamics simulated by the TCMT and SSI can be carried out via the reservoir matrix elements instead. Using the full reservoir density computed via the SSI, N nm (t) may be obtained using the basis modes and Eq. (7b). Comparisons of the real and imaginary parts of the matrix elements are shown in Fig. 8 (b) for the TCMT (solid lines) and SSI (dashed black). Again, excellent agreement is found. Finally, we mention that for γ R γ c , the reservoir evolution still adiabatically follows the condensate dynamics, and in this regime the TCMT can be used to compute the full reservoir density as well. To provide a comparison valid for more general cases, we have analyzed the reservoir matrix elements here instead.

Appendix F: Weak γR regime
The dynamics we consider in this paper, and the simplified stability analysis, hold for the case where γ R γ c , namely the regime where the reservoir dynamics may be adiabatically eliminated. In the opposite, dynamical reservoir regime, where γ R γ c , we have found that complex dynamical ef- fects may arise, as found in recent numerical studies [33,40]. In particular, strongly multimode behaviour emerges for high pumping powers and -more importantly for the physics considered here -for strong nonlinearities.
To illustrate further the kinds of multimode dynamics present in the dynamical reservoir regime, we present in Fig. 9 (a) plots of the polariton number (integrated condensate density) in the left trap as a function of time using the TCMT, and an exact solution using an SSI, as in the main text, Fig. 4. Here, we take γ R = 0.1γ c ; the remaining parameters are summarized in the figure caption. We find complex dynamical features that agree quite well between both methods. The emergence of these features is clarified via the modal description provided by the TCMT. In Fig. 9 (b), the inset shows a spec- tral decomposition of the condensate wavefunction, plotting the normalized modal weights, an a1 . Only peaks corresponding to the five modes with greatest spectral weight are shown; a total of eleven modes is needed to reach the agreement shown in Fig. 9 (a). The first two modes are simply the left and right trap modes, ϕ l and ϕ r respectively. The two modes with next greatest spectral weight, modes number 3 and 5, are depicted in Fig. 9 (b). Hence, in this regime where reservoir relaxation is slow relative to polariton loss, condensate dynamics can typically become very complicated, and may involve multiple interacting modes, as found here.

Appendix G: Neglecting exciton diffusion
The reservoir dynamical equation in Eq. (1a) typically includes an additional term incorporating the diffusion of excitons, where D is introduced as the dimensionless diffusion constant.
In this paper, we have neglected this term, owing to this diffusion constant being relatively small for excitons due to their heavy mass; in typical systems [34] that we are considering here, D ≈ 10 −3 , which corresponds to actual diffusion rates on the order of 10 cm/s 2 .
However, in the present case where reservoir depletion plays an important role in determining synchronization dynamics, one may reasonably ask whether even relatively weak diffusion could strongly affect the observed physics. To verify that neglecting the diffusion term is valid, we first consider a perturbative solution in D; the D = 0 reservoir dynamical equation can be solved when γ R γ c , i.e. in the regime where the reservoir adiabatically follows the condensate evolution, to yield: The superscript (0) indicates results computed with the diffusion term neglected (D = 0). If we now expand the full reservoir density for nonzero D as a power series in D, n R = n R to first order in D is easily obtained: In deriving the above, we have neglected the finite extent of the pump. Also, any modifications of |Ψ| 2 due to exciton diffusion will lead to contributions that are second order in D; hence these are neglected. Again, in the regime of adiabatic elimination, the above equation may be solved: The effect of this correction to the diffusionless value n  We compute the above quantity for a typical synchronized solution, and find v ≈ 0.01D. This is a very small quantity, for which we expect almost no change with the inclusion of the diffusion term.
To further confirm the results of the above analysis, we perform SSI simulations of the gGPE while retaining the diffusion term in Eq. (G1). We compute the numerical results for the three cases indicated in Fig. 4 (a)-(c), and compare with the results for D = 0. The results are computed for D = 10 −3 , and an order of magnitude stronger diffusion, D = 10 −2 . The resulting plots are included in Fig. 10 (a)-(c). We see a negligible difference between the D = 0 and D = 0 case; the bottom panel of Fig. 10 (b) shows a zoomed in version of the desynchronized dynamics in the top panel, highlighting the minute discrepancy. The D = 10 −3 case is barely distinguishable from D = 0, while the use of a stronger diffusion constant also affects the dynamics only marginally. Certainly, no qualitatively different behaviour is observed. These simulations support the omission of the diffusion term in the reservoir dynamics for the systems considered here.