# Stability of solitary waves and vortices in a 2D nonlinear Dirac model

###### Abstract

We explore a prototypical two-dimensional massive model of the nonlinear Dirac type and examine its solitary wave and vortex solutions. In addition to identifying the stationary states, we provide a systematic spectral stability analysis, illustrating the potential of spinor solutions to be neutrally stable in a wide parametric interval of frequencies. Solutions of higher vorticity are generically unstable and split into lower charge vortices in a way that preserves the total vorticity. These conclusions are found not to be restricted to the case of cubic two-dimensional nonlinearities but are found to be extended to the case of quintic nonlinearity, as well as to that of three spatial dimensions. Our results also reveal nontrivial differences with respect to the better understood non-relativistic analogue of the model, namely the nonlinear Schrödinger equation.

Introduction. In the context of dispersive nonlinear wave equations, admittedly the prototypical model that has attracted a wide range of attention in optics, atomic physics, fluid mechanics, condensed matter and mathematical physics is the nonlinear Schrödinger equation (NLS) kivshar (); becbook1 (); becbook2 (); rcg:BEC_book2 (); ablowitz1 (); sulem (); chap01:bourgain (). By comparison, far less attention has been paid to its relativistic analogue, the nonlinear Dirac equation (NLD) nogami (), despite its presence for almost 80 years in the context of high-energy physics jetp.8.260 (); PhysRev.83.326 (); RevModPhys.29.269 (); thirring (); gavri (). This trend is slowly starting to change, arguably, for three principal reasons. Firstly, significant steps have been taken in the nonlinear analysis of stability of such models boucuc (); DEP1 (); DEP2 (); DEP3 (); arXiv:1407.0606 (); linear-a (), especially in the one-dimensional (1d) setting. Secondly, computational advances have enabled a better understanding of the associated solutions and their dynamics sihong_rev (); cooper (); niur_recent (); MR2892774 (); avadhj (). Thirdly, and perhaps most importantly, NLD starts emerging in physical systems which arise in a diverse set of contexts of considerable interest. These contexts include, in particular, bosonic evolution in honeycomb lattices Carr1 (); Carr2 () and a growing class of atomically thin 2d Dirac materials diracmat () such as graphene, silicene, germanene and transition metal dichalcogenides tmdc (). Recently, the physical aspects of nonlinear optics, such as light propagation in honeycomb photorefractive lattices (the so-called photonic graphene) ablowitz3 (); ablowitz4 () have prompted the consideration of intriguing dynamical features, e.g. conical diffraction in 2d honeycomb lattices peleg (). Inclusion of nonlinearity is then quite natural in these models, although in a number of them (e.g., in atomic and optical physics) the nonlinearity does not couple the spinor components.

These physical aspects have also led to a discussion of potential 2d solutions of NLD in Carr1 (); Carr2 (). However, a systematic and definitive characterization of stability and nonlinear dynamical evolution of the prototypical coherent structures in NLD models is still lacking, to the best of our knowledge. The present work is dedicated to offering analytical and numerical insights into these crucial mathematical and physical aspects of higher-dimensional nonlinear Dirac equations bearing in mind the physical relevance and potential observability of such waveforms. As our model of choice, in order to also be able to compare and contrast with the multitude of existing 1d results (e.g. arXiv:1407.0606 (); MR2892774 ()), we select the well-established Soler model Soler () (known in 1d as the Gross–Neveu model Gross1 ()), which is a Dirac equation with scalar self-interaction. Such self-interaction is based on including into the Lagrangian density a function of the quantity (which transforms as a scalar under the Lorentz transformations):

(1) |

where , , , and , are Dirac -matrices satisfying the anticommutation relations , with the Minkowski tensor dewitt (), and . (The Clifford Algebra theory gives the relation between the spatial dimension and spinor components (MR1401125, , Chapter 1, §5.3).) The nonlinearity of the model is generalized in the spirit of cooper (), by using with . The proof of existence of solitary waves in this model (in 3d) is in MR847126 (); MR949625 (); MR1344729 ().

Our results show that the NLD in 2d admits different solutions involving a structure of vorticity in the first spinor component, with the other spinor component bearing a vorticity . We identify such solutions for . While prior stability results have often been inconclusive (particularly in higher dimensions, see, e.g., PhysRevLett.50.1230 ()), our numerical computation of the spectrum of the corresponding linearization operator reveals that only the solutions can be spectrally stable (the spectrum of the linearization contains no eigenvalues with positive real part), and that this stability takes place in a rather wide interval of the frequency of the solitary waves. On the contrary, we find that the states of higher vorticity are generically linearly unstable. Complementing the stability analysis results, our direct dynamical evolution studies show that the unstable higher vorticity solutions break up into lower vorticity waveforms, yet conserving the total vorticity. Importantly, the fundamental solutions are found to be potentially stable in models both with a higher order (quintic) two-dimensional nonlinearity, as well as in higher dimensions (3d) under cubic nonlinearity. These features again reflect differences from the NLS model and as such suggest the particular interest towards a broader and deeper study of NLD models.

An important extension of our stability findings for higher dimensional solutions, is that they remain valid for other types of nonlinearities. These include non-Lorentz-invariant ones such as most notably those arising in atomic Carr1 (); Carr2 (), and optical ablowitz3 (); ablowitz4 () problems. The fundamental difference of those models is that they correspond to massless equations, contrary to the Soler model. For this reason, we have confirmed our stability conclusions by comparison with those emerging from the model for square binary waveguides Biancalana () which leads to a massive nonlinear Dirac equation with the same nonlinearity as in Carr1 (); Carr2 ().

Theoretical Setup. We start from the prototypical 2d nonlinear Dirac equation system, derived from the Lagrangian density (1) with and :

(2) |

where , are the components of the spinor and the nonlinearity is . We note that (Stability of solitary waves and vortices in a 2D nonlinear Dirac model) is a -invariant, translation-invariant Hamiltonian system.

We simplify our analysis by using the polar coordinates, where Eq. (Stability of solitary waves and vortices in a 2D nonlinear Dirac model) takes the form

(3) |

The form of this equation suggests that we look for solutions as with

(4) |

with and real-valued. The value can be cast as the vorticity of the first spinor component.

Once solitary waves have been identified, we explore their stability. This approach has been previously developed in related settings including the multi-component NLS (see e.g. mussli ()), as well as a massless variant of the Dirac equation of hc_rlse (). The presence of a mass in our case allows not only a direct comparison with NLS (when ), but also generates fundamental differences between our results and those of Carr1 (); Carr2 (); hc_rlse (), as discussed below as well.

To examine its spectral stability, we consider a solution in the form of a perturbed solitary wave solution:

(5) |

with a small perturbation. We consider the linearized equation on ,

(6) |

with and with a matrix-valued first order differential operator Supp (). If the spectrum of the linearization operator contains an eigenvalue with , we say that the solitary wave is linearly unstable; in such cases, we resort to dynamical simulations of Eqs. (Stability of solitary waves and vortices in a 2D nonlinear Dirac model) to explore the outcome of the unstable evolution. If there are no such eigenvalues, the solitary wave is called spectrally stable.

A convenient feature of NLS ground states is that the linearization operator at such states, albeit non-selfadjoint, has its point spectrum confined to the real and imaginary axes. This observation is at the base of the VK criterion VaKo (): a linear instability can thus develop when a positive eigenvalue bifurcates from . More precisely, the loss of stability due to the appearance of a pair of positive and a pair of negative eigenvalues follows the jump in size of the Jordan block corresponding to the unitary invariance; this happens when the VK condition is satisfied, with being the charge of a solitary wave.

Crucially, in the NLD case, the spectrum of the linearization at a solitary wave is no longer confined to the real and imaginary axes; the linear stability analysis requires that one studies the whole complex plane. The key observation is that in (6) contains , , , but not ; this allows to perform a detailed study of the spectrum of using the decomposition of spinors into Fourier harmonics corresponding to different Supp ().

In the 3d case, we are not yet able to perform the general spectral analysis, but we studied the part of spectrum in the invariant subspace corresponding to perturbations of the same angular structure as the solitary waves wakano-1966 (), ; this invariant subspace seems most important since it is responsible for the linear instability in the non-relativistic limit which is a consequence of the instability of the 3d cubic NLS.

Numerical results. We have analyzed the existence and stability of solitary waves (, with its first component radially symmetric and the second component having vorticity ) and vortex solutions (, with its components having vortices of order one and two, respectively). Both solitary waves and vortex solutions exist in the frequency interval , a feature critically distinguishing our models from those of Carr1 (); Carr2 (). An intriguing feature of the relevant waveforms is that both the radial profile of the solitary waves and that of the vortices possess a maximum that shifts from to a larger when approaches zero (see Fig. 1), in a way reminiscent of the corresponding 1d solitary wave structures cooper (). Here the relevant state will feature a stationary bright intensity ring. In order to obtain and analyze such coherent structures, we have made use of the numerical methods detailed in Supp (). To confirm the results, we also computed the spectra using the Evans function approach of MR2892774 () adapted to the present problem.

We start by considering the stability of solitary waves in the cubic () case. Figure 2 shows the dependence of the real and imaginary parts of the eigenvalues with respect to the stationary solution frequency . From the spectral dependencies we can deduce several features of the 2d NLD equation: (1) It is known that the 2d NLS equation is charge-critical, and the zero eigenvalues are degenerate sulem (): they have higher algebraic multiplicity. In the NLD case, however, this degeneracy is resolved: in the case, as starts decreasing, two eigenvalues (corresponding to ) start at the origin when and move out of the origin for . The absence of the algebraic degeneracy of the zero eigenvalue prevents solitary waves from NLS-like self-similar blow-up which is possible in charge-critical NLS MR1048692 (). (2) The symmetry and the translation symmetry of the model result in zero eigenvalues with and , respectively (in both and cases). (3) As in the 1d NLD equation, there are also the eigenvalues which are associated with the symmetry of the model ptnlde (). This eigenvalue pair corresponds to , i.e., to a highly excited linearization eigenstate. (4) Contrary to the 1d case, where the solitary waves corresponding to any are spectrally stable, the solitary wave is linearly unstable for because of the emergence of nonzero real part eigenvalues via a Hamiltonian Hopf bifurcation in the spectrum at . Another Hopf bifurcation occurs corresponding to (at ), then yet another one corresponding to .

It is especially interesting that a wide parametric (over frequencies) interval of stability of solitary waves can also be observed in the quintic () NLD case (see Fig. 3); while the quintic NLS solitary waves blow up (even in one dimension), the quintic NLD solitary waves are stable even in two dimensions, except for the interval where the coherent structures experience the same Hopf bifurcation as in the cubic case, and for where an exponential instability created by radial perturbations emerges. Perhaps even more remarkably, the right panel of the Fig. 3 illustrates that this stability of NLD solitons against radial perturbations can be found in suitable frequency intervals even in 3d (see CGG14 () for a discussion of the equations for existence and stability of radial perturbations in 3d). Both of the above cases (quintic 2d and cubic 3d NLD) are charge-supercritical i.e., the charge goes to infinity in the nonrelativistic limit . Contrary to the pure-power supercritical NLS whose solitary waves remain linearly unstable for all frequencies, solitary waves in the Soler model become spectrally stable when drops below some dimension-dependent critical value Supp (). The relevant unstable eigenvalue (associated with and radially-symmetric collapse) is only present as real for , where . This was identified in Soler () as the value at which both the energy and charge of solitary waves have a minimum. Hence, we indeed find that the radially-symmetric collapse-related instability ceases to be present below this critical point. Finally, as regards two dimensions, vortices are unstable for every , because of the presence in the spectrum of quadruplets of complex eigenvalues. These quadruplets emerge (and disappear) for different values of via direct (inverse Hopf) bifurcations; see the right panel of Fig. 2. The spectrum for vortex is quite similar to that of ; for this reason, we do not analyze it further.

In order to analyze the result of instabilities in 2d settings, we have probed the dynamics of unstable solutions directly (see Supp () for details). Prototypical examples of unstable solitary waves and vortices for are shown in Figs. 4 and 5. As can be observed, the solitary waves spontaneously amplify perturbations breaking the radial symmetry in their density and, as a result, become elliptical and rotate around the center of the circular density of the original solitary wave in line with the expected amplification of the unstable eigenmode. On the other hand, the vortices split into three smaller ones. Let us mention that in the latter case, the first spinor component splits into structures without angular dependence, whereas the second component splits into corresponding ones with angular dependence , in accordance with the ansatz of Eq. (4). This preserves the total vorticity across the two components, as is also shown in Fig. 5. Along a similar vein, the instability of an vortex eventually leads to the emergence of five pairs, again preserving the total vorticity. Finally, we have analyzed the outcome of the instabilities caused by radially-symmetric perturbations in the case for (see Fig. 5). We can observe the typical behavior of such solutions, i.e. the density width (and amplitude) oscillate leading to a “breathing” structure, but there is no collapse. This phenomenology is reminiscent of the 1d case nldebook ().

Conclusions and Future Challenges. We have illustrated that solitary waves of vorticity in one spinor component and in the other are spectrally stable within a large parametric interval, suggesting their physical relevance. In that connection, we highlight that although our models of choice may bear a particular nonlinearity, our results suggest that under different nonlinearities including the more physically relevant ones of e.g., Carr1 (); Carr2 () and massive models Biancalana () still bear stable solitary waves for a suitable wide parametric range of frequencies. Thus, the conclusion of higher dimensional stability is more general than the specifics of our particular nonlinearity and hence of broad interest. We also showcased the significant difference of NLD from the focusing NLS equation, where solitary waves are linearly unstable in the charge-supercritical cases. When the NLD solutions were found to be unstable, their dynamical evolution suggested breathing oscillations in the case and splitting into lower charge configurations for and .

It is of interest to extend present considerations to numerous settings. From a mathematical physics perspective, it would be useful to explore further the 3d stability and associated dynamics. This is especially timely given that the 3d analogue of photonic graphene has been experimentally realized very recently soljacic (). Admittedly, the latter setting does not feature a mass in the model, thus the generalization of NLD models such as those appearing in the works of Carr1 (); Carr2 (); hc_rlse () would be particularly important there. It would also be of interest to compare more systematically the present findings with models associated with different nonlinearities, including the case of honeycomb lattices in atomic and optical media or, e.g., those stemming from wave resonances in low-contrast photonic crystals agueev ().

###### Acknowledgements.

The research of A. Comech was carried out at the Institute for Information Transmission Problems of the Russian Academy of Sciences at the expense of the Russian Foundation for Sciences (project 14-50-00150). A.S. was supported by the U. S. Department of Energy. P.G.K. gratefully acknowledges the support of the US-NSF under the award DMS-1312856, and of the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096).## References

- (1) Yu.S. Kivshar and G.P. Agrawal, Optical solitons: from fibers to photonic crystals, Academic Press (San Diego, 2003).
- (2) L.P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, Oxford University Press (Oxford, 2003).
- (3) C.J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press (Cambridge, 2002).
- (4) P.G. Kevrekidis, D.J. Frantzeskakis, and R. Carretero-González, The defocusing nonlinear Schrödinger equation: from dark solitons and vortices to vortex rings. SIAM, Philadelphia (2015).
- (5) M.J. Ablowitz, B. Prinari, and A.D. Trubatch, Discrete and Continuous Nonlinear Schrödinger Systems, Cambridge University Press (Cambridge, 2004).
- (6) C. Sulem and P.L. Sulem, The Nonlinear Schrödinger Equation, Springer-Verlag (New York, 1999).
- (7) J. Bourgain, Global Solutions of Nonlinear Schrödinger Equations, American Mathematical Society (Providence, 1999).
- (8) F.M. Toyama, Y. Hosono, B. Ilyas, and Y. Nogami, J. Phys. A: Math. Gen. 27, 3139 (1994).
- (9) D. D. Ivanenko, Zh. Éksp. Teor. Fiz 8, 260 (1938).
- (10) R. Finkelstein, R. LeLevier, and M. Ruderman, Phys. Rev. 83, 326 (1951).
- (11) W. Heisenberg, Rev. Mod. Phys. 29, 269 (1957).
- (12) W. Thirring, Ann. Phys. 3, 91 (1958).
- (13) S.Y. Lee, T. K. Kuo, and A Gavrielides, Phys. Rev. D 12, 2249 (1975).
- (14) N. Boussaïd and S. Cuccagna, Comm. PDE 37, 1001 (2012).
- (15) D.E. Pelinovsky and Y. Shimabukuro, Lett. Math. Phys. 104, 21 (2014).
- (16) A. Contreras, D.E. Pelinovsky and Y. Shimabukuro, Comm. PDE 41, 227 (2016).
- (17) D.E. Pelinovsky and A. Stefanov, J. Math. Phys. 53, 073705 (2012).
- (18) A. Comech, T. Phan, and A. Stefanov, arXiv:1407.0606.
- (19) N. Boussaïd and A. Comech, arXiv:1211.3336.
- (20) J. Xu, S. Shao and H. Tang, J. Comp. Phys 245, 131 (2013).
- (21) S. Shao, N.R. Quintero, F.G. Mertens, F. Cooper, A. Khare, and A. Saxena, Phys. Rev. E 90, 032915 (2014).
- (22) F. Cooper, A. Khare, B. Mihaila, and A. Saxena, Phys. Rev. E 82, 036604 (2010).
- (23) G. Berkolaiko and A. Comech, Math. Model. Nat. Phenom. 7, 13 (2012).
- (24) J. Cuevas-Maraver, P.G. Kevrekidis, and A. Saxena, J. Phys. A 48, 055204 (2015).
- (25) L. Haddad, K. O’Hara, and L.D. Carr, Phys. Rev. A 91, 043609 (2015).
- (26) L. Haddad and L.D. Carr. New J. Phys. 17, 113011 (2015).
- (27) T. O. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Adv. Phys. 63, 1 (2014).
- (28) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
- (29) M. J. Ablowitz, S.D. Nixon, and Y. Zhu, Phys. Rev. A 79, 053830 (2009).
- (30) M. J. Ablowitz and Y. Zhu, Phys. Rev. A 82, 013840 (2010).
- (31) O. Peleg, G. Bartal, B. Freedman, O. Manela, M. Segev, and D.N. Christodoulides, Phys. Rev. Lett. 98, 103901 (2007).
- (32) M. Soler, Phys. Rev. D 1, 2766 (1970).
- (33) D.J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974).
- (34) B. de Wit and J. Smith, Field theory in particle physics, North Holland Physics Publishing (New York, 1986).
- (35) B.V. Fedosov, Index theorems, in: Partial differential equations, VIII Encyclopaedia Math. Sci., 65, 155. Springer Verlag (Berlin, 1996).
- (36) T. Cazenave and L. Vázquez, Comm. Math. Phys. 105, 35 (1986).
- (37) F. Merle, J. Differential Equations 74, 50 (1988).
- (38) M.J. Esteban and E. Sere, Commun. Math. Phys. 171, 323 (1995).
- (39) A. Alvarez and M. Soler, Phys. Rev. Lett. 50, 1230 (1983).
- (40) N.G. Vakhitov and A.A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973).
- (41) T.X. Tran, X.N. Nguyen, and F. Biancalana, Phys. Rev. A 91, 023814 (2015).
- (42) Z.H. Musslimani, M. Segev, D.N. Christodoulides, and M. Soljačić, Phys. Rev. Lett. 84, 1164 (2000); J. Yang and D.E. Pelinovsky Phys. Rev. E 67, 016608 (2003).
- (43) L. H. Haddad and L.D. Carr, EPL 94, 56002 (2011); L. H. Haddad and L.D. Carr, New J. Phys. 17, 063034 (2015).
- (44) See the Supplementary Material [url], which includes Refs. Boyd (), Trefethen (), Herring () numrec ().
- (45) M. Wakano, Progr. Theoret. Phys. 35, 1117 (1966).
- (46) F. Merle, Comm. Math. Phys. 129, 223 (1990).
- (47) J. Cuevas-Maraver, P.G. Kevrekidis, A. Saxena, F. Cooper, A. Khare, A. Comech, and C.M. Bender, J. Sel. Top. Quant. Electr. 22, 5000109 (2016).
- (48) A. Comech, M. Guan and S. Gustafson, Ann. I. H. PoincarÃ© â AN 31, 639 (2014).
- (49) J. Cuevas-Maraver, P.G. Kevrekidis, A. Saxena, F. Cooper and F.G. Mertens, in: Ordinary and Partial Differential Equations, Chapter 4. Nova Science Publishers (New York, 2015).
- (50) L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J.D. Joannopoulos, M. Soljac̆ic, Science 349, 622 (2015).
- (51) D. Agueev and D. Pelinovsky, SIAM J. Appl. Math. 65, 1101 (2005).
- (52) J.P. Boyd, Chebyshev and Fourier spectral methods, 2nd. Edition. Dover (2001).
- (53) L.N. Trefethen, Spectral methods in MATLAB. SIAM, Philadelphia (2000).
- (54) G. Herring, L.D. Carr, R. Carretero-González, P.G. Kevrekidis, and D.J. Frantzeskakis. Phys. Rev. A 77, 023625 (2008).
- (55) W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd. Edition. Cambridge University Press (2007).

Stability of solitary waves and vortices in a 2D nonlinear Dirac model. Supplementary material

## Appendix A Equivalence between the Soler and massive Thirring models in 2D

An alternative model is the massive Thirring model (MTM) with vector self-interaction Supp_thirring (), where the nonlinear term in the Lagrangian is based on the scalar quantity built from the Lorentz vector which represents the charge-current density. We mention that the MTM in one spatial dimension is completely integrable, and in two spatial dimensions is equivalent to the Soler model.

We demonstrate below that in two spatial dimensions the Soler and the MTM models coincide.

In the massive Thirring model with vector self-interaction, the nonlinear term in the Lagrangian is based on the scalar quantity built from the Lorentz vector which represents the charge-current density.

Using the Dirac matrices , , , and using the identity

valid for any , we compute that in 2d these two scalar quantities coincide:

## Appendix B Numerical methods for the study of the existence of stationary solutions

### b.1 Brief summary of spectral methods

Prior to explaining the numerical methods used for calculating stationary solutions, we will proceed to present a summary of spectral methods needed for dealing with derivatives in continuum settings. For a detailed discussion on these methods, the reader is directed to Supp_Boyd () and references therein.

Spectral methods arise due to the necessity of calculating spatial derivatives with higher accuracy than that given by finite difference methods. To this aim, a differentiation matrix must be given together with collocation (i.e. grid) points , which are not necessarily equi-spaced. Thus, if the spectral derivative of a function needs to be calculated, it can be cast as:

(S1) |

where and . If and the boundary conditions are periodic, the Fourier collocation can be used. In this case,

(S2) |

with even. The differentiation matrix is

(S3) |

Notice that doing the multiplication is equivalent to performing the following pair of Discrete Fourier Transform applications:

(S4) |

with and denoting, respectively, the direct and inverse discrete Fourier transform Supp_Trefethen (). The vector wavenumber is defined as:

(S5) |

The computation of the direct and inverse discrete Fourier transforms, which is useful in simulations, can be accomplished by the Fast Fourier Transform. In what follows, however, the differentiation matrix is used for finding the Jacobian and stability matrices. Notice that the grid for finite differences discretization is the same as in the Fourier collocation; and, in addition, there is a differentiation matrix for the finite differences method, i.e.

(S6) |

with being Kronecker’s delta. It can be observed from the above discussion that in the Fourier spectral method, the banded differentiation matrix of the finite differences method is substituted by a dense matrix, or, in other words, a nearest-neighbor interaction is exchanged for into a long-range one. The lack of sparsity of differentiation matrices is one of the drawbacks of spectral methods, especially when having to diagonalize large systems. However, they have the advantage of needing (a considerably) smaller number of grid points for getting the same accuracy as with finite difference methods.

For fixed (Dirichlet) boundary conditions, the Chebyshev spectral methods are the most suitable ones. There are several collocation schemes, the Gauss-Lobato being the most extensively used:

(S7) |

with being even or odd. The differentiation matrix is

(S8) |

The significant drawback of Chebyshev collocation is that the discretization matrix possesses a great number of spurious eigenvalues or outliers. They are approximately equal to . These outliers also have a significant non-zero real part, which increases when grows. This fact naturally reduces the efficiency of the method when performing numerical time-integration. However, it presents a higher spectral accuracy than the Fourier collocation method (see e.g. Supp_Book ()).

Several modifications must be introduced when applying spectral methods to polar coordinates. They basically rely on overcoming the difficulty of not having Dirichlet boundary conditions at and the singularity of the equations at that point. In addition, in the case of the Dirac equation, the spinor components can be either symmetric or anti-symmetric in their radial dependence, so the method described in Supp_Trefethen (); Herring () must be modified accordingly. As shown in the previously mentioned references, the radial derivative of a general function can be expressed as:

(S9) |

Notice that in this case, the collocation points must be taken as

(S10) |

but only the first points are taken so that the domain of the radial coordinate does not include . Analogously the differentiation matrix would possess now components, but only the upper half of the matrix, of size is used.

If the function that must be derived is symmetric or anti-symmetric, i.e. , with the upper (lower) sign corresponding to the (anti-)symmetric function, equation (S9) can be written as:

(S11) |

Thus, the differentiation matrix has a different form depending whether is symmetric or anti-symmetric:

(S12) |

with defined as in (S11).

### b.2 Existence of stationary solutions

Having established the basic features of spectral methods, we have enough tools for undertaking the numerical analysis of stationary states for the Nonlinear Dirac Equation. We recall that the dynamical equations in polar coordinates read as:

(S13) |

and stationary solutions can be written in the form with

(S14) |

where and could be chosen real-valued.

Thus, the equations for the stationary solutions read:

(S15) |

with .

There are no analytical solutions available for this system. For this reason, numerical methods must be used to this aim. Among them we have chosen to use fixed point methods, as the Newton-Raphson one Supp_numrec (), which requires the transformation of the set of two coupled ordinary differential equations (B.2) into a set of algebraic equations; this is performed by defining the set of collocation points , and transforming the derivatives into multiplication of the differentiation matrices and times the vectors and , respectively, being and as explained in the previous section. Thus, the discrete version of (B.2) reads:

(S16) |

It is important to notice that matrices and correspond to either and , depending on the symmetry of and , which, at the same time, depends on the value of the vorticity . If is even, then () is (anti-)symmetric and (). On the contrary, if is odd, then () is (anti-)symmetric and ().

In order to find the roots of the vector function , an analytical expression of the Jacobian matrix

(S17) |

must be introduced, with the derivatives expressed by means of spectral methods and the matrix must be evaluated at the corresponding grid points. The roots of , , are found by successive application of until convergence is attained. In our case, we have fixed as convergence condition that .

## Appendix C Linear stability of stationary solutions

To find the spectrum of linearization at a solitary wave in two spatial dimensions, we consider a solution in the form of a perturbed solitary wave solution:

(S18) |

with and corresponding to a small perturbation. The linearized equation on has the form

with a matrix-valued first order differential operator

(S19) |

where

with , , and with , evaluated at . To find the spectrum of the operator , we consider it in the space of -valued functions. The key observation which facilitates a computation of the spectrum is that the explicit form (S19) of contains , , , but not . As a consequence, is invariant in the spaces which correspond to the Fourier decomposition with respect to ,

the restriction of to each such subspace is given by

(S20) |

and this allows us to compute the spectrum of as the union of spectra of the one-dimensional spectral problems:

(S21) |

The spectrum is found by evaluating the functions appearing therein at the collocations points and substituting the partial derivatives by the corresponding differentiation matrix. At this point, one must be very cautious because, as also occurred with the Jacobian, there will be two different differentiation matrices in our problem. Thus, the representation of matrix will be:

In the 3d case, the NLD solitary waves were initially found in Supp_wakano-1966 () to be of the form with and with , satisfying

(S22) |

To study the linearization operator in the invariant space which has the same angular dependence as the solitary waves, we consider the perturbed solutions in the form

The linearized equation on is similar to (S19):

## Appendix D Critical frequency in 2D and 3D configurations

The figure below shows the critical frequencies for radially-symmetric collapse, as a function of the exponent associated with the nonlinearity in our generalized NLD model. For , the NLD solitary waves are linearly unstable. Below the linear instability disappears. For , with being the system dimension, there is no linear instability for .

## Appendix E Dynamics

Finally, we briefly discuss how the dynamics of Eq. (B.2) is simulated. In this case, Chebyshev spectral methods are not the most suitable ones, because of the presence of many outliers (see Supp_Book ()). In addition, as it was demonstrated in Supp_jpa (), finite difference methods also raise a number of concerns for the Dirac equation. Thus, the possibility that appears to us to be optimal presently is to use Fourier spectral methods, that works fairly well as long as the frequency is not close to zero.

Consequently, periodic boundary conditions must be supplied to our problem. This is less straightforward when working in polar coordinates in the domain . For this reason, we opt to work with a purely 2D problem in rectangular coordinates in the domain . The simulations of the paper have been performed with a Dormand-Prince numerical integrator using such a spectral collocation scheme with the aid of Fast Fourier Transforms (S4).

## References

- (56) W. Thirring, Ann. Phys. 3, 91 (1958).
- (57) J.P. Boyd, Chebyshev and Fourier spectral methods, 2nd. Edition. Dover (2001).
- (58) J. Cuevas-Maraver, P.G. Kevrekidis, A. Saxena, F. Cooper, and F. Mertens, in Ordinary and Partial Differential Equations (Nova Publishers, New York, 2015)
- (59) L.N. Trefethen, Spectral methods in MATLAB. SIAM, Philadelphia (2000).
- (60) G. Herring, L.D. Carr, R. Carretero-González, P.G. Kevrekidis, and D.J. Frantzeskakis. Phys. Rev. A 77, 023625 (2008).
- (61) M. Wakano, Prog. Theor. Phys. 35, 1117 (1966).
- (62) W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd. Edition. Cambridge University Press (2007).
- (63) J. Cuevas-Maraver, P.G. Kevrekidis, and A. Saxena, J. Phys. A: Math. Theor. 48, 055204 (2015).