Speed-of-light pulses in a massless nonlinear dirac equation

PHYSICAL REVIEW E 100, 022210 (2019) Speed-of-light pulses in a massless nonlinear Dirac equation Jesús Cuevas-Maraver ,1P. G. Kevrekidis,2Franz G. Mertens,3and Avadh Saxena4 1Grupo de Fisica No Lineal, Departamento de Fisica Aplicada I, Universidad de Sevilla, Escuela Politécnica Superior, C/Virgen de África, 7, 41011-Sevilla, Spain and Instituto de Matemáticas de la Universidad de Sevilla (IMUS), Edificio Celestino Mutis, Avenida Reina Mercedes s/n, 41012 Sevilla, Spain 2Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-4515, USA 3Physikalisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany 4Center for Nonlinear Studies and Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Received 8 May 2018; revised manuscript received 7 June 2019; published 13 August 2019) In this work, we explore a massless nonlinear Dirac equation, i.e., a nonlinear Weyl equation. We study the dynamics of its pulse solutions and find that a localized one-hump initial condition splits into a localized twohump pulse, while an associated phase structure emerges in suitable components of the spinor field. For times larger than a transient time tsthis pulse moves with the speed of light, effectively featuring linear wave dynamics and maintaining its shape (both in two and three dimensions). We show that for the considered nonlinearity, this pulse represents an exact solution of the nonlinear equation. Finally, we briefly comment on the generalization of the results to a broader class of nonlinearities. DOI: 10.1103/PhysRevE.100.022210 I. INTRODUCTION In the present work, we are interested in the massless analog of the Dirac equation. This setting, under the form of the so-called Weyl equation, has seen a surge of interest in recent years due to a variety of newly discovered materials called Weyl semimetals [1], e.g., NbAs and TaP [2,3]. These harbor chiral quasiparticles called Weyl fermions and possess topological surface states [4]. Weyl fermions exhibit linear dispersion, just like graphene, but are massless. In fact, Weyl semimetals are the three-dimensional (3D) analogs of graphene with broken spatial inversion or time reversal symmetry. In the Brillouin zone of such materials, linear dispersion arises around certain nodes, the so-called Weyl points, which always occur in pairs. In addition, in Weyl semimetals regions described by different Chern numbers are connected by unclosed lines, the so-called Fermi arcs [5], which can be experimentally observed using angle-resolved photoemission spectroscopy. The Fermi arc starts from one Weyl point and ends at the other one with opposite chirality. The Weyl points (or nodes) are essentially monopoles of the quantized Berry flux in the crystal momentum or reciprocal space. Photonic counterparts of Weyl semimetals have been observed in double-gyroid structures using angle-resolved microwave transmission measurements [6]. When time reversal and spatial inversion symmetries coexist in such a material, a pair of degenerate Weyl points may exist resulting in a Dirac semimetal [1,5], e.g., Cd3As2[7] and Na3Bi [8]. In parallel to these developments of chiefly linear Weyl physics, there has been an explosion of interest in the phenomenology of the nonlinear version of the Dirac equation and its solitary waves; a recent survey of the pertinent phenomenology can be found in [9]. While the relevant model in its massive Thirring form [10] was of interest to integrable systems and its Gross-Neveu/Soler form [11,12]led to extensive studies in solitary waves and their stability [9], arguably, part of the recent appeal of the model has been due to its applicability to a number of relevant physical setups. Among these, we note the dynamical evolution of Bose-Einstein condensates in the presence of honeycomb optical lattices [13–16], as well as the analogous propagation of light in honeycomb photorefractive lattices, the socalled photonic graphene [17–19]. These, in turn, motivated numerical and theoretical studies on the properties of these models and revealed crucial differences from their nonlinear Schrödinger cousins, including, e.g., the potential absence of the collapse instability for suitable parametric intervals in two-dimensional systems [20]. Some of the relevant applications of “massless Dirac physics” such as the optical ones of [6] constitute settings where the tuning of optical intensity may lead to the controllable introduction of nonlinearity. It is worthwhile to also note another piece of relevant motivation in the context of atomic Bose-Einstein condensates and their matter waves (but at the discrete rather than at the continuum level) in the work of [21]. It is the introduction of such a model blending the underlying linear massless 3D Dirac operator [although we also consider the two-dimensional (2D) analog thereof as well] and a cubic nonlinearity that we explore in the present setting, considering it as a nonlinear Weyl equation. Given the extensive number of corresponding studies at the Dirac level, but also its properties under Lorentz transformations and remarkable phenomenology reported below, we select the Gross-Neveu/Soler type of nonlinearity to formulate a massless nonlinear Dirac equation (massless NLDE, for short) that may be a starting point for exploring the interplay of nonlinearity with such linear operators in various contexts. The analytical (and numerical) results obtained herein can operate as a guide for examining other nonlinearities including the more relevant for atomic 2470-0045/2019/100(2)/022210(6) 022210-1 ©2019 American Physical Society JESÚS CUEVAS-MARAVER et al. PHYSICAL REVIEW E 100, 022210 (2019) condensates cases of a Kerr-type nonlinearity [21]. In that vein, at the end of the present exposition, we briefly touch upon the generalization of our findings to the case of the second nonlinearity. Our presentation and main results are as follows. First, we formulate the 3D massless NLDE equation and present some of its principal properties in three spatial dimensions, including most notably the observation that pulse-like initial data split into a two-humped ring density structure (acquiring a suitable phase in some of the spinor components). Beyond a transient time, the resulting density excitation is found to propagate at the speed of light. We demonstrate that, as a consequence, the resulting waveforms satisfy an effective 3D wave equation which is analytically solvable via suitable transformations. To corroborate these findings, we also examine the corresponding 2D case, and demonstrate the generic nature of the relevant phenomenology. It is important to point out here that our results bear fundamental differences from the recent one-dimensional (1D) corresponding study of [22]. In particular, here (a) there is a single radial density structure (as opposed to two pulses in the 1D case); (b) there appears a phase (vorticity) profile that we discuss below, and finally (c) the density decays with the distance, features that are particular to the higher-dimensional settings. Upon elucidating these traits, we summarize our main results, offer a number of remarks regarding other nonlinearities and applications, and propose a number of associated directions for future study. II. MASSLESS NLDE MODEL One of our principal motivations for utilizing the Gross- Neveu/Soler nonlinearity is that the associated quantity ψψ transforms as a scalar under the Lorentz transformation; nevertheless, we will also refer to the Kerr nonlinearity case in Sec. IV. The Gross-Neveu/Soler case suggests the corresponding massless Lagrangian density LNLDE =¯ ψ(iγμ∂μ)ψ+F(¯ ψψ),(1) where ψ(x,t)∈CN,x∈Rn, and γμ,0⩽μ⩽n,areN×N Dirac γ-matrices satisfying the anticommutation relations {γμ,γν}=2ημν, with ημν the Minkowski tensor [23], and ¯ ψ=ψ†γ0. We start by considering the 3D massless NLDE case in which the spinors have N=4 components. The relevant equation derived from the Lagrangian density of Eq. (1), in Cartesian coordinates assumes the form i∂tψ1=−i[(∂x−i∂y)ψ4+∂zψ3]−f(¯ ψψ)ψ1, i∂tψ2=−i[(∂x+i∂y)ψ3−∂zψ4]−f(¯ ψψ)ψ2,(2) i∂tψ3=−i[(∂x−i∂y)ψ2+∂zψ1]+f(¯ ψψ)ψ3, i∂tψ4=−i[(∂x+i∂y)ψ1−∂zψ2]+f(¯ ψψ)ψ4, with F(¯ ψψ)=f(¯ ψψ), which here is chosen as f(¯ ψψ)= g(|ψ1|2+|ψ2|2−|ψ3|2−|ψ4|2); we use the value of the prefactor g=1 herein. We note that Eq. (2)isaU(1), as well as translation-invariant, Hamiltonian system. Notice also that, contrary to the linear case, the transformation of the 4-spinor ψinto 2-spinors with left and right chirality does not decouple the equations. While here it can be thought that FIG. 1. Snapshots showing the evolution of an initial hump in the 3D massless nonlinear Dirac equation (NLDE). Top (middle) row shows an isosurface for 0.25 times the maximum of |ψ1( r,t)|2+ |ψ2( r,t)|2(|ψ3( r,t)|2+|ψ4( r,t)|2) at different values of time t;an octant of the sphere was removed to get a better visualization of the two-humped nature of the resulting structure. The bottom row shows the phase of ψ4( r,t) with emerging vorticity. In each picture, the axes cover the range [−40,40]. the nonlinearity plays the role of an “effective mass” in that sense (i.e., of avoiding decoupling), nevertheless recall that in the present speed-of-light limit, the mass is vanishing (hence this analogy is, in some sense, problematic). We employed the Wakano ansatz [24] ψ( r,0) =φ( r)=⎡ ⎢ ⎢ ⎢ ⎣ u(r) 0 iv(r) cos θ iv(r)sinθeiϕ ⎤ ⎥ ⎥ ⎥ ⎦ (3) to initialize the massless NLDE equation and the spherical frame will be useful in our analytical considerations below. Nevertheless, for our numerical solution of Eq. (2), we use a Fourier-series based pseudo-spectral method in Cartesian coordinates adapting the method used in [20] to the 4-spinor case in 3D (see Appendix A). The total mass stemming from the integration over space of the mass density ρ(r,t)=|ψ1(r,t)|2+|ψ2(r,t)|2+|ψ3(r,t)|2+|ψ4(r,t)|2 (4) is a conserved quantity of the model. We now integrate the massless NLDE model of Eq. (2)for typical pulse-like initial data of the form u(r)=1 2sech r 2,v(r)=0.(5) Notice that we found similar results for other forms of such initial data (e.g., Gaussian, etc.), and also by taking nonzero v(r). Figure 1shows snapshots of (isocontour density, as well as phase profiles of) both |ψ1( r,t)|2+|ψ2( r,t)|2and |ψ3( r,t)|2+|ψ4( r,t)|2, i.e., the densities of two component pairs. Note that the initial localized hump at the first spinor component transforms into a spherical shell two-hump structure that expands with time, whereas the initially null third and fourth spinor components also transform into a similar pattern, with the latter displaying vorticity, in line with the Wakano ansatz of Eq. (3); the second spinor component remains null (within machine precision) during the dynamical evolution. Figure 2shows the density at different times; notice 022210-2 SPEED-OF-LIGHT PULSES IN A MASSLESS NONLINEAR … PHYSICAL REVIEW E 100, 022210 (2019) FIG. 2. Density ρ(r,t) at different values of time in the 3D massless NLDE. the persistence of the two local maxima over the propagation time. The left panel of Fig. 3shows the position of the local density maximum; from this it is evident that the “ring” expands asymptotically with speed 1 beyond a transient time, i.e., for t>ts. The right panel of Fig. 3shows the quantity δ(t)=dn rf[ψ( r,t)],(6) which tends to zero for t>ts, with nbeing the number of spatial dimensions of the system. This quantity is a measure of the nonlinearity of the system during time evolution. It is thus clear that the nonlinear term becomes effectively “deactivated” for t>ts(in line also with the 1D massless Dirac case findings of [22]). As a result, the emerging two-humped pulses propagate at the “speed of light,” effectively satisfying the linear 3D wave equation, given the spontaneous vanishing of the nonlinear term f. In other words, given our observation that f→0, it is straightforward to show that each spinor U=ψ1,2,3,4satisfies the linear 3D wave equation of the form 1 c2∂2 t−∂2 x−∂2 y−∂2 zU=0,(7) FIG. 3. The left panel shows the position of the leftmost density local maximum r0for the 3D massless NLDE; the dashed red line corresponds to a slope 1 line to which r0(t) tends asymptotically. The right panel shows the evolution of δ(t)[seeEq.(6)], showcasing its asymptotic vanishing. Notice that we take r0=0 when there is a single local maximum. FIG. 4. Snapshots showing the evolution of an initial hump in the 2D massless NLDE. Top (middle) row shows the value of |ψ1( r,t)|2(|ψ2( r,t)|2) at different values of time t. Bottom row shows the development of a vortical phase structure within ψ2( r,t). In each picture, the axes cover the range [−80,80]. Recall that in our case c=1. Transformation w=rU can factor out the curvature term (2/r)∂rUand effectively restore a 1D wave equation in the radial variable, ultimately retrieving the full solution in the form U(r,t)=1 r[h(1)(r−ct)+h(2)(r+ct)].(8) In our simulations h(2) =0 and for t⩾tsthe four spinor components ψi=Ui=1 rh(1) i(x−ct) with four functions (i= 1,...,4) produce the two-hump structure seen in the density, per Eq. (4), in Fig. 2. Remarkably, this two-hump structure is a unique, previously undiscovered feature which differs qualitatively from the two-hump structure that was observed in the 1D massless NLDE [22]. In the latter the initial pulse splits symmetrically into two equal humps which move in opposite directions with the speed of light. III. 2D MASSLESS NLDE For the sake of comparison to the 3D case and a better understanding of the massless NLDE in general, next we study the 2D massless NLDE, which also showcases the generality of our findings. In two dimensions the γ-matrices are defined as γ0=σ3and γj=σ3σjwith j=1,2, where σ1,σ2, and σ3are the Pauli matrices. Explicitly, γ1=σ3σ1= iσ2and γ2=σ3σ2. In this context, the simplest case example of interest derived from the Lagrangian density in Eq. (1) can involve solely two spinor components according to the dynamical equations [20] i∂tψ1=−(i∂x+∂y)ψ2−f(¯ ψψ)ψ1,(9) i∂tψ2=−(i∂x−∂y)ψ1+f(¯ ψψ)ψ2, where ψ1,ψ2are the components of the spinor ψ∈C2and the nonlinearity is f=g(|ψ1|2−|ψ2|2). As in the 3D case, Eq. (9)isaU(1), as well as translation-invariant, Hamiltonian system. In Fig. 4, we once again explored the evolutionary dynamics of the 2D analog of the massless NLDE, initializing with a single humped waveform. The relevant results are also generic in their nature within the class of such initial data. We observe here too that a two-humped structure spontaneously emerges in a “ring” form (for the density), with the second 022210-3 JESÚS CUEVAS-MARAVER et al. PHYSICAL REVIEW E 100, 022210 (2019) component also featuring a phase profile, associated with the presence of vorticity in this spinor component. Similar to the 3D case, and showcasing the generality of our observations, we find that for t⩾ts≃20, the pulses propagate with constant speed, namely the speed of light, and the nonlinearity once again is made to vanish due to |ψ1|2=|ψ2|2, leading to an effectively linear dynamics. In a calculation similar to the above 3D case, given in detail in Appendix B(where both standing and traveling wave solutions of the 2D massless NLDE are explored, as applicable in the case of t⩾ts), we find that the effective dynamics for t⩾tsamounts to ∂tψ1+(∂x−i∂y)ψ2=0,(10) ∂tψ2+(∂x+i∂y)ψ1=0.(11) Combining the two equations (by taking, e.g., a timederivative of the first and substituting in the second), we obtain a 2D wave equation for both ψ1and ψ2, which, in turn, leads to the following expression for the density: ρ(r,t)=1 r|f(r−t)|2.(12) Constant factors are omitted here because ψ1and ψ2are solutions of effective linear equations. Notice the important 1/reffect, induced by the presence of the curvature also in the 2D system; such a term would be absent in a massless 1D Dirac setting [22]. IV. CONCLUSIONS, EXTENSIONS, AND FUTURE WORK Motivated by related studies in the photonic realm [6] (where nonlinearity can be naturally introduced) and in the matter wave realm [21] (where nonlinearity is effectively present due to interparticle interaction), we introduce a prototypical massless nonlinear Dirac equation in three dimensions, i.e., a nonlinear Weyl equation, and examine its analog in two dimensions. Considering the case invariant under Lorentz transformations, we utilize the Gross-Neveu/Soler nonlinearity in this study. We obtain pulse solutions of massless NLDE and their time evolution. Beyond a transient time t> ts, we find that these pulses move with the speed of light and satisfy an effectively linear (and explicitly solvable) wave equation. In the process, the role of curvature in the evolution of these pulses, as well as their two-humped structure and spontaneous phase development in suitable components, are also elucidated. Our results provide insight into the localization and dynamics of massless Dirac fields in the presence of nonlinearity. However, they also pose important questions that are especially relevant to address in future studies. In particular, from the theoretical standpoint, while the Gross-Neveu/Soler nonlinearity is of interest given its symmetry properties, in optical and atomic Dirac settings a nonlinearity involving solely |ψi|2ψiin the equation for the ith spinor component (i.e., a Kerr effect solely in each component from its own selfaction) is naturally of interest. It is then particularly relevant to separately explore the latter situation in both three and two dimensions. Remarkably our observations suggest that in these settings also, despite the variation of the nonlinearity, a similar phenomenology is observed. Namely, long-lived pulses appear to propagate outward at the speed of light for the class of initial data considered herein (given the similarity of the pertinent figures to the Soler case considered herein, we do not show them here). Moreover, it would be particularly interesting in experimentally realized photonic crystals, or perhaps in theoretically proposed atomic settings, to explore the possibility of observing this intriguing interplay of linear phenomena (some of which are discussed above) and nonlinearity. ACKNOWLEDGMENTS J.C.-M. thanks financial support from MAT2016-79866- R project (AEI/FEDER, UE). P.G.K. acknowledges support from NSF-PHY-1602994. F.G.M. acknowledges the hospitality of the Center for Nonlinear Studies and Theoretical Division at LANL. This work was supported, in part, by the US Department of Energy. We acknowledge the useful comments of Andrew Comech. APPENDIX A: NUMERICAL METHODS We briefly describe in this Appendix the numerical methods employed for integrating equations (2). For a more complete description of such methods, please refer to [9]. The first step to follow is to implement a grid and a method for discretizing the spatial derivatives of the partial differential equations (PDEs) of Eq. (2) and, consequently, transforming them into a set of coupled ordinary differential equations. Finite difference methods do not usually work well in nonlinear Dirac equations. Instead, one must make use of spectral methods. In our case, as we are dealing with numerical integrations of PDEs and the pulse tends to infinity in an exponential way, a well-suited choice is the Fourier-series based pseudo-spectral method. Such a method requires the use of periodic boundary conditions and an equispaced grid. The implementation is quite simple, as it basically consists of performing direct and inverse Fourier transforms. That is, if we denote by U(x,y,z) any of the spinor components ψ1,2,3,4(x,y,z) then the derivative with respect to, e.g., direction xis given by ∂xU(x,y,z)=F−1 x{ikFx[U(x,y,z)]}, where Fxand F−1 xdenote the direct and inverse, respectively, 1D Fourier transform in the direction x. Notice that U(x,y,z) actually represents an M×M×Marray which only takes values at the grid points. Because of this, Fourier transforms can be accomplished by means of the fast Fourier transform (FFT) and k≡{km}is a vector with Mcomponents given by km=mπ/Lfor m<Mand kM=0. Once we defined our set of ordinary differential equations, the second ingredient is the numerical integrator. In our case, the integrator we prefer to use for simplicity and accuracy is the Dormand-Prince [25] algorithm. All the above schemes were implemented using MATLAB in a desktop PC with 8 Gb of RAM. In our particular case, the domain was a box of size (−L,L]×(−L,L]×(−L,L] with L=48 and the lattice discretization parameter was h= 0.8; with these data, the number of grid points per dimension 022210-4 SPEED-OF-LIGHT PULSES IN A MASSLESS NONLINEAR … PHYSICAL REVIEW E 100, 022210 (2019) was M=120 so that the total number of ordinary differential equations to integrate was equal to 4M3∼7×106. Reaching the value of time t=30 takes around a day. We finally remark that one of the advantages of using spectral methods is that it works even with a rough discretization and that the value of the wave function at different values of (x,y,z) out of the grid can be accurately attained by using spline interpolation. APPENDIX B: TWO-DIMENSIONAL MASSLESS NONLINEAR DIRAC CASE For the 2D case, we can simplify the relevant analysis by using the polar coordinates, where the equations take the form i∂tψ1=−e−iθi∂r+∂θ rψ2−f(ψ1,ψ 2)ψ1, i∂tψ2=−eiθi∂r−∂θ rψ1+f(ψ1,ψ 2)ψ2.(B1) We performed simulations of the relevant 2D analog of the massless NLDE model with the following initial conditions, in line with those used in [20] (for the massive case), ψ( r,0) =φ( r)=u(r)eiSθ iv(r)ei(S+1)θ,(B2) with Sbeing the vorticity and, just as in the 3D case, u(r)=1 2sech r 2,v(r)=0.(B3) The two-humped ring nature of the resulting dynamics (as well as the spontaneous emergence of vorticity) is revealed in the figure shown in the main text. Here, for completeness, we show in Fig. 5the density ρ(r,t) at different times. Also, the left panel of Fig. 6shows the position of the maximum of the density of the ring; from this figure, it is evident that the ring expands asymptotically with the speed of light (as in three dimensions). The right panel of the figure shows δ(t) defined in the main text. The vanishing of this quantity once again indicates the spontaneous “self-annihilation” of the nonlinear terms. Thus, for t>tsthis structure remains the same for all FIG. 5. Density ρ(r,t) at different values of time in the 2D massless nonlinear Dirac equation. FIG. 6. Left panel shows the position of the density maximum r0for the massless 2D nonlinear Dirac equation; dashed red line corresponds to a slope 1 line to which r0(t) tends asymptotically. Right panel shows the evolution of δ(t) and its effective vanishing beyond a transient time. times, but its density is reduced by an r-dependent factor and it moves with the speed of light. Therefore, this structure is supposed to be a solution of the linear 2D form of the equation. In the following we will show that it is indeed the case by solving exactly the associated linear PDEs. We use the following Ansätze: ψ1=e−iθ 2˜ ψ1(r,t),ψ 2=eiθ 2˜ ψ2(r,t),(B4) and obtain the following relations: ∂t˜ ψ1=− ∂r+1 2r˜ ψ2,∂ t˜ ψ2=− ∂r+1 2r˜ ψ1.(B5) First we consider time-independent solutions: ˜ ψ1(r,t)=u(r),˜ ψ2(r,t)=iv(r),(B6) and obtain v+1 2rv=0,u+1 2ru=0.(B7) We find that v=c2e−1 2rdr =c2e−1 2ln r=c2 √r,(B8) with arbitrary c2. In a similar way we have u=c1/√rwith c1 being an arbitrary constant. Thus, =(ψ1,ψ 2)T=1 √rc1e−iθ 2,ic2eiθ 2T.(B9) Denoting ˜u=˜ ψ1+˜ ψ2and ˜v=˜ ψ1−˜ ψ2we decouple Eq. (B5)as ∂t˜u=− ∂r+1 2r˜u,∂ t˜v=+ ∂r+1 2r˜v.(B10) Next, we consider stationary solutions by employing the ansatz ˜u(r,t)=e−iωtu(r),˜v(r,t)=e−iωtv(r),(B11) and obtain u(r)=c1 eiωr √r,v(r)=c2 e−iωr √r.(B12) Finally, we get the solutions ψ1=1 2e−iωte−iθ/2[u(r)+v(r)],(B13) ψ2=1 2e−iωte+iθ/2[u(r)−v(r)]. 022210-5 JESÚS CUEVAS-MARAVER et al. PHYSICAL REVIEW E 100, 022210 (2019) Next we consider traveling wave solutions. We take the decoupled Eq. (B10) and make the ansatz ˜u=1 rβh(1)(r−ct),(B14) with a solution c=1, β=1/2, where h(1)(r−ct)isan arbitrary function, and ˜v=1 rβh(2)(r−ct),(B15) with a solution c=−1, β=1/2. That is, we have two solutions [with h(1) and h(2) arbitrary]: ˜uI=1 √rh(1)(r−ct),c=1,˜vI=0, ˜uII =0,˜vII =1 √rh(2)(r−ct),c=−1.(B16) Thus, the solutions of the original equations [with h(1) and h(2) arbitrary] are ψ1,I=1 2e−iθ/21 √rh(1)(r−t), ψ2,I=1 2e+iθ/21 √rh(1)(r−t). (B17) ψ1,II =1 2e−iθ/21 √rh(2)(r+t), ψ2,II =1 2e+iθ/21 √rh(2)(r+t). [1] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90,015001 (2018). [2] S.-Y. Xu et al.,Science 349,613 (2015). [3]B.Q.Lvet al.,Phys.Rev.X5,031013 (2015). [4]A.A.BurkovandL.Balents,Phys. Rev. Lett. 107,127205 (2011). [5] O. Vafek and A. Vishwanath, Ann. Rev. Cond. Mat. Phys. 5,83 (2014). [6] L. Lu, Z. Wang, D. Ye, L. Ran, J. D. Joannopoulos, and M. Soljacic, Science 349,622 (2015). [7] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88,125427 (2013). [8] Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys.Rev.B85,195320 (2012). [9]J.Cuevas-Maraver,N.Boussaid,A.Comech,R.Lan,P.G. Kevrekidis, and A. Saxena, Solitary Waves in the Nonlinear Dirac Equation, Springer Series in Understanding Complex Systems, edited by V. Carmona et al., Nonlinear Systems Vol. 1, (Springer, New York, 2018), p. 89; arXiv:1707.01946. [10] W. E. Thirring, Ann. Phys. (NY) 3,91 (1958). [11] D. J. Gross and A. Neveu, Phys. Rev. D 10,3235 (1974). [12] M. Soler, Phys. Rev. D 1,2766 (1970). [13] L. H. Haddad and L. D. Carr, Physica D 238,1413 (2009). [14] L. H. Haddad and L. D. Carr, Europhys. Lett. 94,56002 (2011). [15] L. H. Haddad and L. D. Carr, New J. Phys. 17,113011 (2015). [16] L. H. Haddad, K. M. O’Hara, and L. D. Carr, Phys. Rev. A 91, 043609 (2015). [17] O. Peleg, G. Bartal, B. Freedman, O. Manela, M. Segev, and D. N. Christodoulides, Phys. Rev. Lett. 98,103901 (2007). [18] M. J. Ablowitz, S. D. Nixon, and Y. Zhu, Phys. Rev. A 79, 053830 (2009). [19] M. J. Ablowitz and Y. Zhu, Phys. Rev. A 82,013840 (2010). [20] J. Cuevas-Maraver, P. G. Kevrekidis, A. Saxena, A. Comech, and R. Lan, Phys. Rev. Lett. 116,214101 (2016). [21] C. Shang, Y. Zheng, and B. A. Malomed, Phys. Rev. A 97, 043602 (2018). [22] N. R. Quintero, F. G. Mertens, F. Cooper, A. Saxena, and A. R. Bishop, Phys. Rev. E 96,052219 (2017). [23] B. de Wit and J. Smith, Field Theory in Particle Physics (North Holland, New York, 1986). [24] M. Wakano, Progr. Theor. Phys. 35,1117 (1966). [25] E. Hairer, P. S. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems (Springer, Berlin, 2008). 022210-6