Nonlinear gap modes and compactons in a lattice model for spin-orbit coupled exciton-polaritons in zigzag chains

We consider a system of generalized coupled Discrete Nonlinear Schr\"{o}dinger (DNLS) equations, derived as a tight-binding model from the Gross-Pitaevskii-type equations describing a zigzag chain of weakly coupled condensates of exciton-polaritons with spin-orbit (TE-TM) coupling. We focus on the simplest case when the angles for the links in the zigzag chain are $\pm \pi/4$ with respect to the chain axis, and the basis (Wannier) functions are cylindrically symmetric (zero orbital angular momenta). We analyze the properties of the fundamental nonlinear localized solutions, with particular interest in the discrete gap solitons appearing due to the simultaneous presence of spin-orbit coupling and zigzag geometry, opening a gap in the linear dispersion relation. In particular, their linear stability is analyzed. We also find that the linear dispersion relation becomes exactly flat at particular parameter values, and obtain corresponding compact solutions localized on two neighboring sites, with spin-up and spin-down parts $\pi/2$ out of phase at each site. The continuation of these compact modes into exponentially decaying gap modes for generic parameter values is studied numerically, and regions of stability are found to exist in the lower or upper half of the gap, depending on the type of gap modes.


I. INTRODUCTION
Planar semiconductor microcavities operating in the exciton-polariton regime have become a paradigm model for experimental and theoretical studies of nonlinear and quantum properties of light-matter interaction [1]. A major advantage of these systems is that they are solid state devices, that operate in a wide diapason of temperatures between few Kelvins and up to the room conditions. Interaction between the polaritons is much stronger than for pure photons, that lowers power requirements for creating conditions when polariton dynamics can be effectively controlled with the external light sources [2]. Microcavities can also be readily structured to create a variety of potential energy landscapes reproducing lattice structures known in studies of electrons in condensed matter on more practical scales of tens of microns. Thus polaritons can be controlled using band gap and zone engineering [3]. Through their peculiar spin properties and sensitivity to the applied magnetic field, polaritons in structured microcavities have been shown to have a number of topological properties [4]. Thus polariton based devices have a competitive edge over their photononly counterparts through their relatively low nonlinear thresholds and possibility to create micron-scale topological devices. A combination of these two aspects has been recently used to demonstrate a variety of nonlinear topological effects in polariton systems, see, e.g., [5] and references therein.
As a specific example, a polariton BEC in a zigzag chain of polariton micropillars with photonic spin-orbit coupling, originating in the splitting of optical cavity modes with TE and TM polarization, was proposed in Ref. [6]. The simultaneous presence of zigzag geometry and polarization dependent tunneling was shown to yield topologically protected edge states, and in the presence of homogeneous pumping and nonlinear interactions the creation of polarization domain walls through the Kibble-Zurek mechanism, analogous to the Su-Schrieffer-Heeger solitons in polymers, was numerically observed [6]. Of crucial importance is the spin-orbit induced opening of a central gap in the linear dispersion relation. As we will show in this work, the existence of a gap, together with the option of tuning the linear dispersion towards flatness at specific parameter values, also leads to nonlinear strongly localized modes in the bulk (intrinsically localized modes) with properties depending crucially on the relative strength of interaction between polaritons of opposite and equal spin.
The starting point is the following set of two coupled continuous Gross-Pitaevskii equations [7]: These equations describe exciton-polaritons with circularly polarized light-component, where Ψ + corresponds to left (positive spin) and Ψ − to right (negative spin) polarization. Polaritons interact mainly through their excitonic part, and interactions between polaritons with identical polarization are generally repulsive (here normalized to +1), while interactions between those of opposite spins often are weaker and attractive. A typical value is a −0.05 [8], but may range between roughly −1 a 0, and may possibly be also repulsive, or attractive with a magnitude stronger than the self-interaction [9]. Since the exciton-components of the polariton wave functions typically are localized within small spatial regions, the interactions are assumed to be local (point interactions) in this mean-field description. Ω describes the Zeeman-splitting between spin-up and spin-down polaritons in presence of an external magnetic field; in this work we put Ω = 0.
Of main interest here is the term proportional to β: it arises due to different properties associated with polaritons whose photonic components, as expressed in a suitable basis of linear polarization, have TE resp TM polarizations (or, alternatively, longitudinal/transversal w.r.t. the propagation direction (k-vector)). It is commonly described in terms of different effective masses of the lower polariton branches for TE and TM components, β ∝ m −1 T E −m −1 T M , whose ratio typically may be of the order m T E /m T M ≈ 0.85 − 0.95 (see e.g. supplemental material of [10]), although in principle β could have arbitrary sign. Expressed in a basis of circular polarization (spinor basis) as in (1) (Ψ ± = Ψ x ∓ iΨ y ), this TE/TM energy splitting can be interpreted as a spin-orbit splitting, since the dynamics of the two spin (polarization) components couple in a different way to the orbital part of the other component (via derivatives in x and y of the mean-field wave function in (1)).
In this work, we choose the potential V (x, y) as a zigzag potential along the x-direction, considering this geometry as the simplest generalization of a straight 1D chain which yields non-trivial geometrical effects of the spin-orbit coupling between polaritons localized at neighboring potential minima. As an example potential, we may choose e.g.: as illustrated in Fig. 1. Here d is the distance between potential mininma, and 2N is the total number of potential wells in the chain. The geometry is essentially the same as for the coupled micropillars in [6], with all angles for the links between neighboring minima being ±45 • with respect to the x-axis. Evidently one may easily generalize to arbitrary angles, or more complicated expressions for zigzag potentials which may be realized in various experimental settings e.g. with optical lattices [11]. In order to motivate a tight-binding approximation, we assume V 0 1. In order to understand the most important effects of the spin-coupling coupling in (1) in a tight-binding framework, we here consider situations where the effects of spin-orbit splitting inside each potential well can be neglected, and only are relevant in the regions of wavefunction overlap between neighboring wells. For the experimental set-up of [10], this should be a good approximation if the spatial modes inside the wells may be approximated with Laguerre-Gauss modes with zero orbital angular momentum (LG ± 00 in the notation of [10], where the two subscripts stand for radial and orbital quantum numbers of the 2D harmonic-oscillator wave function, and the superscript indicates polarization as in (1).) At least for a single cavity of non-interacting polaritons, these modes should be good approximations to the ground state, so let us assume that interactions (nonlinearity) and spin-orbit couplings are sufficiently weak to be treated perturbatively, along with the inter-well overlaps. The approach may be extended to consider also lattices of spin vortices (excited modes) built up from modes with nonzero OAM (e.g. LG ± 0±1 as considered in [10]); however this will introduce some additional complications and will be left for future work.
Moreover, if V 0 1 we may also neglect the effect of next-nearest-neighbor interactions (distances between two wells in the horizontal x-direction is √ 2 times larger than between nearest neighbors). It may then be a good approximation to use, as the basis set for the tight-binding approximation, the Wannier functions for a full 2D square lattice (these issues are discussed and numerically checked for some realization of a zigzag optical lattice in a recent Master thesis [12]). These may resemble (but certainly differ from) [12] the LG individual modes (e.g. Wannier functions typically have radial oscillatory tails, decaying exponentially rather than Gaussian). In any case, we will assume that the basis functions w(x, y) (expressed in Cartesian coordinates) are qualitatively close to the LG 00 modes. Particularly, they will be assumed to be close to cylindrically symmetric (w(x, y) ∼ e −ω(x 2 +y 2 ) in the harmonic approximation). (Note that this assumption would not be valid for spin vortices arising from LG modes with nonzero OAM.) The outline of this paper is as follows. In Sec. II we derive the tight-binding model, discuss its general properties, and illustrate the linear dispersion relation for the case with ±45 • angles which will be the system studied for the rest of this paper. We also in Sec. II E identify a limit where the linear dispersion relation becomes exactly flat, and identify the corresponding fundamental compact solutions. In Sec. III we construct the fundamental nonlinear localized modes in the semi-infinite gaps above or below the linear spectrum, as well as in the mini-gap between the linear dispersion branches, opened up due to the simultaneous presence of spin-orbit coupling and nontrivial geometry. Analytical calculations using perturbation theory from the weak-coupling and flat-band limits for the semi-infinite and mini-gap, respectively, are compared with numerical calculations using a standard Newton scheme. In Sec. IV the linear stability of the different families of nonlinear localized modes is investigated, and some instability scenarios are illustrated with direct dynamical simulations. Finally, some concluding remarks are given in Sec. V.

A. Derivation of the tight-binding model
Under the above assumptions, we may expand: where, relative to the coordinate system of (2), d = d/ √ 2, x = x − d /2, y = y − d . Note that the (Wannier) basis functions are the same for both components, since we have assumed no spin-orbit splitting inside the wells, Ω = 0, and w are basis functions of the linear problem. Note also that an analogous approach was used in [13] to derive lattice equations for the simpler problem of a pure 1D lattice with a standard spin-orbit coupling term (−i∂ x , linear in the spatial derivative) for atomic BEC's in optical lattices; similar models were also studied in [14][15][16]. For simplicity we will assume below that w(x, y) can be chosen real (which is typically the case in absence of OAM; the generalization to modes with nonzero OAM requires complex w(x, y) and will be treated in a separate work).
Inserting the expansion (3) into (1), we obtain for the first component: +V (x, y) n u n w(x − n d , y − (−1) n d /2). (4) Here, in writing the nonlinear term as a simple sum and not a triple, we have neglected overlap between basis functions on different sites in cubic terms in w (assuming strong localization of w).
Multiplying with w (n) ≡ w(x − nd , y − (−1) n d /2), integrating over x and y, using the orthogonality of Wannier functions and neglecting all overlaps beyond nearest neighbors, we obtain from (4) a 1D lattice equation of the following form for the site amplitudes of the spin-up component: Here the coefficients are: On-site energy, linear coupling coefficients, where the second equality is obviously true if w (n) is cylindrically symmetric; nonlinearity coefficient, on-site spin-orbit interaction, which is identically zero if w (n) is cylindrically symmetric (easiest seen in polar coordinates, with w = w(r) only, ω = β 2π 0 dφe −2iφ rdr(w rr − wr r )w = 0) (but generally nonzero if Wannier modes would have OAM); and nearestneighbor spin-orbit interactions (the relevant 'new' terms here), Since tails of w are exponentially small, we may assume all integrals taken over the infinite plane. Explicitly, with a change of origin we may write e.g. the first term in the integral in (6) as w xx (x ∓ d , y − (−1) n d )w(x, y)dxdy, etc. But for the case with w cylindrically symmetric, we may easier evaluate the integral (6) in polar coordinates, centered at site n ± 1. After some elementary trigonometry we then obtain: Letting φ = π 4 ± (−1) n φ, this can be expressed as where α n ≡ (−1) n π 4 are the angles for the links in the zigzag chain with respect to the x-axis, and the integral defining σ is independent of all signs since cos φ is even. Explicitly, for the π/4 zigzag chain we get Proceeding analogously with the second component, we obtain the corresponding lattice equation for the siteamplitudes of the spin-down component: Here, , Γ, γ are identical as for the first component (i.e., we may put = 0 by redefining zero-energy, and γ = 1 (or alternatively Γ = 1) by redefining energy scale). For the on-site spin-orbit interaction, (note opposite sign of third term compared to ω), which is again zero if w is cylindrically symmetric. And finally, for the nearest-neighbor spin-orbit couplings, (again note sign of third term compared to (6). As before, restricting to cylindrically symmetric w yields where the last equality holds since the integral is equivalent to that of (8). Explicitly, for the π/4 zigzag chain i.e., with opposite signs compared to (9). Note that, under the above assumptions (w real and cylindrically symmetric), the integral defining σ is always real. We note that the resulting lattice equations (5), (10), with spin-orbit coefficients given by (9), (13), are not equivalent to the equations studied in [13][14][15][16]. In particular, we comment on the relation between the present model and that of Ref. [16], who considered a diamond chain with angles π/4 and a Rashba-type spin-orbit coupling. The zigzag chain could be considered as e.g. the upper part of the diamond chain, if all amplitudes of the lower strand would vanish. However, because the spin-orbit coupling used in [16] is linear in the spatial derivatives while in this work it is quadratic, the spin-orbit coupling coefficients in [16] have a phase shift of π/2 between diagonal and antidiagonal links, compared to π in (9), (13).

B. General properties of the TB-equations
Let us put = 0 and γ = 1. As above, assuming cylindrically symmetric basis functions, we have ω = ω = 0. With no external magnetic field, Ω = 0. Equations (5) and (10), with spin-orbit coefficients given by (9) and (13), respectively, then become: One may easily show the existence of the "standard" two conserved quantities for DNLS-type models; Norm (Power): and Hamiltonian: Here, {u n , v n } and {iu * n , iv * n } play the role of conjugated coordinates and momenta, respectively (i.e.,u n = ∂H/∂(iu * n ),v n = ∂H/∂(iv * n ), etc.). We may note that the Hamiltonian is similar to the Hamiltonian for the "inter-SOC" chain of Beličev et al. (Eq. (11) in [15]), but differs by the "zigzag" spin-orbit factor (−1) n in the last term. Note that this factor can be removed by performing a "staggering transformation" on the site-amplitudes of the spin-down component: v n = (−1) n v n , transforming the equations of motion (14) into: Thus, this transformation effectively changes the sign of the linear coupling of the spin-down component into Γ → Γ = −Γ (which may be interpreted as a reversal of the "effective mass" of the spin-down polariton in this tight-binding approximation), while the nonlinear and spin-orbit terms for both components become equivalent. Eqs. (17) differ from the equations derived in [13] for the straight chain with standard spin-orbit coupling only through this signchange of Γ for the spin-down component. Note also that Eqs. (17) are invariant under a transformation v n → −v n , n + 1 → n − 1, i.e., an overall change of the relative sign of the spin-up and spin-down components is equivalent to a spatial inversion.

C. Generalization to arbitrary angles
As mentioned, it is straightforward to generalize the derivation of the tight-binding equations to arbitrary bonding angles α = π/4 in the zigzag chain. We just outline the main steps: In (3) and the following, we replace y − (−1) n d /2 with y − (−1) n tan(α)d /2 (having redefined d = d cos α). In (7), (12), π/4 then get replaced by α, as already indicated. In (9) we get e −i2α σ for diagonal links and e i2α σ for antidiagonal, and in (13) we get e i2α σ for diagonal links and e −i2α σ for antidiagonal. Then in the tight-binding equations of motion (14), the last term for the spin-up component gets replaced by +e (−1) n i2α v n+1 − e −(−1) n i2α v n−1 , and for the spin-down component by +e −(−1) n i2α u n+1 − e (−1) n i2α u n−1 . For the rest of this paper we will assume α = π/4 and leave the study of effects of variation of the binding angle to future work.

D. Linear dispersion relation
Let u n = ue i(kn−µt) , v n = ve i[(k+π)n−µt] (i.e., v n = ve i(kn−µt) removing factors (−1) n ), with |u|, |v| 1. Inserting it into (14) (or (17)) and neglecting the nonlinear terms then yields: Thus, as illustrated in Fig. 2 (assuming σ < Γ), the spin-orbit coupling opens up gaps in the linear dispersion relation at k = ±π/2, of width 4σ. Note that in contrast to the models for straight chains studied in [13,15], no external magnetic field is needed to open the gap for the zigzag chain. The gap opening is a consequence of the simultaneous presence of spin-orbit coupling and nontrivial geometry, which was also noted for the more complicated diamond chain in [16].
The amplitude ratios between the components may be obtained as v/u = −Γ cos k∓ √ Γ 2 cos 2 k+σ 2 sin 2 k σ sin k . For weak spin-orbit coupling (σ Γ), the polariton is mainly spin-up (u v) on the lower dispersion branch and spin-down on the upper branch (v u) when Γ cos k > 0, and the opposite when Γ cos k < 0.

E. Flat band and compact modes
We may also note that in the particular case of |Γ| = |σ|, the dispersion relation becomes exactly flat. In this case, there are eigenmodes completely localized on either upper or lower part of the chain, with alternating v n = ±iu n on this part (i.e., v n ≡ u n ≡ 0 either for odd or even n). These modes persist also in the presence of nonlinearity (interactions). With the flat band, it is also possible to construct exact compact solutions localized on two neighboring sites. Explicitly, we get for σ = +Γ: and for σ = −Γ: In both cases, the nonlinear dispersion relation for these compactons yields µ = (1 + a)|A| 2 ∓ 2Γ. We will discuss further properties of these nonlinear compactons (e.g. stability) below.

III. NONLINEAR LOCALIZED MODES
A. Single-site modes above the spectrum in the weak-coupling limit For the case of no spin-orbit coupling (σ = 0 and small Γ 1), analysis of fundamental nonlinear localized solutions (including their linear stability) of (14) was done in [17]. It would be straightforward to redo a similar extensive analysis including also a small σ, but it is not the main aim of this work. We focus here first on discussing the effect of small coupling on polaritons with main localization on a single site n 0 .
In the limit of Γ = σ = 0 ("anticontinuous limit"), stationary solutions of (17) are well known. There are two spin-polarized solutions: (spin-down); and one spin- e iθ , with an arbitrary relative phase θ between the spin components. Comparing the Hamiltonian (16) for these solutions at given norm P , we have H = P 2 /2 for the spin-polarized modes and H = (1 + a)P 2 /4 for the mixed mode, so the mixed mode has lowest energy as long as a < 1.
When µ > 0 does not belong to the linear spectrum (18), we search for continuation of these modes for small but nonzero Γ, σ into nonlinear localized modes with exponentially decaying tails and frequency above the spectrum. (We here assume a > −1; if a < −1 the localized modes arising from the spin-mixed solution will have µ < 0 and thus lie below the spectrum.) We may calculate them explicitly perturbatively to arbitrary order in the two small parameters Γ, σ; here we give only the first-and second-order corrections to the five central sites (amplitudes of other sites will be of higher order): It can be seen from such expressions (extending to higher orders) that amplitudes do decay exponentially above the spectrum, µ > 2Γ. However, for spin-mixed modes with |u n | = |v n |, it is important to remark that even though the second-order corrections in (23) can be obtained for arbitrary relative phases θ, the fourth-order correction to site n 0 can be made consistent with the condition |u n | = |v n | only if Γ 2 σ 2 sin(2θ) = 0. Thus, since a solution with θ = π is equivalent to θ = 0 through spatial reflection in the central site, the only non-equivalent single-site centered spin-mixed modes existing for nonzero Γ and σ have θ = 0, π/2. We also remark that, for a stationary and localized solution, current conservation imposes the general condition: Numerically calculated examples for the spin-up and spin-mixed modes are illustrated in Fig. 3. Note from (23) that, for the spin-mixed mode with θ = π/2, |u n0+1 | 2 + |v n0+1 | 2 = |u n0−1 | 2 + |v n0−1 | 2 , i.e., the reflection symmetry around the central site gets broken on the opposite sublattice (upper or lower part of the chain) if there is a nontrivial phase-shift between the spin-up and spin-down components at the central site. As |Γσ|/µ 2 increases the spatial asymmetry increases (Fig. 3 (d)), until the solution typically bifurcates with an inter-site centered (two-site) mode with equal amplitudes at sites n 0 and n 0 + 1 before reaching the upper band edge at µ = 2Γ.

B. Fundamental gap modes from the flat-band limit
Since the gap in the linear spectrum opened by the spin-orbit coupling at k = ±π/2 appears only when Γ and σ are both nonzero, the standard anticontinuous limit Γ = σ = 0 is not suitable for constructing nonlinear localized modes with frequency inside this gap ("discrete gap solitons"). Instead, we may use the flat-band limit |Γ| = |σ| = 0, where the exact nonlinear compacton modes (19)- (20) can be used as "building blocks" for the continuation procedure. Analogously to above, we may then calculate gap solitons perturbatively in the small parameter |Γ| − |σ|. To be specific, we assume a > −1, Γ ≥ σ > 0, and consider the continuation of a single two-site compacton from the lower flat band µ = −2Γ into the gap. From the limiting solution (19) with the upper sign, we then obtain the lowest-order corrections to six central sites (amplitudes at other sites are of higher order) as: This family of fundamental gap modes (called type I gap modes) can be continued throughout the gap, with a numerical example illustrated in Fig. 4 (a). Profiles of another two types of gap modes numerically found to exist as nonlinear continuation of fundamental compactons, are depicted in Fig. 4 (b,c). Family of gap modes of type II ( Fig.  4 (b)) originates from compact solution which is superposition of two neighboring overlapping in-phase compactons.
On the other hand, type III gap modes evolve in the presence of nonlinearity from superposition of two neighboring overlapping compactons with a π/2 phase difference (Fig. 4 (c)).

IV. LINEAR STABILITY OF NONLINEAR LOCALIZED MODES
Linear stability of the above modes can be checked from the standard eigenvalue problem. If we denote the amplitudes of the exact stationary modes of (17) as {u n }, we may express the perturbed modes as u n = u Inserting into (17) and linearizing, we obtain the following linear system of equations for the perturbation amplitudes {c n , d n , f n , g n }: Linear stability is then equivalent to (26) having no complex eigenvalues. We may easily solve it for the uncoupled modes. Due to the overall gauge invariance of (17) (u n → e iφ u n , v n → e iφ v n ), there are always two eigenvalues at λ = 0. For the spin-polarized modes, the remaining two eigenvalues are at λ = ±(1 − a)µ, while for the spin-mixed mode there is a fourfold degeneracy at λ = 0. The latter is explained by the arbitrary phase difference θ between the u and v components for this mode. To see whether linear stability of the fundamental modes survives switching on the couplings Γ, σ, we first note that the linear spectrum of (26) corresponding to sites with u (0) n ≡ v (0) n ≡ 0 has four branches, at λ ∈ ±[µ − 2Γ, µ − 2σ] and λ ∈ ±[µ + 2σ, µ + 2Γ]. Thus, unless a = 0, 1, or 2, we see immediately that the fundamental spin-polarized modes must remain linearly stable at least for small couplings. The general stability properties for larger Γ and/or σ will be discussed below for the different fundamental modes separately.
A. Spin-polarized modes above the spectrum Typical results from numerical diagonalization of (26) for the family of fundamental spin-polarized modes above the spectrum are shown in Fig. 5. As is seen, these modes are linearly stable in their full regime of existence when a < 1. The magnitude of the frequency of the internal eigenmode arising from local oscillations at the central site lies above the linear spectrum when a < 0 (Fig. 5 (a)) and below the linear spectrum when 0 < a < 1 (Fig. 5 (b)). In both cases, it smoothly joins the band edge as µ → 2Γ (linear limit), without causing any resonances. On the other hand, for a > 1, the Krein signature of this eigenmode will change, as a consequence of the spin-polarized mode now having a lower energy than a spin-mixed mode, and thus it is no longer an energy maximizer for the system. This results in small regimes of weak oscillatory instabilities when the internal mode collides with the linear spectrum for frequencies close to the band edge, as shown in Fig. 5 (c). Here, purely imaginary eigenvalues are represented by green (light gray) triangles, and complex eigenvalues are represented by blue (dark gray) squares and red (middle gray) circles for their real and imaginary parts, respectively.

B. Spin-mixed modes above or below the spectrum
For the fundamental spin-mixed modes continued from (23), the four-fold degeneracy of zero eigenvalues resulting from the relative phase θ is generally broken for non-zero coupling as only modes with integer 2θ/π can be continued, and moreover the structures of modes with θ = 0 and θ = π/2 become non-equivalent. We discuss here first the case θ = 0, and show in Fig. 6 typical results from numerical diagonalization for different values of a. First, for a < −1, as remarked above the spin-mixed modes lie below the linear spectrum (µ < −2Γ), and the pair of eigenvalues originating from λ = 0 in the anticontinuous limit (µ → −∞) generally goes out along the imaginary axis ( Fig. 6 (b)), where it remains. Thus, spin-mixed modes with θ = 0 and a < −1 are generically unstable. On the other hand, when a > −1 the spin-mixed modes lie above the linear spectrum (µ > 2Γ), and for −1 < a < 1 this eigenvalue pair goes out along the real axis ( Fig. 6 (a)). Thus, these modes remain linearly stable for sufficiently large µ (or, equivalently, weak coupling), but become unstable through oscillatory instabilities (complex eigenvalues, see Figs. 6 (c,d)) as they approach the linear band edge with widening tails, causing resonances between the local oscillation mode at the central site and modes arising from oscillations at small-amplitude sites.
An example of the dynamics that may result from the oscillatory instabilities of the spin-mixed modes in this regime is shown in Fig. 7. Note that, after the initial oscillatory dynamics, the solution settles down at the stable fundamental spin-up mode (in this particular case the mode center is also shifted one site to the right).
As illustrated in Figs. 6 (c,d), the stability regime increases for a increasing towards 1, and exactly at a = 1 the spin-mixed states are always stable. However, for a > 1 the eigenvalue pair originating from zero again goes out along the imaginary axis (not shown in Fig. 6) as for a < −1, and thus spin-mixed modes with θ = 0 are generally unstable also for a > 1. In fact, this latter instability can be considered as a stability exchange with the θ = π/2 spin-mixed mode, which, as illustrated in Fig. 8, is generally unstable with purely imaginary eigenvalues for a < 1 (Figs. 8(a,b)) but stable for a > 1 (Fig. 8(c)).

C. Compact modes in the flat-band limit
In the flat-band limit, we may obtain exact analytical expressions for the stability eigenvalues of the single two-site compacton modes. We focus as above on the specific case with Γ = σ > 0 and a > −1, when the nonlinear compacton originating from µ = −2Γ (Eq. (19) with upper sign) enters the mini-gap for increasing µ. For all zero-amplitude sites, the eigenvalues are just those corresponding to the flat-band linear spectrum, λ = ±µ ± 2Γ. For the compacton sites, four eigenvalues correspond to local oscillations obtained by eliminating the surrounding lattice: λ = 0 (doubly degenerate as always) and λ = ±2 2Γ(µ + 4Γ). Since the eigenvalues of these internal modes are always real for Γ > 0 and they do not couple to the rest of the lattice, they do not generate any instability. The remaining eigenvalues describe the modes coupling the perturbed compacton to the surrounding lattice, and are obtained from the subspace with The rather cumbersome result can be expressed as: Oscillatory instabilities are generated if the expression inside the square-root in (27) becomes negative. Since obtaining explicit general expressions for instability intervals in µ would require solving a nontrivial fourth-order equation, we show in Fig. 9 numerical results for the specific parameter values a = ±0.5 and 1.5. As can be seen, the compacton remains stable throughout the mini-gap as long as a ≤ 1 but develops an interval of oscillatory instability in the semi-infinite gap above the spectrum. The instability interval vanishes exactly at a = 1, but then moves into the upper part of the mini-gap for a > 1. Purely imaginary eigenvalues, resulting from the terms outside the square-root in (27) becoming negative, also appear in the semi-infinite gap for a > 1.

D. Gap modes in the mini-gap
For the fundamental (type I) gap mode continued from the single two-site compacton (25) (assuming again a > −1 and Γ > σ > 0), we illustrate in Fig. 10(a) typical results for the numerical stability analysis. As can be seen, as σ decreases from the compacton limit σ = Γ, weak instabilities start to develop mainly close to the two gap edges. A further decrease in σ yields instabilites in most of the upper half of the gap, while the mode remains stable in large parts of the lower half. Comparison with the stability eigenvalues for the exact compacton ( Fig. 9(b)) shows that the instabilities in the upper part of the gap result from resonances between modes corresponding to compacton internal modes (27) and the continuous linear spectrum modes, which get coupled as the tail of the solution gets more extended. (These are seen in Fig. 9(b) as eigenvalue collisions at µ ≈ 0.005 and µ ≈ 0.015, but do not generate any instability in this figure since the corresponding eigenmodes are uncoupled in the exact compacton limit. However, they generate oscillatory instabilities when the exact compacton condition is not fulfilled, as seen in Fig. 10(a).) On the other hand, the instabilities appearing close to the lower gap edge, where the shape of the gap mode is far from compacton-like and closer to a continuum gap soliton (see Fig. 11 (a)) arise from purely imaginary eigenvalues. Direct numerical simulations of the dynamics in this regime ( Fig. 11 (b)) shows that the main outcome of these instabilities is a spatial separation of the spin-up and spin-down components.
As for the type II gap modes that arise in the mini-gap from the superposition of two in-phase neighboring single compactons in the presence of nonlinearity, we obtained pure imaginary eigenvalues in the whole mini-gap, even for the case when value of σ slightly differs from Γ (see Fig. 10 (b)). Here, with further decrease of σ, eigenvalues related to oscillatory instabilities start to occur but only in the upper half of the mini-gap.
On the other hand, instability eigenvalue spectra for type III gap solutions contain only imaginary parts of complex eigenvalues (see Fig. 10 (c)). These instabilities are always present in the lower half of the mini-gap and expand to the upper part as we move further from the compacton limit.

V. CONCLUSIONS
We derived the relevant tight-binding model for a zigzag-shaped chain of spin-orbit coupled exciton-polariton condensates, focusing on the case with basis functions of zero angular momentum and chain angles ±π/4. The simultaneous presence of spin-orbit coupling and nontrivial geometry opens up a gap in the linear dispersion relation, even in absence of external magnetic fields. At particular parameter values, where the strength of the dispersive and spin-orbit nearest-neighbor couplings are equal, the linear dispersion vanishes, leading to two flat bands with associated compact modes localized at two neigboring sites.
We analyzed, numerically and analytically, the existence and stability properties of nonlinear localized modes, as well in the semi-infinite gaps as in the mini-gap of the linear spectrum. The stability of fundamental single-peaked modes in the semi-infinite gaps was found to depend critically on the parameter a describing the relative strength of the nonlinear interaction between polaritons of opposite and identical spin (the latter assumed to be always repulsive). Generally, a spin-mixed mode with phase difference π/2 between spin-up and spin-down components is favoured when a > 1 (cross-interactions repulsive and stronger than self-interactions), while a spin-polarized mode is favoured for a < 1, which is the typical case in most physical setups. However, significant regimes of linear stability were found also for spin-mixed modes with zero phase difference between components when |a| < 1, and for spin-polarized modes when a > 1.
For parameter values yielding a flat linear band, nonlinear compactons appear in continuation of the linear compact modes, in the mini-gap as well as in the semi-infinite gaps. The linear stability eigenvalues for a single two-site compacton were obtained analytically, and shown to result in purely stable compactons inside the mini-gap when a < 1, while regimes of instability were identified in the semi-infinite gaps, and when a > 1 also inside the mini-gap. Continuing compact two-site modes away from the exact flat-band limit yields the exponentially localized fundamental FIG. 10: Imaginary parts of stability eigenvalues for the continuation of: fundamental (I type) (a), type II (b) and type III (c) gap mode inside the mini-gap, when Γ = 0.01 and a = 0.5. Pure imaginary eigenvalues are depicted with green (light gray) triangles, while red (dark gray) circles correspond to imaginary parts of complex eigenvalues. When σ = 0.01, the eigenvalues for fundamental gap mode are those of the compacton illustrated in Fig. 9(b). From bottom to top, σ is decreased to 0.007. Blue vertical dotted lines represent the locations of the lower and upper gap edges. Profiles of the corresponding solutions at σ = 0.007 in the mid-gap (µ = 0) are depicted in Fig. 4.
nonlinear gap modes inside the mini-gap. Several new regimes of instability develop, but the fundamental gap modes typically remain stable in large parts of the lower half of the gap when a < 1. We also found numerically nonlinear continuations of superpositions of two overlapping neighboring compactons (i.e., localized on three sites) with phase difference zero or π/2, where the latter also were found to exhibit significant regimes of linear stability in the mini-gap.
The model studied here may have an experimental implementation with exciton-polaritons in microcavities. Recently, microcavities have been actively investigated as quantum simulators of condensed matter systems. Polaritons have been proposed to simulate XY Hamiltonian [18], topological insulators [19], various types of lattices [20][21][22] among other interesting proposals [23] many of which were realized experimentally. In fact, the quasi one-dimensional zigzag chain considered here may be a more practical system to study the effects of interactions in presence of spinorbit coupling as compared to the full two-dimensional systems mentioned above. A possible realization of the studied system could be using microcavity pillars or tunable open-access microcavities [24]. In the latter ones, large values of TE-TM splitting can be achieved exceeding that of monolithic cavities by a factor of three [10]. Apart from directly controlling the strength of TE-TM splitting by changing parameters of the experimental system such as the offset of the frequency from the center of the stop band of the distributed Bragg reflector [25], one more possibility to control . Note the tendency for the spin-up and spin-down components to localize mainly on odd and even sites, respectively, after t ∼ 10 4 . parameters of the system is provided by using the excited states of the zigzag nodes such as spin vortices which were shown to influence the sign of the coupling strength between the sites in a polaritonic lattice [26].
Finally, we note also the recent realizations of zigzag chains with large tunability for atomic Bose-Einstein condensates [27], opening up the possibility for studying related phenomena involving spin-orbit coupling in a different context.