A population-competition model for analyzing transverse optical patterns including optical control and structural anisotropy

We present a detailed study of a low-dimensional population-competition (PC) model suitable for analysis of the dynamics of certain modulational instability patterns in extended systems. The model is applied to analyze the transverse optical exciton–polariton patterns in semiconductor quantum well microcavities. It is shown that, despite its simplicity, the PC model describes quite well the competitions among various two-spot and hexagonal patterns when four physical parameters, representing density saturation, hexagon stabilization, anisotropy, and switching beam intensity, are varied. The combined effects of the last three parameters are given detailed considerations here. Although the model is developed in the context of semiconductor polariton patterns, its equations have more general applicability, and the results obtained here may benefit the investigation of other pattern-forming systems. The simplicity of the PC model allows us to organize all steady state solutions in a parameter space ‘phase diagram’. Each region in the phase diagram is characterized by the number and type of solutions. The main numerical task is to compute inter-region boundary surfaces, where some steady states either appear, disappear, or change their stability status. The singularity types of the boundary points, given by Catastrophe theory, are shown to provide a simple geometric overview of the boundary surfaces. With all stable and unstable steady states and the phase boundaries delimited and characterized, we have attained a comprehensive understanding of the structure of the four-parameter phase diagram. We analyze this rich structure in detail and show that it provides a transparent and organized interpretation of competitions among various patterns built on the hexagonal state space.


Introduction
Density patterns spontaneously generated by translational symmetry-breaking instabilities are common occurrences in many extended nonlinear systems [1][2][3][4]. Typically, in the presence of small fluctuations, the spatially uniform state in these systems is destabilized by the systems' internal dynamics in favor of states supporting density patterns such as stripes, squares, hexagons and more complex configurations. Coherent nonlinear optics, driven by laser-matter interactions, is among the areas where these modulational patterns have been fruitfully studied. The optical systems displaying these patterns include lasers, [5][6][7][8] Kerr media and optical parametric oscillators, [9][10][11] atomic gases, [12][13][14][15][16] photorefractive crystals, [17] liquid crystal light valves, [18,19] and, most recently, semiconductor quantum well microcavities [20][21][22][23][24][25]. In most cases, phase-conjugate wave-mixing processes drive directional instabilities in the wave field fed by a uni-directional input light beam. section 5. Section 6 concludes with a summary of results. Appendix A contains a derivation of the PC model equations from the polariton field equations.

The PC model
In this section, we define the PC model and explain the relations between this modelʼs ingredients and microscopic exciton-polariton processes. These relations are demonstrated more concretely in a derivation of the PC model from the microscopic polariton equations in appendix A.
The dynamical variables of our PC model are three real, positive functions of time    ( ) , , 1 2 3 , each of which is the square root of the particle density in a finite transverse momentum mode scaled to the pumped particle density in the spatially uniform mode. We restrict our considerations to only one spin/polarization state. As developed below, the variables satisfy the equations ) saturation, db models anisotropy in the system, γ favors the mutual reinforcement of growth in the three modes, and θ imposes a cross-saturation among the modes. This description of the various terms underlines the fact that equations (1)-(3) as laid out above may also be applicable to the analysis of other pattern-forming systems. The information on the underlying microscopic physics is carried only in the PC model parameters. Figure 1 shows the geometrical relations in transverse momentum space among the polariton modes represented by the model variables. A normally incident pump excites polaritons in a zero (transverse) momentum mode represented by the red circle in the center. Phase conjugate scattering induces modulational instabilities that feed a pair of modes of opposite momenta. Assuming the occupation densities in the pair to be equal, we let one  | | i 2 represent the common density. Thus, a PC model state with only one  i larger than zero represents a 'two-spot' polariton pattern, and one with three non-zero  i ʼs represents a 'hexagon' polariton pattern.
Most terms on the right hand sides of equations (1)-(3) represent polariton scattering processes. We show some of these processes diagrammatically in figure 2. It is noted that the terms with coefficient γ involve scattering of two polaritons from the pump mode and one off-axis (finite transverse momentum) mode ( 3 in figure 2) into two neighboring modes ( 1 and  2 in figure 2). Momentum conservation then implies the three off-axis modes are separated by 60°, which is the motivation for the choice of a hexagonal state space. Since γ favors the hexagon pattern, we will also refer to it as 'hexagon stabilizing'. All the scattering processes represented in equations (1)- (3) are discussed in more detail in [21,22]. We also note that the scatterings are coherent (i.e. phase-dependent) processes. It is shown in appendix A how the phases exert their effects in the Since we are mostly interested in characterizing the qualitative behaviors of the solutions, we can further simplify equations (1)-(3) by taking the following scale transformations: In terms of these scaled variables and parameters, the equations are reduced to 1 2   3  3  3  3   3  1   2  2   2  3  1 2 These equations, with four parameters, b db g S , , , , will be analyzed in the following sections. The three-dimensional space with coordinates    ( ) , , 1 2 3 will be called the state space. We will also use a vector notation with unit vectors    {ˆˆˆ} , , 1 2 3 in this space: . In this notation, we can write equations (10)-(12) as It is straightforward to verify that ´= F 0 . The vanishing of the curl implies that the vector F can be written as In fact, this function can be constructed by inspection:  . As will be seen below, defining this potential function facilitates our organizational understanding of the modelʼs phase diagrams through Catastrophe theory.

Steady states and their stability: construction of phase diagrams
In this section and the next, we characterize the physical steady states, i.e. steady states in the sector , of the dynamical model defined by equations (10)- (12). We seek a global overview of how the qualitative behaviors of these solutions depend on the modelʼs four parameters. To this end, we construct a division of the four-dimensional parameter space into regions, over each of which the number of physical steady states and their stability properties stay unchanged. The steady state solutions vary smoothly with the parameters in each region. When an interface boundary between regions is crossed, the number of solutions may change and/or some stable solutions may become unstable and vice versa. We refer to the 4D parameter space, with this structure embedded, or any lower-dimensional section of it, as a 'phase diagram'. In this section, we give a brief account of the mathematical procedure used to construct these phase diagrams. The resulting phase diagrams are discussed in the next section.
We start by recalling for clarity some well-known mathematical facts. A steady state A 0 of the PC model is a solution to the coupled polynomial equation The stability of A 0 is characterized by the stability matrix ( ) M A 0 , which is given as a matrix function of A by if ¹ ¹ ¹ i j k iand is equal to 0 otherwise. The steady state A 0 is stable if and only if all the eigenvalues of ( ) M A 0 are negative. If one or more eigenvalues are larger than zero, A 0 is unstable. A 0 is called a degenerate steady state if ( ) M A 0 is singular, i.e. one or more eigenvalues of ( ) M A 0 are exactly zero. Otherwise, A 0 is non-degenerate. Noting the definition equation (14) of the potential function ( ) V A , one sees that the steady states are the functionʼs critical points (extrema, saddle points), and the stability matrix is the second derivative, or Hessian, matrix of V : . A critical point of V is also called degenerate when the Hessian matrix at that point is singular.
An application of the implicit function theorem [66] implies that the number of steady state solutions to equations (10)- (12) or the stability of the solutions can only change at those points in parameter space where the stability matrix of at least one steady state is singular, i.e. at least one steady state becomes degenerate. Algebraically, this condition is equivalent to requiring that equations (10)-(12) together with admit at least one solution. Thus, the set of all parameter-space points where these equations have (simultaneous) solutions constitute the phase boundaries of the phase diagrams defined above (or include the phase boundaries as a subset). Plotting these phase boundaries amounts, in principle, to testing for the existence of solutions    ( ) , , 1 2 3 to the set of four polynomial equations 1 2 3 for each point in parameter space. An alternative way to solve this problem is to include one of the parameters, which we choose to be the external control S, as a variable. That is, for each set of b db g ( ) , , values that is scanned, one solves the above four equations for solution(s)    ( ) S , , , . The phase boundaries are then constructed by plotting the solution S as a (possibly multivalued) function of b db g ( ) , , . This way, the number of parameters that need to be scanned is reduced by one. We have found that the second way of computing the phase boundaries is the much faster one, and we have used it for obtaining the results in this paper.
Solving four polynomial equations is in general a nontrivial mathematical task, especially when it needs to be done repeatedly over the entire scanned set of b db g ( ) , , values. This root-finding problem is handled with the method of Gröbner basis [67,68] which essentially provides a solution strategy that algebraically converts the task of solving a multivariate polynomial equation set to that of solving a sequence of single-variable polynomial equations. In general, the converted problem is much easier to solve than the original one. Softwares for both the algebraic reduction and the subsequent solution of the sequence of single-variable equations are available, for example, in Mathematica. We give a brief description of this solution strategy in appendix B.
After computing the phase boundaries, we organize them according to the qualitative nature of their degeneracies. For this purpose, we use the results of Catastrophe theory [66,69,70] which classifies degenerate critical points of smooth functions into generic types according to the behaviors, upon small variations of the parameters, of the functions in the neighborhood of the critical points. The phase boundaries' being identified as consisting of degenerate critical points of the potential function equation (14) makes this application possible. The generic types of critical points of ( ) V A that we encounter in this paper are described in appendix D. Reduction in computational load. We conclude this section with a brief comment on the computational advantage gained from using the algorithm based on Gröbner basis described above.
A 'brute force' way to construct the phase boundary structure in the full 4D parameter space would scan over the four parameters b g db ( ) S , , , and, at each parameter point, would need to search for all stationary solutions through numerical methods like Newtonʼs method. Degeneracy would then need to be tested for each solution so obtained. To perform the same task, our algorithm requires scanning over three parameters, b g db ( ) , , , and, at each point in this 3D subspace, one run of the Gröbner basis program gives directly all values of the last parameter S, at which degeneracy occurs as well as the corresponding degenerate solutions in state space.

Steady state phase diagrams of the PC model
In this section, the steady state phase diagrams of the PC model, obtained in the way described in the previous section, are discussed. We present a 'global' characterization of the steady states, in the sense that we examine the stability of all solutions in each phase region and trace the changes in any of them when a phase boundary is crossed. To set the stage, we first note the action of each term in equations (10)- (12). The terms linear in the dynamical variables tend to drive exponential growth and are responsible for destabilizing the trivial solution ). The single-variable cubic terms, with coefficients β or b db + , are self-saturating and favor solutions with more symmetrical distributions of magnitude among the three variables. The other cubic terms are cross-saturating, i.e. when one variable gains strength, it tends to suppress the other variables through the actions of these terms As a result, these terms favor 'winner-take-all' solutions in which one variable dominate over the other two. The quadratic terms, with coefficients γ, are cross-reinforcing, again favoring more symmetrical solutions. The interplay among these tendencies underlies the competition among spatial patterns in our modelʼs physical interpretation. We begin our description of the phase diagram in the limit g db = = = S 0, i.e. along the β-axis (within the interval b < < 0 1) in the four-dimensional parameter space, where our model equations become sufficiently simple to allow all steady states to be obtained analytically. There are exactly eight steady state solutions in the physical sector (    , , 0 ) of the state space, which we label Going into the general volume of the phase diagram, we use as much as possible the same set of labels on the solutions at each parameter space point, the label assignment being determined by tracing each solution over a continuous path in the phase diagram back to the β-axis. As will be seen, this labeling scheme suffers from ambiguities caused by the non-trivial topology of the phase boundaries. A solution in some regions of parameter space may be traced to two different solutions on the β-axis by traversing two 'non-equivalent' paths. Nevertheless, even with these ambiguities, we have found that the labeling scheme is very useful in generating an organizing overview of the solutions.
Our parameter space (b db g S , , , ) is four-dimensional. To help visualization, we fix the value of β and construct the phase diagram in the 3D space spanned by the other parameters. To reveal more details, we discuss various 2D sections of this phase diagram before displaying the 3D version of the phase boundary surfaces. At the end, we briefly indicate how the 3D phase diagram evolves when β varies between 0 and 1.

Neighborhood of
0, 1: spontaneously broken symmetry limit As mentioned above, we start with characterizing the solutions for parameter points on the β-axis, with g db = = = S 0. Models similar to this limit of our model have been used in many fields like laser physics [71] and biology [72]. Figure 3 shows all the solutions in this case in the physical (    , , 0 ) sector of state space, which can be grouped into 4 types: figure 3, which represents the unstable uniform state of the polariton field.
, are a permutation of 1, 2 and 3. Marked as points figure 3, these solutions are stable for b < < 0 1 and correspond to 2-spot transverse patterns in the far field.
, are a permutation of 1, 2 and 3. Marked as points figure 3, these solutions are unstable and represent 4-spot patterns in the far field. figure 3, this solution is unstable for b < < 0 1 and represents a symmetric hexagonal pattern in the far field.
With db = = S 0, the equations of motion equations (10)- (12) are completely symmetric in the three state variables. The above results say that when self-saturation (β) is weaker than cross-saturation (q = 1), the solution that possesses the symmetry of the model equations (H 4 ) is unstable, while the solutions that break the symmetry (T T T , , 1 2 3 ) are stable. For this reason, we call the present limit the 'spontaneously broken symmetry limit'. For b > 1, T T T , , 1 2 3 are destabilized in favor of the symmetric solution H 4 . The parameter region for b > 1 is, however, not of physical interest for the polariton fluids to which we are applying our model.

2D phase diagram for
We now proceed to describe the structure of the phase boundaries in parameter space and characterize the stationary solutions in each phase region. In this subsection we perform this task for the two-dimensional phase diagram, with coordinates (γ, db), defined by setting b = = S 0.5, 0. As stated in section 3, a phase boundary is a set of parameter space points at each of which one or more stationary solutions become degenerate. Within each phase region demarcated by the boundaries, all the steady states continuously and individually evolve as parameters vary, keeping their stability characteristics.
The simplification brought about by the condition S = 0 enables all steady states and phase boundaries to be found through elementary algebraic steps. Their analytic expressions are listed in appendix C for arbitrary fixed β. Using these, we plot the (g db , ) phase diagram in figure 4 for b = 0.5. In each region, we list the physical steady states grouping them into stable and unstable categories. The naming of the states is determined by tracing the states from the region continuously back to the point g db = ( ) ( ) , 0 ,0, where the states are described in the previous subsection. Throughout the phase diagram in figure 4, the solutions = T i , 1, 2, 3 i remain twospot, i.e., for T i , only  i is nonzero. Where db ¹ 0 the symmetry between  1 and the other two variables is explicitly broken. Among the four hexagon solutions, H 4 and H 1 observe the remaining symmetry   = 2 3 , while H 2 and H 3 do not.
The phase boundaries are labeled I-VI. As mentioned above, their analytical forms are given in appendix C, and the generic types of degenerate critical points that are present-folds, cusps, elliptic umbilics-are described in appendix D. In the following, we describe the qualitative change(s) in the solutions when the phase boundaries are crossed in the direction of increasing γ: Curve I. Each point on this curve is a fold singularity where unstable H 4 and H 1 meet and both disappear on the other side.
Curve II. This is again a curve of folds, but H 4 and H 1 reappear here, H 4 remaining unstable while H 1 changes from unstable to stable.
The arrows indicate the directions of change of these states when γ increases from 0. Right: occupied modes (circles) in transverse momentum space in the instability spatial patterns corresponding to the steady states shown in the left panel. Curve V. Cusps where two independent occurrences take place in parallel. In one, the unstable H 2 , the stable T 2 , and an unstable solution from the non-physical sector meet. Only T 2 survives on the other side but it turns unstable. In the second, H 3 , T 3 , and a non-physical sector solution meet, with T 3 surviving and turning unstable.
Curve VI. Cusps where H 1 and T 1 and an unstable solution from the non-physical sector meet, with only T 1 surviving on the other side but becoming unstable. , where curves V and VI cross, no higher singularity than cusps result because the two curves represent degenerate critical points that are separated in state space.
Another 2D section of the phase diagram which can be constructed 'by hand' is one where we set the asymmetry db to 0. The derivation of the solutions and the phase boundaries go along similar lines as that given in appendix C for the case S = 0 and will not be shown. The resulting (γ , S) phase diagram for b = 0.5 is plotted in figure 5. The state labeling scheme follows that of figure 4 except that the 'H' states are no longer individually designated by subscripts. This is because tracing each of these states continuously back to one state in the spontaneously broken symmetry limit (figure 3) can only be done unambiguously on a 2D surface such as the section shown in figure 4. In the higher-dimensional parameter space, it is possible to start with an 'H' state at g db = = = S 0, trace a continuous closed loop in parameter space, and return to the same parameter point but to a different 'H' state. So a consistent labeling scheme for the 'H' states would be much more complicated than that used in figure 4. Not having pursued this general labeling scheme, we simply record the number of 'H' states in each region. The solutions T i , however, can still be labeled consistently with the scheme in figure 4. The phase boundaries are labeled I-IV, of which boundary IV contains cusp points while the other three are fold lines. The elliptic umbilic point D4 in figure 4 also appears here.

Other 2D sections
In this subsection we discuss two 2D phase diagrams for which none of the parameters are set to zero. The solutions and phase boundaries in these cases cannot be obtained by hand. As explained in section 3, we use the method of Gröbner bases to compute the phase boundaries.  Figure 6 shows a (γ , S) phase diagram with b = 0.5 and asymmetry db = 0.18. Geometrically, then, this 2D section is parallel to the one in figure 5 with db shifted from 0 to the positive side. This shift, reducing selfsaturation, tends to boost  1 at the expense of the other two variables. Consequently, the symmetry between  1 and  3 is broken, and the solutions T 1 and T 3 do not reach degeneracy at the same fold curve, as they do along curve I in figure 5. In effect, curve I splits into two for nonzero db, becoming I′ and I″ in figure 6. T 1 and an H solution are removed as solutions at the higher I″, while T 3 is annihilated at the lower I′. Physically, the advantage given to  1 requires a stronger control beam in direction 2 to destabilize T 1 . The behaviors of solutions across the boundary IV′ is qualitatively the same as those across curve IV in figure 5. Curve II from figure 5 evolves into curve II′ and II″, and curve III into III′ and III″. Basically, as db increases from 0, curves II and III are pushed up in the (γ , S) plane, separated at the elliptic umbilic point, and pulled apart. The apices where II′ and II″ and where III′ and III″ meet are cusp points. It is easily seen from equation (10) that any solution for which  = 0 1 is not affected by changes in db, which explains the fact that the phase boundaries I′ and IV′ occupy the same positions as curves I and IV respectively in figure 5.
Finally, figure 7 shows a 2D phase diagram for b = 0.5 and a fixed g = 0.5. Dominating this diagram is a curved-edge triangle with three fold curves as edges and a cusp point at each vertex. In the 3D picture constructed below, it will be shown that this triangle and the curves II and III (including II′, II″, III′, and III″) are parts of one two-dimensional structure of surfaces.

Phase boundary surfaces in 3D
Although the 2D phase diagrams allow us to view the evolution of the solutions in detail, they may also leave the impression that the phase boundaries form a complicated tangle of lines. For organizational insight, it would be useful to obtain an overall picture of the phase boundaries in a higher dimensional parameter (sub)space. For this purpose, we scan β sparsely and γ and db densely, and for each set of values of b g db , , we solve by Gröbner basis methods, as explained in section 3, for all solutions    ( ) S , , , . The components S of the solutions are extracted, and, when considered as functions of b g db , , , they make up all the boundary 'surfaces' in the 4D parameter space. We put quotation marks around the word 'surfaces' because these point sets are mostly three-dimensional.
To represent these phase boundaries graphically, we fix the value of β at 0.5 and plot the surfaces' sections in (S, db, γ) space in two separate diagrams, figures 8 and 9, for clarityʼs sake. Figure 8 shows a set of surfaces arranged as two triangular pyramids, their apices meeting at a point. The smooth parts of the surfaces contain fold points and the pyramidal edges, where surfaces meet, are cusps. This structure is characteristic of singularity surfaces surrounding an elliptic umbilic singularity [66,69] which is the apex-meeting point (appendix D). Located at g db = ( ) ( ) S , , 0.35, 0, 0 , this singularity point is labeled as  D4 in the 2D diagrams. This elliptic umbilic structure organizes boundary curves in the 2D diagrams: the curves I-IV in figure 4, curves II and III in figure 5, curves II′, II″, III′, III″ in figure 6, and curves I, II, III in figure 7 are different 2D sections of this structure. Figure 9 displays the remainder of the boundary surfaces (for b = 0.5). The red surface contains cusps while the blue surfaces contain folds. The boundary curves in the 2D diagrams which are not associated with the elliptic umbilic structure are sections of these surfaces. Even though it may not be possible to tune the material parameters in the polariton simulation models to access all regions in our 4D parameter space, we have seen that the 3D perspective that we take here generates useful organizational understanding and a simple geometrical picture of the global phase boundary structure.

Interpreting the simulation results of polariton models
In this section, we show several examples of how polariton simulation results [21,22] are qualitatively reproduced in the PC model as dynamical evolutions essentially governed by the structure of the phase diagrams presented in section 4. For this, the dynamical variables    ( ) , , 1 2 3 are traced in time.

Hysteresis at transitions driven by asymmetry in the absence of external control
In simulations with polariton models with hexagonal six-mode state spaces [21], hysteresis behaviors were observed in discontinuous transitions between almost symmetric hexagonal and two-spot patterns (e.g. figure 10 in [21]) when the anisotropy favoring one direction is varied. Through the theory development in section 2, the conditions for these polariton simulations can be qualitatively translated to those for the PC model. The corresponding settings are b < 1, S = 0, γ sufficiently large, and varying db in the phase diagram. We show here that in this parameter region, hysteresis behaviors similar to those in the polariton simulations are indeed obtained. We illustrate this with the particular case of g = 0.88, b = 0.5, S = 0, for which the path of increasing db can be traced in the phase diagram in figure 4. Figure 10(a) shows a zoomedin view of the regions in figure 4 traversed by this path, which crosses two phase boundaries when going in either direction. Below the upper boundary db = » C 0.177 2 , the hexagon state H 1 is stable, whereas above the lower boundary db = » C 0.073 1 , the two-spot state T 1 is stable. Since a bistable region exists between the two curves, db < < C C 1 2 , one expects to see hysteresis behaviors in the adiabatic evolutions of the system when db is changed gradually. Figure 10(b) displays the hysteresis in the solutions to the timedependent model equations under these conditions. On the left of the figure, the system follows the hexagon solution branch H 1 , while on the right it follows the two-spot branch T 1 . In the bistable region, the system stays on the branch it was on when entering the region until this branch disappears or turns unstable at the phase boundary at the far end of the region, where the system jumps to the other solution branch. As stated above, these behaviors reproduce those seen in polariton simulations under the corresponding physical conditions.
Although the system behaves similarly at the two phase boundaries, we recall that the two carry different types of degenerate critical points. The lower boundary db = C 1 , belonging to curve VI in figure 4, is a cusp singularity, where one solution, an unstable T 1 , bifurcates into three: T 1 remains and turns stable, and two unstable solution, H 4 and one in the non-physical sector, are created. As db increases further, H 4 migrates towards H 1 and eventually these two solutions meet and disappear at db = C 2 , which is a fold singularity on curve II in figure 4. It is worth noting that, even though, due to its instability, the influence of H 4 on the systemʼs time evolution is not obvious, this solution is crucial in bringing about the phase boundaries and hence the attendant abrupt pattern changes.

Transitions driven by the control beam
Adiabatic evolutions of the polariton patterns driven by gradual changes in the control beam intensity were also studied in simulations in [21], where continuous evolutions of the patterns punctuated by abrupt jumps also appear. We will demonstrate here that similar behaviors are also obtained in the present model in the parameter region considered in the previous subsection. Figure 11 shows the steady states reached by adiabatic evolution for db = 0.18, g = 0.88 as the control strength S increases from zero. When the control source is off, i.e. S = 0, the anisotropy db favoring  1 produces  3 ) for increasing control strength S, with fixed g = 0.88 and db = 0.18. At S = 0, the anisotropy db favors a two-spot state with only  1 being non-zero. As S increases, the state undergoes two abrupt transitions, first to a hexagon and then to a two-spot state with only  2 being non-zero. As in figure 10(b), each discontinuous jump carries a hysteresis loop, but the parts of the loops for decreasing S are not shown. The path of evolution in the parameter phase diagram is traced by the blue directed line in figure 6. a two-spot state T 1 in which only  1 is nonzero. When sufficiently strong, the control beam biases the system to the two-spot state T 2 . Figure 11 shows that, in its transition under increasing S from T 1 to T 2 , the system passes through two abrupt transitions, the intermediate state being a hexagon. The discontinuous jumps among steady states, like those considered in the previous subsection, can be more systematically understood as the systemʼs crossings of boundaries in the phase diagrams studied in section 4. Figure 6 shows the g ( ) S , phase diagram for an anisotropy (db = 0.18) matching that set in figure 11. The phase trajectory relevant to figure 11 is the vertical i exist, among which only T 1 is stable. As S increases, the states O and T 3 move out of the physical region. In panel (b), the system has just crossed phase boundary III″, where the hexagonal states H 4 (stable) and H 2 (unstable) are created. As S increases further, H 2 and T 1 approach each other (panel (c)) and eventually meet and annihilate each other when the system crosses the fold boundary I″, leaving only H 4 and T 2 (panel (d)). It is at I″ that the system jumps from T 1 to H 4 in figure 11. Then at crossing the cusp IV′, the unstable T 2 bifurcates into three states: a stable T 2 , an unstable H 2 , and a state ¢ T 2 which does not come into the physical region (panel (e)). Then H 4 and H 2 approach and annihilate each other when the boundary III′ is crossed (panel (f)). At III′, the system jumps from H 4 to T 2 in figure 11.
This example shows again the role played by unstable steady states in forming the phase boundaries even though the effects of these states are not obvious when one follows, as in figure 11, only the evolution of the system among the stable patterns. To further illustrate the impact of the quadratic scattering processes, we set g = 0 and trace the steady states through increasing S in figure 13. The corresponding phase diagram trajectory is a vertical line going upwards from the point ( ) 0, 0 in figure 6. A reading of figures 6 and 13 shows, as in previous examples, that qualitative changes in the steady state structure happen at phase boundaries. Briefly, H 2 and H 4 meet and disappear at boundary II′, H 1 and T 3 disappear at I′, and finally T 1 and H 3 disappear at I″, leaving T 2 as the only existing steady state. Figures 12 and 13 can model the switching of the state from T 1 , favored by the anisotropy db, to T 2 by applying the control beam S, favoring  2 . The difference in setting is, we recall, the quadratic scattering strength γ is set to 0.88 for figure 12 and to zero for figure 13. One important resulting difference is that when the control beam is off (S = 0), T 1 is the only stable state in figure 12 while all three two-spot states are stable in figure 13. Physically, if the system is initially prepared in T 1 and a sufficiently strong control beam is applied, the system switches to T 2 in both cases. But if the control beam is subsequently turned off, the system will return to T 1 when quadratic scatterings are acting (figure 12) but will stay in T 2 when these scatterings are absent ( figure 13). Another difference is that for g = 0, since T 1 remains stable as S increases until it disappears when the phase boundary I″ is crossed (between panels (e) and (f) in figure 13), the jump from T 1 to T 2 happens in one step, unlike in figure 11 where the transition goes through an intermediate hexagon state.

Ultrafast switching of patterns by the control beam
Ultrafast switching of patterns have been obtained in theoretical and experimental studies [20][21][22] of quantumwell microcavity polaritons. We show here that, in the approximate parameter regime that we have considered, the solutions of the time-dependent PC model equations (10)- (12) reproduce the ultrafast switching behavior obtained in some polariton model simulations. We have already studied, in the previous subsection, the implications of the steady state phase diagram on adiabatic switching between T 1 and T 2 . We relax the adiabaticity assumption here. Figure 14 shows a switching process by applying a control beam favoring  2 that has a 'square-wave' time dependence. Equations (10)- (12) are solved with the initial state T 1 . At 6 ns, the control beam is switched on, its strength S taking about 0.5 ns to reach its maximum value S max . The system passes through a transient hexagon state, where all three  i are nonzero, before settling into the two-spot state T 2 . The control beam is switched off at 12 ns, causing the system to revert to T 1 . As can be seen, this reversible switching process is repeatable in time. These behaviors are can be well understood by the phase diagram structure studied in figures 6, 11, and 12. For the parameter values in figure 14, with = S S max , the hexagon is not a stable steady state but is close to being one, causing the state trajectory to be drawn to it for 2 to 3 ns before proceeding to its destination.
The features of the PC model solutions described in this section were observed in polariton simulations [20,21]. These features are analyzed here as parts of the parameter phase diagram of steady states. It would be interesting to investigate whether other parts of the phase diagram are accessed by realistic polariton simulations.

Summary and outlook
We have developed a three-variable PC model to gain physical insight into the dynamics of competition among various extended instability-driven patterns in semiconductor quantum well microcavities. The patterns covered by the PC model are limited to those constructed by six or fewer single-spin, momentum space modes, but their competitions as four physical parameters are varied present a rather rich structure. To better understand the competition dynamics, we have constructed a phase diagram of steady state solutions in parameter space. Each phase region is characterized by the number and stability types of steady states. The boundaries separating the regions, where these solution attributes change, govern the qualitative changes of the systemʼs evolution. The phase diagram has helped in interpreting the pattern competitions seen in more detailed polariton simulations.
The PC model can be extended to analyze polariton systems with two spin states being populated, which carry pattern dynamics that differ in some important ways from those in single-spin systems [73][74][75][76].  Hartree-Fock approximation. E k,inc is the external input light field at momentum k, and t c is its coupling to the cavity field. The state space of the equations is first truncated to six modes lying on the vertices of a regular hexagon centered at the origin of the | | k -space, whose radius is determined by the simulation results. After removing the fast oscillatory factor w e t i p , the reduced field equations read as  , the dimensionless coefficients appearing in (A.8) and (A.9) can be expressed as follows: with ( ) X n indicating the nth order derivative of X evaluated at the initial locked phase, and, for reference, A.20 We also define For a completely symmetric (isotropic) system, it is straightforward to see that a a = i , g g = i , b b = i and q q = i are independent of i and hence the dynamics of the system is governed by 0, , , 0 ; (iv) solve each equation for x 3 and retain those solutions ( ) x x x , , 1 2 3 that solve all equations. At each step, the task is to solve for all solutions to a single-variable polynomial equation, which can be achieved with standard numerical, or where possible, analytic means. As stated above, the set of roots obtained in this way for the Gröbner basis is by construction the same as that for the original starting polynomial set.
We note that the Gröbner basis is not unique. It depends, for example, on the way the unknowns (( ) in the above example) are chosen to be ordered in generating the basis. The Gröbner basis method is used in this paper to compute the boundaries in the PC model phase diagram. For this problem, we use available programs in Mathematica for both the construction of the Gröbner basis polynomials and the back substitution solution for their roots.
i are the right hand sides of equations (10)-(12) respectively. Multiplying equation (C.2) by  2 and equation (C.3) by  3 and subtracting one equation from the other, one finds after some algebra implying that the hexagon solutions must satisfy either the condition   = ( ) C.5 2 3 or C . 6 With each solution of this equation, one uses equation (C.2) again to obtain We label this group of solutions H 1 and H 4 by continuity from the solutions on the β-axis (with g db = = = S 0). The last group of solutions are obtained by using equation (C.6) with equations (C.1) and (C.2) to eliminate  2 and  3 , resulting in a linear equation for  1 2 , yielding the answer With  1 known,  2 ,  3 are given by the solutions of the quadratic equation which is plotted as curve V. For T 2 or T 3 , the condition yields which is plotted as curve VI. Both curves are made up of cusp points.
, has zero eigenvalue(s). In simple cases, elementary Catastrophe theory [66,69,70] classifies degenerate critical points into forms according to how the critical point structure changes when small perturbations are added to the function (V in our case) around the degenerate critical point. When applied to the function V, we find that the phase boundaries are parameter values where the degenerate critical points fall into three types: folds, cusps, and elliptic umbilics. We briefly summarize in this appendix some characteristics of these singularities relevant to our discussions.
A fold point in parameter space is where the stability matrix of one critical point has one zero eigenvalue, and a shift of the state space origin to the critical point and a smooth variable transformation