Externally Controlled Lotka-Volterra Dynamics in a Linearly Polarized Polariton Fluid
EExternally Controlled Lotka-Volterra Dynamicsin a Linearly Polarized Polariton Fluid
Matthias Pukrop and Stefan Schumacher
1, 2 Department of Physics and CeOPP, Universit¨at Paderborn,Warburger Straße 100, 33098 Paderborn, Germany College of Optical Sciences, University of Arizona, Tucson, AZ 85721, USA (Dated: January 15, 2021)Spontaneous formation of transverse patterns is ubiquitous in nonlinear dynamical systems of allkinds. An aspect of particular interest is the active control of such patterns. In nonlinear optical sys-tems this can be used for all-optical switching with transistor-like performance, for example realizedwith polaritons in a planar quantum-well semiconductor microcavity. Here we focus on a specificconfiguration which takes advantage of the intricate polarization dependencies in the interactingoptically driven polariton system. Besides detailed numerical simulations of the coupled light-fieldexciton dynamics, in the present paper we focus on the derivation of a simplified population com-petition model giving detailed insight into the underlying mechanisms from a nonlinear dynamicalsystems perspective. We show that such a model takes the form of a generalized Lotka-Volterra sys-tem for two competing populations explicitly including a source term that enables external control.We present a comprehensive analysis both of the existence and stability of stationary states in theparameter space spanned by spatial anisotropy and external control strength. We also constructphase boundaries in non-trivial regions and characterize emerging bifurcations. The populationcompetition model reproduces all key features of the switching observed in full numerical simula-tions of the rather complex semiconductor system and at the same time is simple enough for a fullyanalytical understanding.
I. INTRODUCTION
The demand for integrated optoelectronic devices inoptical communication networks has resulted in an in-crease of research activities targeted at functional all-optical components. For example, a wide range of dif-ferent approaches has been proposed to realize efficientall-optical switches exploiting the nonlinear optical prop-erties of different material systems, including organicphotonic crystals [1], rubidium atomic-vapor cells [2],or GaAs semiconductor microcavities as in Refs. [3–5]and Refs. [6–8]. The latter utilize the optical controlof transverse optical patterns to achieve transistor-likeswitching performance [2, 9]. Spontaneous formation ofspatially extended patterns has been intensively studiedduring the past decades [10] with applications to differ-ent areas of science including formation of sand ripplesand desert sand dunes [11], animal coat patterns suchas zebra stripes [12, 13], or geographic patterns in para-sitic insect populations [14]. With applications to opticalswitching, however, the quest for efficient external con-trol of these patterns naturally arises, but in many caseshas not been explored in detail. In the present work weinvestigate a specific example of all-optical switching ofpolariton patterns in a semiconductor quantum-well mi-crocavity system. In contrast to our previous work [9] wegive a detailed analysis of the optical switching dynam-ics from a nonlinear dynamical systems perspective. Tothis end we derive a simplified mode competition modelthat governs the essentials of the system dynamics but issimple enough to fully characterize the possible station-ary solutions and phase-space singularities in the relevantparameter space. This simplified model has the very gen- eral mathematical form of a generalized Lotka-Volterrasystem, with the addition of an inhomogeneity for ex-ternal control. In the main part of the paper (cf. Sec.III), we present a complete steady state analysis of thisnonlinear system. The solution space we obtain in depen-dence of spatial anisotropy and external control strengthis of very general nature and may be similarly realized inother systems where external control of population com-petition is studied such as chemical reactions or in thelife sciences.As our specific example, here we study planar semicon-ductor microcavities with a strong coupling between thecavity field and the exciton polarization that gives rise tothe formation of exciton polaritons [15]. These quasipar-ticles consist of a photonic and an excitonic part and arecharacterized by a normal-mode splitting in the disper-sion relation, on the lower branch leading to long coher-ence times and strong nonlinear interactions. The latterare driven by four-wave mixing processes of coherent po-lariton fields which can also be interpreted as polariton-polariton scattering, mediated through exciton-excitonscattering. For certain excitation conditions, small spa-tially varying density fluctuations can experience hugegrowth in particular modes due to the intrinsic feedbackmechanisms driven by four-wave mixing. This causesspatially homogeneous density distributions to becomeunstable such that the system’s symmetry is sponta-neously broken. This results in the formation of station-ary patterns, directly observable in the far field emissionfrom the microcavity. Figure 1(a) shows the excitationgeometry used with the pump at normal incidence (zeroin-plane momentum) and finite off-axis ( k (cid:54) =0) signals dueto amplified fluctuations. Figure 1(b) shows the normal- a r X i v : . [ c ond - m a t . o t h e r] J a n FIG. 1. (Color online) (a) Sketch of a planar quantum-well semiconductor microcavity with continuous-wave pumpat normal incidence (on-axis) and optional off-axis controlbeam. (b) Corresponding polariton dispersion relation inthe normal-mode splitting regime. The dominant four-wavemixing process is indicated, driving the pairwise scatteringof pump-induced polaritons onto the lower polariton branchwith opposite in-plane momenta for signal and idler modes. mode splitting of the dispersion relation into a lower(LPB) and upper polariton branch (UPB) alongside thebare exciton and cavity dispersions. Phase-matching con-ditions determine the efficiencies of the different scatter-ing processes, and therefore determine the possible pat-tern geometry. For scalar (only one circular polarizationstate) polariton fields hexagon patterns are favored [16–18]. Extending the model with polarization dependencea complex interplay of the TE/TM cavity-mode split-ting in the linear regime and the spin-dependent exciton-exciton interactions in the nonlinear regime arises. Theresulting polarization-induced spatial anisotropy [19] de-termines the possible unstable modes and stable patternsformed. For linearly polarized pump excitation slightlyabove instability threshold, cross-linearly polarized res-onant modes parallel and perpendicular to the pump’spolarization plane constitute the basic instabilities [20],resulting in two-spot or four-spot patterns. Making useof this spatially anisotropic polariton-polariton scatter-ing a reversible optical switching for orthogonal two-spotpatterns can be realized that is triggered by a weak ex-ternal control beam [9]. The advance of taking the po-larization dependence into account, compared to scalarswitching schemes [6–8], is the explicit breaking of therotational symmetry, offering to control the orientationof the initial pattern via a linearly polarized pump and toprovide an automatic back-switching mechanism. It alsomakes background-free detection possible, since the pat-tern switching takes place in the cross-polarized channel.
II. ORTHOGONAL SWITCHING OFTWO-SPOT PATTERNS
A detailed numerical investigation of the present sys-tem for the following excitation conditions was alreadydiscussed in Ref. [9]. The system is excited (driven) byan x -linearly polarized continuous wave pump at nor-mal incidence with flat-top Gaussian shape in the QWplane and an intensity slightly above the off-axis insta- FIG. 2. (Color online) Schematic illustration of the dom-inant contributions to the in-plane scattering of polaritonsonto the elastic circles defined by the TE and TM modes ofthe lower polariton branch. The pump indicated in the centeris x -polarized and slightly above off-axis instability threshold.Parallel and perpendicular to the pump polarization direction,the scattering can occur in either TE or TM mode (left panel).For pumping spectrally well below the exciton resonance, thescattering with polarization orthogonal to the pump’s is pre-ferred (center panel). With spatial anisotropy (which is partlyinduced by the linear pump polarization) scattering along thepump’s polarization direction onto the TE mode can domi-nate (right panel). bility threshold. In this case spontaneous breaking ofspatial and polarization symmetry is observed. Signalsare formed by resonant polariton scattering onto the elas-tic circle defined by the polariton dispersion (cf. Fig.1). As can also be understood based on a linear sta-bility analysis [20, 21], the scattering occurs predomi-nantly in four spatial directions oriented orthogonal andparallel to the pump polarization plane, respectively. Ingeneral, in any spatial direction signals can form eitherin the TE or in the TM mode as illustrated in Fig.2. However, for the pump polaritons scattered to finite k through Coulomb interaction, due to the underlyingspin-dependent exciton-exciton interactions, the scatter-ing probability is higher for a polarization state orthogo-nal to the pump polarization state [21]. Therefore, in thespatial direction parallel to the pump polarization state,the scattered signals preferentially form in the TE mode,and in the TM mode for scattering orthogonal to thepump polarization state (cf. Fig. 2). Out of these foursignals, a stationary two-spot pattern (only one modepair with opposite in-plane momenta) can be preparedby introducing some anisotropy favouring one directionover the other. Alongside a slight polarization inducedspatial anisotropy that is due to a slightly higher densityof states for the TE modes [20], an additional anisotropycan be introduced by tilting the pump beam slightly awayfrom normal incidence (cf. Fig. 2). In the nonlinear sys-tem studied, at sufficiently high densities of the favoredtwo-spot pattern, cross-saturation processes will lead tothe extinction of the signals in the other orthogonal direc-tion. This allows us to single out a two-spot pattern T which will serve as the initial state in the switching pro-cess. Here, we start with the two-spot pattern oriented inthe x -direction (direction of intrinsic anisotropy). Thisselection process is schematically visualized in Fig. 2.If now a weak (compared to the pump intensity) y -polarized control beam with an in-plane momentum onthe elastic circle and spatially orthogonal to the initialpattern is applied, stimulated scattering leads to popula-tion revival of the corresponding two-spot pattern. Dueto the cross-saturation effect, the initial pattern is desta-bilized and finally switched off completely while the or-thogonal pattern reaches a steady state T (cf. Fig. 3).After the control beam has been switched off again, theanisotropy leads to re-emergence of the initial T patternwhile the T state vanishes again. If the control beam istoo weak so that complete switching may not be possible,the system will remain in a stationary four-spot patternstate F. Based on this scheme, in Ref. [9] transistor-like reversible switching was demonstrated including asystematic study of switching times, minimum controlpower needed, and achievable gain. This was done bynumerical simulations of the nonlinear set of equationsof motion governing the coherent coupled light field andexciton dynamics in the microcavity system in the two-dimensional QW plane in real space,i (cid:126) ∂ t E ± =( H c − i γ c ) E ± + H ± E ∓ − Ω p ± + E ± pump i (cid:126) ∂ t p ± =( (cid:15) e − i γ e ) p ± − Ω(1 − α psf | p ± | ) E ± + T ++ | p ± | p ± + T + − | p ∓ | p ± . (1)Here, the index ± denotes the different com-ponents in the circular polarization basis, H c = (cid:126) ω − (cid:126) ( m TM + m TE )( ∂ x + ∂ y ) is the cavityHamiltonian and (cid:15) e the flat exciton dispersion. γ c and γ e represent the photon and exciton lossrates, Ω describes the photon-exciton coupling, and H ± = − (cid:126) ( m TM − m TE )( ∂ x ∓ i ∂ y ) couples the two po-larization components due to TE/TM-splitting. Thecubic nonlinearities consist of a phase-space filling term α psf , repulsive interaction T ++ for excitons with parallelspins, and attractive interaction T + − for excitons withopposite spins.Based on Eqs. (1) with an x -linearly polarized contin-uous wave pump beam at k ≈ | E | , in k -space withoutcontrol beam such that the stationary pattern T formsin panel (a) and with the control beam on, switching tothe stationary pattern T in panel (b). The control beamhas the same frequency as the pump beam (cl. Fig. 1).System parameters used and details of the calculationsare given in Appendix A. Upon switching off the controlbeam the switching action is reversed and the patternreturns to its original state T . Switching times in therange of ∼
100 ps are achievable [9].The focus of the present paper is not on the full numer-ical simulations of Eqs. (1) and resulting switching per-formance. Rather, starting from Eqs. (1) we will derivea simplified population competition model for selectedmodes in k -space (details of the derivation are given inAppendix B) to provide further insight into the under-lying phase space singularities that dictate the globalbehavior and solution space of the nonlinear dynamicalsystem studied. In a similar fashion this approach was FIG. 3. (Color online) (a) Stationary T and (b) T state ofthe polariton pattern switch as obtained from the numericalsolution of Eqs. (1) for x -linearly polarized cw pump at k ≈ y -polarization component of | E | in k -space in arbitrary units (a) before and (b) after switching on acontinuous wave control beam. Switching off the control beamresults in reversal of the switching process such that aftersufficiently long times the system returns to the stationaryT pattern. previously applied to the switching between subsets of ahexagonal pattern [22]. We will systematically analyzethe existence and stability properties of possible steadystates in dependence of the strength of the different in-volved physical processes for a typical orthogonal switch-ing setup comparable to the one introduced above. Tothis end, we will construct phase boundaries in represen-tative regions of parameter space and characterize rele-vant bifurcations. This will lead to a general understand-ing of crucial parameter dependencies such as the ratiobetween the control beam strength and the anisotropy,and also of the coexistence of solutions in certain pa-rameter regions and hysteresis behavior. Based on thesimplified model, we will also be able to show that thepolariton dynamics and optical switching phenomena dis-cussed in the present paper can mathematically be un-derstood based on an generalized Lotka-Volterra modelincluding an external control parameter. III. POPULATION COMPETITION MODEL
The simplified population competition model discussedin the remainder of the present paper can be derived fromthe set of equations of motion in Eq. (1) as detailed inAppendix B. It reads˙ A = α A − β A − θ A A ˙ A = α A − β A − θ A A + S. (2)Here A and A are real-valued amplitudes of the twoelementary states of the system, i.e., the two orthogonalpolariton mode pairs in k -space. In case of stationary so-lutions the phases in the full model (1) become “locked”,allowing us to remove them as dynamical variables fromthe equations which is suitable for the following steadystate analysis. The six dimensionless real-valued pos- FIG. 4. (Color online) Illustration of the different polariton-polariton scattering processes in k -space for the two elemen-tary mode-pairs A and A in Eq. (2). The dashed circle rep-resents ring in k -space on which efficient scattering is possible. Green solid ( red open ) circles mark the incoming (outgoing)modes. Solid ( dashed ) arrows mark their corresponding mo-menta. The α i describe the growth process of the two modepairs (basic instability leading to pattern formation) and the β i and θ i represent the self- and cross-saturation process. itive parameters, α i , β i , θ i are directly related to themain (up to third order) polariton-polariton scatteringprocesses of the original system as illustrated in Fig. 4.They can be calculated from the physical parameters ofthe full model (cf. Appendix C). These parameters areintrinsically different for A and A due to the polar-ization dependence and anisotropy. The linear processrepresenting growth of the resonant modes is describedby α i . Saturation processes are represented by the cu-bic terms which can be divided into self-saturation ( β i )and cross-saturation ( θ i ). The external control is de-scribed by the inhomogeneity S . If we rewrite Eqs. (2)as ( ˙ A , ˙ A ) T ≡ ( f , f ) T = f , steady states are charac-terized by f = . Four qualitatively different stationarysolutions are possible: i) two-spot pattern T with A (cid:54) =0and A =0, ii) two-spot pattern T with A =0 and A (cid:54) =0,iii) four-spot pattern F with A (cid:54) =0 and A (cid:54) =0, and iv) thetrivial solution with A =0 and A =0. Here we are onlyinterested in physical solutions A i ≥
0, and therefore thestate space is R ≥ , i.e. the first quadrant of the ( A , A )-plane. The linear term ( α ) leads to exponential growthof the corresponding mode pair and the self-saturation( β ) has a stabilizing effect. For these two processes, theequations (2) are automatically decoupled and may besolved separately, resulting in a stable four-spot pattern.However, the cross-saturation terms ( θ ) couple the twomodes and tend to suppress a particular mode pair, fa-voring the other mode pair, resulting in a two-spot pat-tern. Additionally, the external control acts as a sourceterm for A . The PC model thus describes the dynamicalcompetition between the three types of possible station-ary patterns: T , T and F. A phase in the PC modelis defined as a set consisting of the number of steadystates and their stability properties. These sets, i.e. the phase, are functions of the seven parameters and a phaseboundary in parameter space indicates the change in thenumber of steady states and/or their stability. This canonly happen at points where at least one eigenvalue ofthe corresponding Jacobian matrix J ≡ (cid:0) ∂ A j f i (cid:1) is zero[23], which is equivalent to the condition det J =0. Wethus find phase boundaries in parameter space at pointssatisfying { f , det J } = . A. Homogeneous Case S = 0 The homogeneous case of Eq. (2) ( S = 0) has theform of a generalized Lotka-Volterra (GLV) model [24]with cubic nonlinearities. The transformation A i → A i ≡ ˜ A i yields the usual Lotka-Volterra (LV) equations[25] with quadratic nonlinearities while conserving the( A , A ) phase space structure in the positive quadrant[26]. Hence, for S = 0 the phase portraits of (2) are topo-logically equivalent to those of the following well-knownsystem ∂ t ˜ A i = ˜ A i r i − (cid:88) j =1 c ij ˜ A j (3)with i = 1 ,
2, growth rate vector r = 2( α , α ), and com-munity matrix C = 2 (cid:18) β θ θ β (cid:19) , (4)which describes self- and cross-interaction. This LVmodel (3) for interspecific competition has been studiedin many different contexts, e.g. ecology [27], chemistry[28], economics [29], physics [30] and it was shown [25, 31]that its dynamical behavior is limited to three casesdepending on the parameters: i) Coexistence regimefor sgn(det C )=+1 (larger self-saturation), ii) bistabil-ity regime for sgn(det C )= − S = 0 the systemis solvable analytically we obtain explicit expressions forall steady states and phase boundaries: i) T : A = (cid:113) α β , A =0, stable for α α > β θ , ii) T : A =0, A = (cid:113) α β , stablefor α α < θ β , iii) F: A = (cid:113) α θ − α β θ θ − β β , A = (cid:113) α θ − α β θ θ − β β ,only exists for β θ ≶ α α ≶ θ β and is stable for θ θ <β β , iv)trivial solution: A = A =0. We already see that the four-spot solution does not exist in the entire S =0 plane incontrast to the two-spot pattern. The trivial solution also FIG. 5. (Color online) Phase boundaries in ( δα, S ) param-eter space showing regions of different stable and unstablesteady states in dependence of the anisotropy and externalcontrol. (a) Larger self-saturation β β >θ θ . (b) Largercross-saturation β β <θ θ . Green ( dark gray ) / red ( lightgray ) letters mark stable / unstable steady states. Lines S ∗ and S ∗ show the phase boundaries (5) and (6). exists everywhere but is always unstable since the eigen-values of its Jacobian α , α are positive. This solutionwill not be listed hereinafter. For a systematic discussionwe introduce a variable anisotropy parameter δα in thefirst equation of (2), α → α + δα , and set α = α =1. Itfavors A / for δα ≷
0, respectively. The resulting phaseboundaries are shown in Fig. 5 where the homogeneouscase ( S =0) is included in the extended region at the bot-tom of each plot. It shows the structure of the usual LVmodel. Green ( dark gray ) / red ( light gray ) letters markstable / unstable steady states. For large anisotropy | δα | only stable and unstable two-spot patterns are possible.This dominance regime does no longer depend on theinteraction parameters β i and θ i . The middle regionon the other hand is divided into two cases: A stablefour-spot pattern for larger self-saturation ( β β >θ θ ,coexistence, upper panel) and two simultaneously stable FIG. 6. (Color online) Phase plane flow for cases correspond-ing to Fig. 5. Black (white) dots mark stable (unstable)steady states. Black lines show representative orbits. Thethick black line in 3b and 6 is the stable manifold of the sad-dle point which separates two basins of attraction. two-spot patterns for larger cross-saturation ( β β <θ θ ,bistability, lower panel).The first four cases in Fig. 6 show the correspondingflow given by representative trajectories in state space forthe homogeneous case. The system’s dynamics are un-ambiguous in cases 1, 2, and 3a where only one attractorexists which determines the long time behavior. How-ever, in case 3b two attractors exist and the dynamicsnow depend on the system’s history leading to a hys-teresis effect. The two basins of attraction are markedwith different colors blue ( dark gray ) and orange ( lightgray ). They are separated by the stable manifold of thesaddle point F (often called separatrix ), defined by theset of points ( A , A ) which satisfy ( A , A ) → F for t → ∞ . In contrast to this, the line connecting all threesteady states is the unstable manifold of F, consisting ofpoints which satisfy ( A , A ) → F for t → −∞ . Anyinitial point will stay in its region ( blue or orange ) andend up at the corresponding attractor (T or T ). Thesystem thus shows hysteresis behavior which might pre-vent complete back-switching, and therefore should beavoided for switching purposes, for example by increas-ing the anisotropy.The Hartman Grobman theorem [23] ensures that thebehavior near hyperbolic equilibrium points is completelydetermined by its linearization. This no longer holds fornon-hyperbolic fixed points which are characterized bythe existence of at least one zero real-part eigenvalue ofthe corresponding Jacobian. Therefore, the behavior atthe phase boundaries in parameter space can not be an-alyzed via linear stability analysis. Instead one can usecenter manifold theory and normal forms [23] in order todetermine the occurring bifurcations. They are transcrit-ical for the usual LV model but different in the case ofcubic nonlinearities due to the additional Z × Z sym-metry. For both two-spot solutions symmetric pitchfork(PF) bifurcations occur at the phase boundaries whenthe anisotropy parameter is changed. They are super-critical for β β >θ θ (larger self-saturation) leading tothe coexistence regime with a stable F pattern and sub-critical for β β <θ θ (larger cross-saturation) leadingto the bistability regime with stable T and T pattern.Pitchfork bifurcations are typical for dynamical systemswith inversion symmetry, here A i → − A i . Since we areonly interested in positive solutions in the first quadrantof the phase plane, we either observe a transition from astable T solution to an unstable T and a stable F solu-tion (supercritical) or the same transition with reversedstability (subcritical).In conclusion, all phase portraits of the homogeneouscase can be completely reduced to the results of the well-known Lotka-Volterra model for two competitive species.Due to the additional inversion symmetry we obtainpitchfork bifurcations at the phase boundaries. The ho-mogeneous case is crucial for the polariton switch be-cause it describes the initial pattern formation and back-switching process in absence of the control beam. Bothcan only work reliably in regions where only the T pat-tern exists as a single attractor. Therefore, a sufficientminimum anisotropy in favor of A is needed. We havealso seen that if the cross-saturation dominates over theself-saturation, the coexistence of A and A is desta-bilized resulting in the extinction of A i and survivalof A i +1 . A strong interspecific competition thus pre-vents the coexistence. Otherwise, if the self-saturation is stronger than the cross-saturation, the dominating in-traspecific competition promotes coexistence. The pa-rameters used for the full numerical simulations in Sec.II lead to the larger self-saturation case and sufficientanisotropy (cf. Appendix C), providing reliable initialpattern formation and back-switching, since the F pat-tern is unstable. B. Inhomogeneous Case
S > The inhomogeneous ( S>
0) case of the PC model (2)describes the actual switching process induced by an ex-ternal control beam. This term can also be motivatedfor other systems where the GLV model is commonlyused, e.g. to include constant migration/harvesting inthe description of ecological systems or constant influxin chemical reactions. Therefore, it can be interpreted asan extension of the generalized Lotka-Volterra dynam-ics and the following analysis is of very general interestbut has not been investigated before. The influence ofconstant terms in the usual LV model with quadraticnonlinearities was investigated in Ref. [32] from a purelymathematical perspective, but in general inhomogeneouspopulation competition models did not receive much at-tention in the past. Here, we analyze and apply the GLVmodel with a constant inhomogeneity to the case of theorthogonal switching of two-spot polariton patterns.Also including inhomogeneity, all steady states andphase boundaries can still be determined analytically dueto the system’s simplicity. However, we note that in gen-eral solving multivariate polynomial equation systems isa difficult task which can be simplified using algebraicmethods (e.g. Gr¨obner Basis [33, 34]), as was shown inRef. [22]. The non-zero source term S in (2) breaks theinversion symmetry for A and prevents a T solution.This leaves us with only two competing patterns (calledT & F ) this time around. However, as we will see belowit includes the possibility of an additional qualitativelydifferent four-spot pattern solution F . The control pa-rameter S acts as a constant source term for the A modepair, leading to linear growth. Supporting A means si-multaneously suppressing A due to the cross-saturationeffect, leading to a more asymmetric F solution in fa-vor of A . A second four-spot pattern, F , which is infavor of A replaces T . The phase boundaries are givenexplicitly by S ∗ ( δα ) = β (cid:114) δαθ − (cid:114) δαθ , (5)and additionally, for the case of larger cross-saturation( β β <θ θ ) S ∗ ( δα ) = 6 (cid:115) δα ) θ − β ] β ( θ θ − β β ) . (6)The explicit expression for the steady states T , F ,and F are not given here due to their excessive length. FIG. 7. (Color online) Evolution of steady states for A and A in parameter space. Green (s) / red (u) surfaces belongto stable / unstable solutions. Black lines correspond to thephase boundaries projected on the ( δα, S ) plane (a) Largerself-saturation: one of the two pitchfork bifurcations on the S =0 line remains for S>
0, whereas the other one vanishescompletely. (b) Larger cross-saturation: one of the two pitch-fork bifurcations on the S =0 line remains for S>
0, the otherone is replaced by a saddle-node (SN) bifurcation.
The phase boundaries in the ( δα, S ) parameter space canbe seen in Fig. 5 for both the case of higher self- andhigher cross-saturation. Starting from the threshold val-ues δα ∗ , of the homogeneous case, we find continuouslychanging phase boundaries ( S ∗ , ) to higher anisotropyvalues for increasing control S . In the case of higher self-saturation there is only one phase boundary S ∗ due tothe absence of a T state. For sufficiently strong control S there is only the stable steady state T . This regionrepresents the desired outcome of a successful switchingprocess. Its phase boundary determines the minimumcontrol strength required to achieve switching for a givenanisotropy value. For lower values of S , a stable four-spot pattern arises and depending on the ratio θ θ β β asecond unstable four-spot pattern occurs together with astability transition of T .Again, we can draw the flow in the ( A , A ) phaseplane and mark the basins of attraction as shown inFig. 6, cases 4, 5, and 6. For θ θ >β β a region with twoattractors occurs similar to the S =0 case, but this timewith a stable F state instead of T , which again implieshysteresis behavior. For increasing anisotropy, the unsta-ble F state will approach T along the unstable manifolduntil they meet and an unstable T and stable F stateemerge. This corresponds to the approach of the twophase boundaries S ∗ and S ∗ and the annihilation of S ∗ , which can only happen for parameter values 1 < θ θ β β < .Otherwise the two boundaries diverge. Thus, the case ofhigher cross saturation is divided into two subcases, de-fined by either the survival or the vanishing of the middleregion with two attractors.Another important effect of the inhomogeneity S is theexplicit breaking of the Z symmetry in the time evolu-tion of A . In the case of larger cross-saturation a cusppoint (C) arises in the two-dimensional bifurcation dia-gram in Fig. 5(b). The pitchfork bifurcation at S =0 isreplaced by a saddle-node bifurcation for S (cid:54) =0, resultingin the creation/destruction of a stable-unstable pair ofF-solutions while varying one of the parameters. In thecusp point two saddle-node phase boundaries, S ∗ and − S ∗ (not shown), meet tangentially. The other pitch-fork bifurcation ( δα ∗ ) remains ( S ∗ ) in presence of theinhomogeneity since the Z symmetry of A is still con-served. This is the only bifurcation in the case of largerself-saturation in Fig. 5(a) for S (cid:54) =0. Figure 7 shows theevolution of the equilibrium surfaces for A and A in the( δα, S ) parameter space and the corresponding bifurca-tions. The pitchfork bifurcations are visible on the S =0line. In the case of higher cross-saturation (lower panel)one of the pitchfork bifurcations unfolds into a saddle-node bifurcation with increasing S . In the case of higherself-saturation (upper panel) the supercritical pitchforkbifurcation remains stable with increasing S .In conclusion, the main difference for the inhomoge-neous case is the absence of a steady T state and theexistence of an additional stable F state in the case ofhigher cross-saturation. This corresponds to the unfold-ing of one of the subcritical pitchfork bifurcations intoa saddle-node bifurcation. The control term S thus pre-vents extinction of A and promotes coexistence respec-tively. Furthermore, the F -T -bistability region van-ishes for high anisotropy values, if the parameter condi-tion 1 < θ θ β β < is satisfied, resulting in a single remain-ing phase boundary S ∗ similar to the case of higher self-saturation. The occurring saddle-node and subcriticalpitchfork bifurcations are problematic since they bothimply sudden vanishing of a stable fixed point, mean-ing the system undergoes an abrupt transition to an-other stable fixed point, i.e. hysteresis is possible. Thiscan only happen in the bistability regime (larger cross-saturation). The numerical simulation of the switchingprocess presented in Section II takes places in the largerself-saturation regime (model parameters calculated inAppendix C), which is advantageous for the switchingpurpose, since no hysteresis can occur. However, if theexternal control beam is not sufficiently strong, we ob-serve a stable F pattern in the numerical simulations,as predicted by the PC model. C. Remarks
So far, we have discussed steady states, their stabil-ity properties, and bifurcations occurring in the solutionspace of the population competition model in Eq. (2).From a dynamical perspective, we observe critical slow-ing down [35, 36] near the phase boundaries due to thecontinuously vanishing real part of the Jacobian’s eigen-value responsible for the bifurcation. This correspondsto the divergence of switching times observed in the nu-merical simulations in Ref. [9]. For example, approach-ing S ∗ , for β β ≷ θ θ from above results in divergenceof the switching time. Similarly, approaching δα ∗ , for β β ≷ θ θ from the left side results in divergence of theback-switching time. Hence, to achieve favorable perfor-mance, switching should be done for parameters suffi-ciently far away from the phase boundaries.We note that in general a nonlinear dynamical sys-tem can have periodic solutions which are characterizedby closed orbits in the phase plane. Here, we use Du-lac’s Criterion [37, p. 202] to rule out any periodic so-lutions. A simplified version reads: The existence of afunction g ( A , A ) with the property that ∇· ( g f ) is signdefinite in the entire considered state space, rules out anyclosed orbits in this area. If we choose g = A A , we obtain ∇ · ( g f )= − β A A + β A A + SA A ), and therefore closedorbits in the positive quadrant are impossible. Anotherobservation is that for θ = θ ≡ θ (symmetric coupling) wecan write the system in Eq. (2) as a gradient field f = ∇ V with the potential function V = (cid:88) i =1 α i A i − β i A i − θ A i A i +1 + S i A i (7)with A i +2 = A i , S =0, and S = S . Closed orbits are im-possible in gradient systems [37, p. 199], however, forthe general case θ (cid:54) = θ (asymmetric coupling) this argu-ment is no longer applicable. Also, assuming the gradi-ent system defined by the potential function (7) allowsus to use the language of catastrophe theory and observethat two of Thom’s seven elementary catastrophes occur,namely folds (denoted as A in Arnold’s notation) andcusps ( A ). We further note that our detailed investiga-tion above is limited to destabilizing linearity, stabiliz-ing nonlinearities, and positive control parameter. Theseconditions match the numerical and experimental obser-vations for the physical system under investigation here.For this case we are able to completely characterize allpossible steady states, their corresponding bifurcations,and rule out any other bifurcations involving double zeroand pure imaginary eigenvalues of the Jacobian. IV. CONCLUSIONS
We have presented a detailed analysis of an all-opticalswitching concept based on transverse patterns in aninteracting polariton system in a planar quantum-wellsemiconductor microcavity. From the relevant equationsbased on a microscopic semiconductor theory here we de-rived a simplified population competition (PC) model de-scribing the system dynamics restricted to selected modes in k -space. Interestingly, the resulting rather simple PCmodel shows all key features of the system dynamicsalso observed in the full numerical simulations in theparameter range of interest here. In addition to whatcan be learned from the numerical simulations, the PCmodel enables us to systematically identify phase bound-aries in parameter space and singularities governing theglobal dynamics of the nonlinear system. Interestinglythe rather complicated original system of interacting mi-crocavity polaritons, for the switching phenomenon stud-ied here can be completely characterized by only sevenflow portraits in a two-dimensional state space. Themodel derived has the very generic form of a generalizedLotka-Volterra (GLV) system extended with an inhomo-geneity term to achieve external control. Such a systemhas not been investigated before. Considering the wide-spread use of GLV systems the understanding obtained inthe present work is of very general nature and will be sim-ilarly applicable also to other fields where external con-trol of population competition is studied such as chemicalreactions or population competition in life-sciences. V. ACKNOWLEDGMENTS
We gratefully acknowledge funding from the DeutscheForschungsgemeinschaft through project SCHU 1980/5-2and through the Heisenberg program (No. 270619725).We further acknowledge valuable discussions with RolfBinder and a grant for computing time at the PaderbornCenter for Parallel Computing (PC ). Appendix A: Numerical Details
We used the following parameters (appropriate forGaAs systems) to numerically simulate the switching pre-sented in section II: m TE =1 . · m TM =0 .
215 meVps µ m − ,Ω=6 . γ c =0 . γ e =0 . α psf =5 . · − µ m ; T ++ = − T + − =5 . · − meV µ m , (cid:15) e = 1 .
497 eV. Coherentcw excitation 2 . I pump =2 . · · I probe = 54 kWcm − . The equations (1) were solvedon a finite 2D grid in real space using a 4th order Runge-Kuttamethod with variable step size. In Section III we use parame-ters β β /θ θ =1 .
14 for the case of larger self-saturation and β β /θ θ =0 .
69 for the case of larger cross-saturation. All bi-furcations were determined analytically for the homogeneouscase and numerically with the help of MATCONT, a matlabsoftware package, for the inhomogeneous case.
Appendix B: Derivation of the PC Model
We follow qualitatively the derivation of the hexagon PCmodel [22], but here we are including polarization effects forlinearly polarized excitation, and therefore considering a dif-ferent reduced k -space. We start with the coupled equationsof motion (1) for the exciton polarization and the cavity fieldand transform them into k -space and in the linear polarization basis. Introducing the reduced k -space consisting of modes { k ≡ , k , k , k , k } with relations k = − k and k = − k (cf. Figure 8) results in 20 equations. We consider an x -FIG. 8. Definition of the reduced k -space.polarized pump. This leads to the following selection rules[19]: • xy -polarized probe with k (cid:107) e x excites TMTE -mode • xy -polarized probe with k ⊥ e x excites TETM -modeWe assume that all dynamical quantities oscillate with pump frequency ω . Removing the phase factor e − i ωt yields ashift of the dispersion by − (cid:126) ω . We consider all phase-matchedscattering processes uo to third order within the reduced k-space. In agreement with the linear stability analysis of thecorresponding system reporting a D ∼ = Z × Z symmetry[20], we assume an equal excitation of opposite modes in thereduced k -space, i.e. p k = p k ≡ p and p k = p k ≡ p (analo-gously for E ). We further assume the pump induced densi-ties E x/y and p x/y at k = 0 to be constant. Using a cyclicdefinition for the two mode pairs p j = p j +2 with j =1 , q = k − k (cid:48) − k (cid:48)(cid:48) ) scatteringprocesses ( q , k (cid:48) , k (cid:48)(cid:48) ) for each mode pair in Eq. (1):( j, , , (0 , j, , (0 , , j ) , × ( j, j, j ) , × ( j, j +1 , j +1) . (B1)In our setup the pump is x -polarized. Therefore, terms ∝ p y or ∝ E y are omitted. Furthermore, we assume the co-linearlypolarized off-axis modes to be very small as for these modesthe instability threshold is not reached, i.e. p xj (cid:28) p yj and E xj (cid:28) E yj . Hence, we also neglect terms ∝ p xj or ∝ E xj . Thisleaves us with the following four equations for the y -polarizedmode pairs:i (cid:126) ∂ t p yj = ( (cid:15) j − (cid:126) ω − i γ e ) p yj + α PSF Ω[ − p y ∗ j p x E x + p x ∗ p yj E x + | p x | E yj + 3 p y ∗ j p yj E yj + 2 p y ∗ j p yj +1 E yj +1 ]+ ( T ++ + T + − )[3 p y ∗ j p yj p yj + 2 p y ∗ j p y j +1 ] − ( T ++ − T + − )[ p y ∗ j p x ]+ T ++ | p x | p yj − Ω E yj (B2)i (cid:126) ∂ t E yj = ( (cid:126) ω yj − (cid:126) ω − i γ c ) E yj − Ω p yj + E y pump ,j . (B3)If we further assume that the time evolution of E followsadiabatically the evolution of p and set ∂ t E j ≈
0, we can write E j = Ω (cid:126) ω j − δ j − (cid:126) ω − i γ c p j ≡ λ j α PSF
Ω e i θ j p j (B4)with j =0 , ,
2. The pump field was also set to zero E pump ,j =0and will be added later manually. We define θ j as the phasebetween E yj and p yj and the ratio of their amplitudes is givenby λ j α PSF Ω . Here, δ j includes anisotropy effects due to higherdensity of states for the TE modes and tilting of the pump. So the parameters are given by λ j e i θ j = α PSF Ω (cid:126) ω yj − δ j − (cid:126) ω − i γ c ) for j = 1 , λ e i θ = α PSF Ω (cid:126) ω x − (cid:126) ω − i γ c ) for j = 0 . (B6)Finally, we obtain two equations for the two elementary statesof the system. We also factorize the exciton field into phaseand magnitude, i.e. p yj =˜ p yj e i ϕ j , and split the results into sep-arate equations of motion for magnitude and phase: ∂ t ˜ p yj = L j ˜ p yj + (cid:88) k =1 C jk ˜ p y k ˜ p yj (B7) ∂ t ϕ j = K j + (cid:88) k =1 D jk ˜ p y k . (B8)If we define φ j ≡ ϕ j − ϕ , the coefficients are then given by (cid:126) L j = − γ e + λ ˜ p x o (sin( θ ) − sin( θ − φ j )) + λ j sin( θ j )(˜ p x − α PSF ) − ( T ++ − T + − )˜ p x sin( − φ j ) (cid:126) C jk = (cid:40) λ j sin( θ j ) , j = k λ k sin( θ k + 2 φ k − φ j ) + ( T ++ + T + − )sin(2 φ k − φ j ) , j (cid:54) = k (cid:126) K j = (cid:15) j − (cid:126) ω − T ++ ˜ p x (cos( θ ) − cos( θ ) − φ j ) − λ j cos( θ j )(˜ p x − α PSF ) − ( T ++ − T + − )˜ p x cos( − φ j ) (cid:126) D jk = (cid:40) − λ j cos( θ j ) − ( T ++ + T + − ) , j = k − λ k cos( θ k + 2 φ k − φ j ) − ( T ++ + T + − )cos(2 φ k − φ j ) , j (cid:54) = k . (B9)Analogously to Ref. [22], we remove the phases as dynamicalvariables by assuming locked phases and linearization. Wedefine the time-dependent phase as φ j ( t ) ≡ δφ j ( t ) + φ (0) j where the locked phases satisfy K j ( φ (0) j )=0 and δφ j ( t ) is a smalldeviation. Expanding equations (B7) and (B8) up to firstorder in δφ j around φ (0) j and neglecting terms ∝ δφ j ˜ p y k leadsto an explicit expression for δφ j which can be substitutedback to obtain ∂ ˜ p yj = (cid:34) L j ( φ (0) j ) + (cid:88) k =1 (cid:18) C jk ( φ (0) j , φ (0) k ) − D jk ( φ (0) j ,φ (0) k ) L (cid:48) j ( φ (0) j ) K (cid:48) j ( φ (0) j ) (cid:19) ˜ p y k (cid:35) ˜ p yj . (B10)We rewrite these two equations in a shorter form: ∂ t ˜ p y = ˜ α ˜ p y − ˜ β ˜ p y − ˜ θ ˜ p y ˜ p y ∂ t ˜ p y = ˜ α ˜ p y − ˜ β ˜ p y − ˜ θ ˜ p y ˜ p y (B11)We replace ˜ p yj with a product of a dimensionless quantity ˆ˜ p yj and a characteristic quantity ˜ p yj,c which carries the originaldimension, i.e. ˜ p yj =ˆ˜ p yj · ˜ p yj,c . We do the same for the indepen-dent time variable t =ˆ t · t c , so that the derivative changes to ∂ ˜ p yj ∂t = ˜ p yj,c t c ∂ ˆ˜ p yj ∂ ˆ t . The characteristic values are chosen in a waythat the corresponding dimensionless quantities are of magni-tude 1. With the definitions ˆ t ≡ t and ˆ˜ p yj ≡ A j , we can finallywrite down the population competition model as ∂ t A = α A − β A − θ A A ∂ t A = α A − β A − θ A A + S (B12)where we have manually added a control parameter S for the A mode pair. Appendix C: Calculation of Model Parameters
The model parameters can be calculated from the physi-cal parameters via equation (B10). Their values, especiallythe signs, depend on the specific choice of the locked phases.The system’s physical behavior observed in the full numeri-cal simulations suggest destabilizing linear terms and stabi-lizing nonlinear terms in the PC model. Chosing phases sat-isfying this condition, characteristic values ˜ p yj,c =1 µ m − and t c =1 ps, and anisotropy effects δ =0 . δ =0 meV,leads to the following model parameters for the switching sim-ulation presented in section II: α =0 . α =0 . β =0 . β =0 . θ =0 . θ =0 . Nature Photonics , 2:185, 2008.[2] Andrew M. C. Dawes, Lucas Illing, Susan M. Clark, andDaniel J. Gauthier. All-optical switching in rubidiumvapor.
Science , 308:672–674, 2005.[3] Dario Ballarini, Milena De Giorgi, Emiliano Cancellieri,Romuald Houdr´e, Elisabeth Giacobino, Roberto Cin-golani, Alberto Bramati, Giuseppe Gigli, and DanieleSanvitto. All-optical polariton transistor.
Nature com- munications , 4:1778, 2013.[4] Alberto Amo, TCH Liew, Claire Adrados, RomualdHoudr´e, Elisabeth Giacobino, AV Kavokin, and A Bra-mati. Exciton–polariton spin switches.
Nature Photonics ,4:361, 2010.[5] Daniele Sanvitto and St´ephane K´ena-Cohen. The roadtowards polaritonic devices.
Nature materials , 15:1061,2016.[6] R. Kheradmand, M. Sahrai, H. Tajalli, G. Tissoni, andL. A. Lugiato. All optical switching in semiconductor mi-croresonators based on pattern selection.
The European Physical Journal D , 47:107–112, 2008.[7] Schumacher Stefan, Kwong N. H., Binder R., andSmirl Arthur L. Low intensity directional switching oflight in semiconductor microcavities. physica status so-lidi (RRL) - Rapid Research Letters , 3:10–12, 2009.[8] Luk, M. H. and Tse, Y. C. and Kwong, N. H. and Leung,P. T. and Lewandowski, P. and Binder, R. and Schu-macher, S. Transverse optical instability patterns in semi-conductor microcavities: Polariton scattering and low-intensity all-optical switching.
Phys. Rev. B , 87:205307,2013.[9] Przemyslaw Lewandowski, Samuel M. H. Luk, ChrisK. P. Chan, P. T. Leung, N. H. Kwong, Rolf Binder, andStefan Schumacher. Directional optical switching andtransistor functionality using optical parametric oscilla-tion in a spinor polariton fluid.
Opt. Express , 25:31056–31063, 2017.[10] M. C. Cross and P. C. Hohenberg. Pattern formationoutside of equilibrium.
Rev. Mod. Phys. , 65:851, 1993.[11] C. Bowman and A.C. Newell. Natural patterns andwavelets.
Reviews of Modern Physics , 70:289 – 302, 1998.[12] H. Meinhardt.
Models of Biological Pattern Formation .Academic Press, London, 1982.[13] J.D. Murray.
Mathematical Biology - II: Spatial Modelsand Biomedical Applications . Springer, New York, 2ndedition, 2003.[14] M.P. Hassell, H.N. Comins, and R. M. May. Spatial struc-ture and chaos in insect population dynamics.
Nature ,353:255 – 258, 1991.[15] Alexey Kavokin, Jeremy J. Baumberg, GuillaumeMalpuech, and Fabrice P. Laussy.
Microcavities . OxfordUniversity Press, Inc., New York, NY, USA, 2008.[16] Vincenzo Ardizzone, Przemyslaw Lewandowski, M. H.Luk, Y. C. Tse, N. H. Kwong, Andreas L¨ucke, Marco Ab-barchi, Emmanuel Baudin, Elisabeth Galopin, Jacque-line Bloch, Aristide Lemaitre, P. T. Leung, PhilippeRoussignol, Rolf Binder, Jerome Tignon, and StefanSchumacher. Formation and control of turing patternsin a coherent quantum fluid.
Scientific Reports , 3:3016,2013.[17] Egorov, O. A. and Werner, A. and Liew, T. C. H. andOstrovskaya, E. A. and Lederer, F. Motion of patternsin polariton quantum fluids with spin-orbit interaction.
Physical Review B , 89:235302, 2014.[18] Saito, Hiroki and Aioi, Tomohiko and Kadokura,Tsuyoshi Order-disorder oscillations in exciton-polaritonsuperfluids.
Physical review letters , 110:026401, 2013.[19] Lewandowski, Przemyslaw and Lafont, Ombline andBaudin, Emmanuel and Chan, Chris K. P. and Leung, P.T. and Luk, Samuel M. H. and Galopin, Elisabeth andLemaˆıtre, Aristide and Bloch, Jacqueline and Tignon,Jerome and Roussignol, Philippe and Kwong, N. H. andBinder, Rolf and Schumacher, Stefan Polarization de-pendence of nonlinear wave mixing of spinor polaritonsin semiconductor microcavities.
Phys. Rev. B , 94:045308,2016.[20] Stefan Schumacher. Spatial anisotropy of polariton am-plification in planar semiconductor microcavities inducedby polarization anisotropy.
Phys. Rev. B , 77:073302, 2008.[21] Schumacher, S. and Kwong, N. H. and Binder, R. In-fluence of exciton-exciton correlations on the polariza-tion characteristics of polariton amplification in semicon-ductor microcavities.
Physical Review B , 76(24):245324,2007.[22] Y C Tse, Chris K P Chan, M H Luk, N H Kwong, P TLeung, R Binder, and Stefan Schumacher. A population-competition model for analyzing transverse optical pat-terns including optical control and structural anisotropy.
New Journal of Physics , 17:083054, 2015.[23] S. Wiggins.
Introduction to Applied Nonlinear Dynami-cal Systems and Chaos . Texts in Applied Mathematics.Springer New York, 2006.[24] L´eon Brenig. Complete factorisation and analytic so-lutions of generalized lotka-volterra equations.
PhysicsLetters A , 133:378–382, 1988.[25] Josef Hofbauer, Karl Sigmund, et al.
The theory of evo-lution and dynamical systems: mathematical aspects ofselection . 1988.[26] Benito Hern´andez-Bermejo and V´ıctor Fair´en. Lotka-volterra representation of general nonlinear systems.
Mathematical biosciences , 140:1–32, 1997.[27] Peter J Wangersky. Lotka-volterra population models.
Annual Review of Ecology and Systematics , 9:189–218,1978.[28] Roger H. Hering. Oscillations in Lotka-Volterra systemsof chemical reactions.
Journal of Mathematical Chem-istry , 5:197–202, 1990.[29] Richard M Goodwin. A growth cycle. In
Essays in eco-nomic dynamics , pages 165–170. Springer, 1982.[30] Willis E Lamb Jr. Theory of an optical maser.
PhysicalReview , 134:A1429, 1964.[31] Joan Roughgarden.
Theory of population genetics andevolutionary ecology: an introduction . 1979.[32] Kie Van Ivanky Saputra, Lennaert van Veen, and GillesReinout Willem Quispel. The saddle-node-transcriticalbifurcation in a population model with constant rate har-vesting.
Discrete & Continuous Dynamical Systems - B ,14:233–250, 2010.[33] Bruno Buchberger.
Ein Algorithmus zum Auffindender Basiselemente des Restklassenringes nach einemnulldimensionalen Polynomideal . PhD thesis, Universityof Innsbruck, 1965.[34] Bruno Buchberger. Ein algorithmisches Kriterium f¨ur dieL¨osbarkeit eines algebraischen Gleichungssystems.
Ae-quationes mathematicae , 4:374–383, 1970.[35] Marten Scheffer, Jordi Bascompte, William A Brock,Victor Brovkin, Stephen R Carpenter, Vasilis Dakos,Hermann Held, Egbert H Van Nes, Max Rietkerk, andGeorge Sugihara. Early-warning signals for critical tran-sitions.
Nature , 461:53, 2009.[36] Jorge R Tredicce, Gian Luca Lippi, Paul Mandel, BasileCharasse, Aude Chevalier, and B Picqu´e. Critical slow-ing down at a bifurcation.
American Journal of Physics ,72:799–809, 2004.[37] Steven H. Strogatz.