Modulational instability and solitary waves in polariton topological insulators

Optical microcavities supporting exciton-polariton quasi-particles offer one of the most powerful platforms for investigation of rapidly developing area of topological photonics in general, and of photonic topological insulators in particular. Energy bands of the microcavity polariton graphene are readily controlled by magnetic field and influenced by the spin-orbit coupling effects, a combination leading to formation of linear unidirectional edge states in polariton topological insulators as predicted very recently. In this work we depart from the linear limit of non-interacting polaritons and predict instabilities of the nonlinear topological edge states resulting in formation of the localized topological quasi-solitons, which are exceptionally robust and immune to backscattering wavepackets propagating along the graphene lattice edge. Our results provide a background for experimental studies of nonlinear polariton topological insulators and can influence other subareas of photonics and condensed matter physics, where nonlinearities and spin-orbit effects are often important and utilized for applications.


Introduction
Topological insulators and topologically protected edge states attract nowadays unprecedented attention in diverse areas of science, including solid-state physics, acoustics, matter waves, graphene-based applications and photonics, see, e.g., [1,2] for recent reviews. Topological insulators are characterized by the presence of the complete band-gap in the bulk of the material, like in a usual insulator, while at the same time they admit in-gap edge states propagating at the surface, where conduction of electrons becomes possible in the presence of magnetic field. These edge states are generally allowed if materials placed in contact have bulk Hamiltonians characterized by different topological invariants -Chern numbers. Apart from their conductive properties, edge states in topological insulators are immune to backscattering through their topological protection [1,2]. Topological edge states were shown to exist in HgTe quantum wells, BiSb alloys, and some other materials, where their unidirectional character and immunity to backscattering were also confirmed [1,2]. Spin-orbit interactions for electrons is the key phenomenon underpinning existence of the topological insulator phase, whose physics is closely linked with quantum Hall effect and integer Hall conductance [1,2].
Electromagnetic topological edge states have also been under intense investigation [3]. They were predicted and observed in gyromagnetic photonic crystals with pronounced Faraday effect in microwave range [4,5], in arrays of coupled resonators [6,7], and in metamaterial superlattices [8]. One of the most spectacular realizations of unidirectional edge states at optical frequencies was reported in honeycomb arrays of helical waveguides [9].
Electronic and photonic edge states mentioned above are purely linear. Though nonlinearities associated with either optical transitions or with inter-particle interactions are inherent in the majority of optical systems, the investigation of their impact on non-topological and topological edge states is in its infancy at this moment. Thus, non-topological edge states in photonic graphene have been recently studied in [10], the edge states were shown to exist in the presence of nonlinearity, while an attempt to observe soliton effects gave initial localization subsequently accompanied by noticeable radiation into the bulk, indicating coupling to the extended modes of the lattice. Nonlinearity was shown to strongly affect transmission and reflection of edge modes at the corners of graphene-like photonic structures [11]. Photonic graphene stripes of small width (ribbons), where the edge effects are mixed with the bulk dispersion, were considered in [12] and various non-topological nonlinear localized states were found, which bear more from the solitons in the bulk of the lattice.
Even more rare results on nonlinear topological states that were obtained so far are connected with scalar optical models describing evolution of excitations in arrays of twisted waveguides. Bulk nonlinear modes (i.e. modes located in the depth of periodic structure and not on its surface) of topological insulator made of helical waveguides were obtained in [13]. Traveling topological states in helical arrays were constructed in [14], but only in the frames of simplified discrete model. Dynamical excitation of their very well localized continuous counterparts considered in the unpublished Ref. [15] illustrates that longitudinal oscillations of waveguides introduce strong radiative losses leading to notable reduction of peak power already after traversing of ten sites of the structure. We should also mention here a paper linking topology and nonlinearity in a one dimensional dimer chain [16].
Phenomenology of topological insulators and edge states was recently transferred to a rapidly developing domain of exciton-polaritons in microcavities [17,18]. Main advantages of polaritons include sufficiently strong spin-orbit coupling originating in the cavity induced TE-TM splitting of the polariton energy levels [19,20], established technology of the microcavity structuring into arbitrary lattice potentials [19,21], and very strong nonlinear interactions of polaritons through their excitonic component. The latter was used for recent demonstrations of superfluidity [22,23], generation of dark quasi-solitons and vortices [24][25][26][27], bright spatial and temporal solitons [28][29][30], and other effects. The observed polariton effects with linear and nonlinear lattice potentials include one- [31] and twodimensional [32,33] gap polariton solitons, visualization of Dirac cones [34] and flat bands [35], and visualization of non-topological edge states [21]. Recently, it has been shown theoretically that attractive nonlinear interaction between polaritons with opposite spins can compensate and exceed Zeeman energy shifts due to magnetic field and thereby lead to the inversion of the propagation direction of the edge states [36]. Apart from this result the nonlinear effects with topological polariton edge states remain unexplored.
Thus, to the best of our knowledge, robust (long-living) nonlinear topological edge states confined in the direction parallel to the surface were never demonstrated in the frames of continuous physical models and in real-world systems, where unidirectionality and topological protection are provided by physical effects different from basic temporal variation of the underlying potential that always introduces undesirable losses. Moreover, compact topologically protected nonlinear edge states were never addressed in nowadays rapidly developing and open for experimental exploration spinor systems, such as spin-orbit coupled polariton and Bose-Einstein condensates, that may feature much richer dynamics and provide principally new tools for control of the state of the system in comparison with conventional scalar settings. It should be stressed that realization of topological insulators supporting surface transport of localized nonspreading excitations over considerable time intervals is a problem of fundamental physical importance, whose solution may pave the way to a number of practical applications.
We show here one such system is represented by the interacting polariton topological insulator. We use continuous model for spin-orbit coupled polariton condensate in honeycomb arrays accounting for spin-dependent interactions and Zeeman splitting in the external magnetic field to demonstrate a variety of new nonlinear topologically protected edge states and to perform their rigorous stability analysis, for the first time for this class of topologically protected nonlinear modes. Extended periodic edge states are found in the exact form as truly stationary nonlinear solutions bifurcating from linear periodic edge states. We show that their modulation instability results in splitting of the extended nonlinear edge states into sets of fully localized edge quasi-solitons travelling along the interface over very long time intervals without notable deformations. Such localized quasisolitons appear to be very robust objects on any practical time scales. Even though they radiate as they move along the interface, the rate of radiation can be exceptionally small in the proper range of parameters ensuring that such states can traverse huge distances along the surface of the material without being destroyed or scattered. The approximate expression for the shape of quasi-solitons is derived.

Topological edge states in the linear regime
Polariton condensate in the lattice potential in the presence of the external magnetic field can be described by a system of coupled Gross-Pitaevskii equations for the spin-positive and spin-negative components y  of the spinor polariton wavefunction Quasi-conservative nonlinear dynamics has been observed in several experiments with exciton polaritons, see, e.g. [19,24,30,31], and used in many theoretical studies, see, e.g. [17][18][19]37]. Following this trend we have also chosen to work here in the idealized conservative limit, since the very fact of existence of unidirectional edge states is not connected with presence of losses, so we do not include them to have most transparent picture of the phenomenon. Here the relation between the wavefunction components in the circular polarization basis and those associated with TE (subscript y ) and TM (subscript x ) polarizations is given by The b -term describes spin-orbit coupling originating from the TE-TM energy splitting of the cavity photons, which translates into slightly different effective masses , -. Accounting for the spinor nature of the system and transforming into the basis of circular polarizations results in the spinorbit term raised to the second power, see, e.g., [37]. Parameter z e is the Zeeman energy splitting of spin + and spin -polaritons proportional to the externally applied magnetic field. 0 ( , ) x y e  describes potential energy landscape in the microcavity [19,21,32,33]. In our case ( , ) x y  is a honeycomb lattice that is cut in the x direction. The distance between the lattice sites is a , so that . 0 e is the scaling coefficient with the dimension of energy, while  is the dimensionless function. Amplitudes y  can be assumed dimensionless and scaled in a way that the nonlinear energy shift achieved for the unit polariton density equals 0 e . This scaling well reflects the physically realistic situation, when energies, associated with the lattice potential, Zeeman effect and nonlinearity have the same order of magnitude [17][18][19].
0.05 s = -is the strength of the weak attractive interaction of polaritons with the opposite spins [38]. Local potential minima in ( , ) x y  are described by Gaussian functions (unless stated otherwise), so that both of them are an order of magnitude less than the energy shift induced by the lattice. Note, that the parabolic approximation for the polariton energy momentum dependence dictated by the Gross-Pitaevskii approximation is well obeyed providing that all other energy shifts in our model are less than the width of the low polariton branch by a factor of 2-3 or more, which is well satisfied if we realistically assume that the latter is 10 meV  wide. With these normalizations the dimensionless version of Eq. (1) can be written as where we retained old notations for scaled coordinates , x y and evolution time t .
For these parameters and for , 0 b W ¹ the simplest linear mode of an isolated local minimum of the potential is characterized by the presence of vortex with topological charge 2 in y + component with ring-like density profile and by trivial constant phase of bell-shaped ycomponent. Inter-estingly, next mode carries vortex in the ycomponent, while y + component has trivial phase. The appearance of charge-2 vortex in one of the components is a consequence of spin-orbit coupling and it can be observed only for 0 b ¹ . If the potential has two well-separated minima, the above mentioned charge-2 vortices carried by one of the components and located around each potential minima split into two charge-1 vortices. When two potential minima become very close, one observes linear modes with complex phase distributions with more than 4 phase singularities in one of the components.
. Left and right columns show y + and y -, respectively.
(c) and for the lattice with bearded edges at Black lines correspond to modes residing in the bulk of the lattice, while color lines indicate edge states. Topologically protected edge states in (b)-(d) with opposite slopes of ( ) k e belong to the opposite edges of the lattice. In all cases 0.5 W = .
There are three types of termination of extended honeycomb arrays corresponding to the zigzag (Fig. 1, top left), bearded (Fig. 1, top right), and armchair edges. The armchair cut usually does not support the edge states [9], so we do not consider it here. We first address the spectrum of linear modes that are periodic along the array edge and localized along the xaxis on both sides from the edge. These modes are the Bloch functions 3 a is the y -period of the potential. The corresponding momentum giving the width of the Brillouin zone is . We calculated Bloch functions numerically using a unit cell containing 36 potential minima (top row of Fig. 1 depicts three such unit cells, i.e. three periods of potential along the y -axis). Representative spectra for the lattice with zigzag and bearded edges are shown in Fig. 2 in the form of the energy- Due to spinor character of the model the spectrum consists of two groups of bands. At , 0 b W = these two groups are fully degenerate, and correspond to the Bloch modes with The inclusion of Zeeman splitting lifts this degeneracy and results in vertical splitting of two energy bands by 2W , clearly visible in Fig. 2(a) for zigzag edge. The spectrum in Fig. 2(a) shows two Dirac points at K/ 3 k =  , where first and second bands in each group touch each other. These points are traces of Dirac points in the spectrum of bulk honeycomb lattice. Red branches in Fig. 2(a) correspond to the non-topological edge states appearing due to truncation of the array, while black branches correspond to modes concentrated in the bulk. In the bearded edge case the edge modes appear in the region [ K/3, K/3] k Î -+ between two Dirac points, while for zigzag edge they appear outside this domain, at This picture qualitatively changes when spin-orbit coupling is accounted for, see Figs. 2(b), 2(c) obtained for zigzag edge. It leads to opening of the full gap between the first (upper most) and second bands, and eliminates Dirac points. The width of the gap increases with b . Most importantly, spin-orbit coupling removes degeneracy of the edge states, so that the states residing on the opposite sides of the array acquire opposite group velocities / k e ¶ ¶ , associated unidirectionality, and topological protection features [17]. The degeneracy is removed and edge states become unidirectional only for nonzero values of b (the effect is observable already at 0.01 b  ). Thus, spin-orbit coupling is an absolutely necessary ingredient for observation of all effects associated with unidirectionality that are mentioned below. It is SO-coupling that leads to nonzero Chern number and topological character of unidirectional edge states, both linear ones and bifurcating from them nonlinear modes. Moreover, spin-orbit coupling leads to specific band folding [ Fig. 2(b)] accompanied by opening of the additional lower energy gaps, where topologically protected edge states appear in the intervals of momenta roughly complimentary to the ones where the edge states of the primary gap exist. To the best of our knowledge, the counterparts of these states were not encountered in scalar (non-spinor) honeycomb lattice models with zigzag termination [21] and in the Floquet topological insulators formed by helical waveguides [9].  . Red and green circles correspond to linear modes depicted in Fig. 1. Line colors correspond to colors used in Fig. 2(c) to denote different branches.
Examples of the profiles for the edge states in the top and bottom gaps are shown in Fig. 1 for the same parameters as used in Fig. 2(c) (zigzag edges). Edge states extend into the bulk of the lattice when momentum k approaches the boundary of the existence domain (second and fourth rows), but are well localized for k values returning energies close to the center of the gap (third row). The modes residing at the opposite edges, but having the same energies for different momenta are "mirror images" of each other (compare fifth and second rows). Edge states from the bottom gap associated with magenta line in Fig. 2(c) also can be well-localized (see sixth row in Fig. 1) provided that their energy is far from the gap edges. Note, that the ( ) k e plots for the topological edge states in our system are not antisymmetric with respect to the K/ 2 k = line, as in arrays of twisted waveguides [9]. The degree of asymmetry is controlled by the spin-orbit coupling. The asymmetry of the existence intervals with respect to the K/ 2 k = point and stronger localization around this point can be well seen in the dependence of the edge state width x w on the momentum [ Fig. 3(a)]. The width is defined as are the integral form-factors and the norms for the two spin components defined on one y -period of the structure, respectively.

Nonlinear topological edge states and their modulational instability
Solutions accounting for nonlinearity are sought in the same form as the linear ones ) .
We solved Eq.  Stability analysis of nonlinear edge states ( , ) u x y  was performed by perturbing them with small (1% in amplitude) broadband input noise and direct modeling of their subsequent evolution up to very large times (usually up to ). Nonlinear edge states were found to be unstable (recall that in the presence of periodic potentials modulation instability is possible even in the medium with defocusing nonlinearity [39,40]). The instability is particularly pronounced close to the edge of the band and is strongly suppressed when one approaches a point, where the nonlinear edge states bifurcate from the linear spectrum. Only low-frequency perturbations with huge y -periods substantially exceeding separation between pillars can destabilize such modes and their growth rates are so small that the instability manifests itself at times 3 10 t > far exceeding lifetime of polaritons, so such modes will appear as stable ones in experiments. Modulation instability bandwidth (i.e. the range of frequency of modulations along the y axis that can seed instability) becomes smaller with increase of m . Such bandwidth can be calculated from Eq. (2) using initial conditions ( 0) ( , )[1 cos( )]exp( ) t u xy y i k y y n k   = = + , where n and k are the amplitude and frequency of small perturbations. Such perturbations experience clear exponential growth at the initial stage of instability development as long as k is within modulation instability band, that allows to determine instability growth rate d as a function of k . Fig. 5(d) shows representative ( ) d k dependence that reveals the finite instability bandwidth. Typical dynamics of instability development stimulated by broadband noise is shown in Fig. 6. Instability initially leads to pronounced modulation of the edge state in the y -direction. As time goes this modulation is becoming accompanied by radiation into the depth of the array. However, instead of decay and disintegration, the edge state breaks up into a set of weakly radiating, but fully localized solitons (see the pattern at 960 t = in Fig. 6). This suggests that predominantly repulsive nonlinearity and

Edge quasi-solitons
To develop more regular approach to the question of quasi-soliton existence we rewrite Eq. (2) in the form  and therefore describes spatial shape of the linear Bloch mode with momentum k , and we also take into account that corresponding energy e also depends on quasi-momentum k . Here k is the momentum offset from the carrier soliton momentum k and the amplitude ( , ) A t k is assumed well localized in k . Using Taylor series expansion in k for ( , , ) x y k k +  in the above integral one arrives at the expression for the shape of the edge state wavepacket: We further assume that the spinor  changes with k much slower than the eigenvalue e , that is a valid assumption as simulations show. This allows to keep only 0 n = term in (5), so that . Finally, we multiply the equation  by †  from the left and integrate it over one period along the y -axis and along the entire x -axis, assuming slow variation of the envelope function ( , ) A y t with y , that allows to remove it from all integrands. This yields the nonlinear Schrödinger equation for the envelope function: where we kept only first two terms proportional to / k e e ¢ = ¶ ¶ and where 0 m e -£ is the energy shift due to repulsive nonlinearity. Note that energy shift in (8)   . Middle row shows peak amplitude aof ycomponent for linear and nonlinear cases. ydistributions in the top and bottom rows correspond to red circles. Notice that edge states in all panels in this figure move upwards. To simultaneously stress the fact that edge state moves in one direction and to provide details of its shape, we show ydistributions within relatively small vertical window, but we also applied identical vertical shift for distributions at time moments 420, 540, 660 t = , where edge state in nonlinear medium remains nearly invariable. The distribution at 0 t = was not shifted. Here Figure 3(b) shows dispersion coefficient e¢¢ as a function of k for different linear edge states: we use the same colors as in Fig. 2(c) to denote different branches. One can see that there exist momentum intervals where dispersion coefficient e¢¢ is positive (effective polariton mass is negative). All such intervals for every linear edge mode give rise to unidirectional edge quasi-solitons. Note, that for every branch there exist a unique k value where dispersion e¢¢ vanishes and where excitation with broad envelope may evolve almost without broadening even in the linear limit [shape distortions in this case will be determined by weaker e¢¢¢ dispersion that was omitted in (7)].
Top row of Fig. 7 shows the central result of this work -evolution of the quasi-soliton constructed using Eq. (5) with the envelope function given by Eq. (8) in the nonlinear topological insulator state. This quasisoliton corresponds to the red branch in Fig. 2(c) and 0.55 K k = . This k value was selected approximately in the middle of the domain (in terms of k ) where dispersion coefficient e¢¢ for linear edge modes from red branch of Fig. 2(c) is positive. Indeed, this branch exists at 0.24 / K 0.62 k < < , while dispersion changes its sign at 0.44 K k » .
The particular momentum value 0.55 K k = at which we demonstrate edge soliton formation is not preferred in comparison with other k values at which 0 e¢¢ > : we were able to generate similar long-living nonlinear edge states for multiple values of k .
One can see that after the initial transient where the peak amplitude decreases due to internal reshaping of the input, the unidirectional quasisoliton forms whose amplitude remains almost constant in time and whose velocity nearly coincides with e¢ (see red curve in the central row with dependence of peak amplitude nlin aof ycomponent). Note, that, during time period shown, the quasi-soliton traverses over 100 periods of the array. One can observe small radiation into the depth of array that gradually reduces soliton amplitude (that is why we call such states quasisolitons), but even on the time scales it is a negligible effect, so that quasi-solitons found here are exceptionally robust at any time scales. Similar unidirectional states were found for other dispersion branches (thus counterparts of soliton from Fig. 7, but from green branch move in the opposite y -direction) and also in arrays with bearded edges.  To confirm that the obtained localized states indeed exist due to nonlinear self-action, we used the same input, but switched off nonlinearity, see the bottom row in Fig. 7. Without nonlinearity we have observed a pronounced asymmetric expansion of the wavepacket, while the peak amplitude lin aof ycomponent was strongly decreasing upon evolution (see black curve in the middle row). We have also checked that linear edge states with 0 e¢¢ < do not give rise to quasi-solitons. To confirm exceptional robustness of quasi-solitons and their immunity to backscattering in nonlinear regime we considered their evolution in the lattice potential that was made finite also along the y -axis (Fig. 8). We used the same quasi-soliton input as in Fig. 7. One can see that the soliton survives even upon passage of several array corners and that it returns to its initial location after making a closed loop along the surface of the array with minimal decrease in peak amplitude of both components. Note, that the radiation into the depth of array is pronounced only along the top and bottom armchair edges that formally do not support edge states and that is why they were made relatively short.
Finally, the absence of backscattering on edge defects is illustrated in Fig. 9, where we removed one of the micropillars on the surface of the array. Before collision with this defect, the quasi-soliton was allowed to evolve over sufficiently long time to ensure that its amplitude has reached a "steady state" value. Note, that the soliton experiences considerable reshaping at the point of collision with the defect, so that its amplitude drops nearly by a factor of 2, but it returns to its initial value immediately after passage of the defect.

Summary
Summarizing, we predicted formation of extended nonlinear edge states in polariton condensates with spin-orbit coupling held in the honeycomb lattice potentials. Nonlinear edge states spontaneously decay into sets of fully localized quasi-solitons that can also be selectively excited by using proper envelope function, derived in this paper. The edge solitons have been shown to be robust with respect to the passage through the intervals of the lattice edges that do not support the edge states in the linear approximation and through other lattice defects.