Nonlinear Dirac equation solitary waves in external fields

PHYSICAL REVIEW E 86, 046602 (2012) Nonlinear Dirac equation solitary waves in external fields Franz G. Mertens,1,*Niurka R. Quintero,2,†Fred Cooper,3,4,‡Avinash Khare,5,§and Avadh Saxena4,| 1Physikalisches Institut, Universit¨ at Bayreuth, D-95440 Bayreuth, Germany 2IMUS and Departamento de Fisica Aplicada I, E.P.S. Universidad de Sevilla, 41011 Sevilla, Spain 3Santa Fe Institute, Santa Fe, New Mexico 87501, USA 4Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 5Indian Institute of Science Education and Research, Pune 411021, India (Received 13 August 2012; published 24 October 2012) We consider nonlinear Dirac equations (NLDE’s) in the 1 +1 dimension with scalar-scalar self-interaction g2 κ+1(¯ )κ+1in the presence of various external electromagnetic fields. We find exact solutions for special external fields and we study the behavior of solitary-wave solutions to the NLDE in the presence of a wide variety of fields in a variational approximation depending on collective coordinates which allows the position, width, and phase of these waves to vary in time. We find that in this approximation the position q(t) of the center of the solitary wave obeys the usual behavior of a relativistic point particle in an external field. For time-independent external fields, we find that the energy of the solitary wave is conserved but not the momentum, which becomes a function of time. We postulate that, similarly to the nonlinear Schr¨ odinger equation (NLSE), a sufficient dynamical condition for instability to arise is that dP(t)/d ˙ q(t)<0. Here P(t) is the momentum of the solitary wave, and ˙ qis the velocity of the center of the wave in the collective coordinate approximation. We found for our choices of external potentials that we always have dP(t)/d ˙ q(t)>0, so, when instabilities do occur, they are due to a different source. We investigate the accuracy of our variational approximation using numerical simulations of the NLDE and find that, when the forcing term is small and we are in a regime where the solitary wave is stable, that the behavior of the solutions of the collective coordinate equations agrees very well with the numerical simulations. We found that the time evolution of the collective coordinates of the solitary wave in our numerical simulations, namely the position of the average charge density and the momentum of the solitary wave, provide good indicators for when the solitary wave first becomes unstable. When these variables stop being smooth functions of time (t), then the solitary wave starts to distort in shape. DOI: 10.1103/PhysRevE.86.046602 PACS number(s): 05.45.Yv, 03.70.+k, 11.10.Lm I. INTRODUCTION Classical solutions of nonlinear field equations have a long history as a model of extended particles [1–3]. In 1970, Soler [3] proposed that the self-interacting 4-Fermi theory was an interesting model for extended fermions. Later, Strauss and Vasquez [4] were able to study the stability of this model under dilatation and found the domain of stability for the Soler solutions. Solitary waves in the 1 +1 dimensional nonlinear Dirac equation (NLDE) have been studied [5,6]in the past in the case of massive Gross-Neveu [7] (with N=1, i.e., just one localized fermion) and massive Thirring [8] models. In those studies, it was found that these equations have solitary-wave solutions for both scalar-scalar (S-S) and vector-vector (V-V) interactions. The interaction between solitary waves of different initial charge was studied in detail for the S-S case in the work of Alvarez and Carreras [9] by Lorentz boosting the static solutions and allowing them to scatter. Recently, we extended the solutions previously found to a more general interaction of the form g2 κ+1(¯ )κ+1 [10]. For the nonrelativistic limit of the NLDE, namely the nonlinear Schr¨ odinger equation (NLSE), there have been *[email protected] †[email protected] ‡[email protected] §[email protected] |a[email protected]v recent studies of the behavior of the forced NLSE. Using a collective coordinate (CC) theory, the authors found [11–14] that a sufficient dynamical condition for instability to arise is that dp(t)/dv < 0.Here p(t) is the normalized canonical momentum p(t)=1 M(t) ∂L ∂˙ q,M(t)=dx(x,t)(x,t)isthe mass, and ˙ q(t)=v(t) is the velocity of the solitary wave. One of the points we will investigate in the paper is whether this dynamical stability criterion is also valid for the NLDE. There has been recent interest in the stability of NLDE with higher-order nonlinearity [15]. Comech (private communication) has been able to prove that for κ=1, the Vakhitov-Kolokolov [16] criterion guarantees linear stability in the nonrelativistic regime of the NLDE equation for solutions of the form (in the rest frame) (x,t)=ψ(x)e−iωt, where ωis less than but approximately equal to the mass parameter min the Dirac equation. He was also able to show linear instability in the same nonrelativistic regime for κ⩾3. This is the first rigorous result for the Dirac equation, but it only applies in the nonrelativistic regime. Here we want to understand if we can determine in the relativistic regime for what values of ωdo the solitary waves become unstable, with and without forcing terms, even when they are stable in the nonrelativistic regime. What we find is that when the solitary waves are only metastable for the unforced problem, the critical time for the solitary wave to become unstable in the forced problem for weak forcing is similar to the critical time in the unforced problem. When the solitary wave maintains its basic shape, the CC equations give a good description of the actual time evolution at all times. This is true for weak 046602-1 1539-3755/2012/86(4)/046602(20) ©2012 American Physical Society MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) ramp potentials, harmonic potentials, and spatially periodic potentials, when ω>ω c, and ωcis the critical value above which the unforced solitary wave is stable. The collective coordinates q(t) and P(t), the position and momentum of the solitary wave, are “smooth” functions of tfor the external potentials we have chosen. Their counterparts in the numerical simulation are the first moment of the charge density and the total momentum of the numerical solution. When the numerical evolution of these counterparts to the collective coordinates start deviating from their CC values, this is a signal that the shape of the solitary wave is beginning to change. This usually rapidly develops into non-smooth behavior of q(t) and P(t) in the numerical solution. This is how we determine the onset of the instability time tcfor the forced NLDE solitary wave. Unfortunately, for the potentials we study, we always obtain dp/d ˙ q>0, which fulfills a necessary condition for stability. Thus, this criterion does not yield a prediction of the instabilities. This paper is organized as follows: In Sec. II,wereview the known exact solutions for the unforced NLDE and discuss the conservation laws that govern their behavior. In Sec. III we extend Bogolubsky’s discussion [17] of the stability of these solitary-wave solutions to changes in the frequency ω for arbitrary nonlinearity parameter κ. In Sec. IV we consider the NLDE in external electromagnetic fields and find particular exact solitary-wave solutions and discuss the stability of these solutions. In Sec. Vwe introduce our variational method based on using for our variational wave functions the exact wave functions for the solitary waves of the unforced problem, with the position, width parameter, and phase of these solutions being promoted to collective coordinates depending on time. We write the relativistic equations for these collective coordinates which are similar to point particle relativistic dynamical equations. The potential the average position of the solitary wave sees is a particular average of the external potential weighted with the charge density. In Sec. VI we postulate our stability criterion for an arbitrary external potential based on just solving the CC equations. This condition is a sufficient condition for instability. In Sec. VII we examine and solve the collective coordinate (CC) equations for three types of potentials—a ramp potential, a harmonic potential, and a spatially periodic potential. We also compare the solution to the CC equations to the numerical simulation of the NLDE equation. We state our conclusions in Sec. VIII. In the appendix we discuss identities that are obeyed by the solutions in the rest frame. II. REVIEW OF EXACT SOLUTIONS TO THE NLDE In this section, we review the exact solutions to the NLDE, using the notation of Ref. [10]. We are interested in solitarywave solution of the NLDE given by (iγμ∂μ−m)+g2(¯ )κ=0.(2.1) These equations can be derived in a standard fashion from the Lagrangian density L=i 2[¯ γμ∂μ−∂μ¯ γμ]−m¯  +g2 κ+1(¯ )κ+1. (2.2) For solitary-wave solutions, the field goes to zero at infinity. It is sufficient to go into the rest frame, since the theory is Lorentz invariant and the moving solution can be obtained by a Lorentz boost. In the rest frame we consider solutions of the form (x,t)=e−iωtψ(x).(2.3) We are interested in bound-state solutions that correspond to positive frequency ω⩾0 and which have energies in the rest frame less than the mass parameter m, i.e., ω<m. In our previous paper [10], we chose the representation γ0=σ3,iγ1=σ1. Here, to make contact with the numerical simulations paper of Alvarez and Carreras [9] we choose instead γ0=σ3;γ1=iσ2.Defining A,B via ψ(x)=A(x) iB(x)=R(x)cos θ isin θ,(2.4) we obtain the following equations for Aand B: dA dx +(m+ω)B−g2(A2−B2)κB=0, (2.5) dB dx +(m−ω)A−g2(A2−B2)κA=0. A first integral of these equations can be obtained by realizing that from energy-momentum conservation we have ∂μTμν =0; Tμν =i 2[¯ γμ∂ν−∂ν¯ γμ]−gμν L.(2.6) Thus, for stationary solutions, T10 =const,T 11 =const.(2.7) Now, using (2.3), we obtain T11 =ωψ†ψ−m¯ ψψ +LI;LI=g2 κ+1(¯ ψψ)κ+1.(2.8) For solitary-wave solutions vanishing at infinity, the constant is zero and we get the useful first integral, T11 =ωψ†ψ−m¯ ψψ +LI=0.(2.9) Multiplying the equation of motion on the left by ¯ and using (2.3) we have that (κ+1)LI=−ωψ†ψ+m¯ ψψ −¯ ψiγ1∂1ψ. (2.10) Therefore, we can rewrite T11 =0as ωκψ†ψ−mκ ¯ ψψ −¯ ψiγ1∂1ψ=0.(2.11) For the Hamiltonian density we obtain H=T00 =i 2[¯ γ1∂x−∂x¯ γ1]+m¯  −LI ≡h1+h2−h3.(2.12) Each of the hiare positive definite. From Eqs. (2.9) and (2.10), one has the relationship κLI=−¯ ψiγ1∂xψ. (2.13) From this, we have h3=1 κh1(2.14) and, in particular for κ=1, H=m¯ ψψ. 046602-2 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) In terms of R,θ one has ¯ ψiγ1∂1ψ=ψ†ψdθ dx.(2.15) This leads to the simple differential equation for θfor solitary waves, dθ dx =−ωκ+mκcos 2θ;ωκ≡κω;mκ=κm. (2.16) The solution is (in this section and what follows we will choose the position of the solitary wave to be initially at x0=0) θ(x)=tan−1(αtanh βκx),(2.17) where α=mκ−ωκ mκ+ωκ1/2 =m−ω m+ω1/2 ,β κ=m2 κ−ω2 κ1/2. (2.18) Thus, we have tan θ(x)=αtanh βκx, sin2θ(x)=(m−ω)sinh 2βκx mcosh 2βκx+ω; cos2θ(x)=(m+ω) cosh2βκx mcosh 2βκx+ω,(2.19) where we have used the identities 1+α2tanh2βkx=mcosh 2βkx+ω m+ωsech2βkx, (2.20) 1−α2tanh2βkx=ωcosh 2βkx+m m+ωsech2βkx. Solving Eq. (2.9) for R2we obtain R2=(κ+1)(mcos 2θ−ω) g2(cos 2θ)κ+11/κ .(2.21) Now we have dθ dx =β2 κ ωκ+mκcosh 2βκx=−ωκ+mκcos 2θ, (2.22) where βκ=m2 κ−ω2 κ=κ√m2−ω2,so cos 2θ=mκ+ωκcosh 2βκx ωκ+mκcosh 2βκx=m+ωcosh 2βκx ω+mcosh 2βκx.(2.23) We can rewrite R2using the right-hand side of Eq. (2.22) as R2=ω+mcosh 2βκx m+ωcosh 2βκx (κ+1)β2 κ g2κ2(m+ωcosh 2βκx)1/κ . (2.24) Using the identities in Eq. (2.20), we obtain the alternative expression R2=1+α2tanh2βκx 1−α2tanh2βκx ×sech2βκx(κ+1)β2 κ g2κ2(m+ω)(1 −α2tanh2βκx)1/κ .(2.25) In particular, for κ=1, R2=2(m−ω) g2 (1 +α2tanh2βx) (1 −α2tanh2βx)2sech2βx (2.26) and A2=R2cos2θ=2 g2 (m2−ω2)(m+ω) cosh2βx (m+ωcosh 2βx)2, (2.27) B2=R2sin2θ=2 g2 (m2−ω2)(m−ω)sinh 2βx (m+ωcosh 2βx)2. For arbitrary κwe have A=(m+ω) cosh2(κβx) m+ωcosh(2κβx)(κ+1)β2 g2(m+ωcosh(2κβx))1 2κ , B=(m−ω)sinh 2(κβx) m+ωcosh(2κβx)(κ+1)β2 g2(m+ωcosh(2κβx))1 2κ . (2.28) Because of Lorentz invariance we can find the solution in a frame moving with velocity vwith respect to the rest frame. The Lorentz boost is given in terms of the rapidity variable η as follows (here c=1): v=tanh η;γ=1 √1−v2=cosh η;sinhη=v √1−v2. (2.29) In the moving frame, the transformation law for spinors tells us that (x,t)=cosh(η/2) sinh(η/2) sinh(η/2) cosh(η/2)  ×0 1[γ(x−vt),γ (t−vx)] 0 2[γ(x−vt),γ (t−vx)],(2.30) since cosh(η/2) =(1 +γ)/2; sinh(η/2) =(γ−1)/2.(2.31) This in component form, 1(x,t)=[cosh(η/2)A(x)+isinh(η/2)B(x)]e−iωt, (2.32) 2(x,t)=[sinh(η/2)A(x)+icosh(η/2)B(x)]e−iωt, where x=γ(x−vt); t=γ(t−vx).(2.33) Note that cosh2(η/2) +sinh2(η/2) =cosh η=γ. A. Conservation laws of the NLDE The Lagrangian is invariant under the transformation of phase →ei, which by Noether’s theorem leads to the conserved current, ∂μjμ(x)=0; jμ=¯ γμ. (2.34) This leads to charge conservation, Q=dx†, (2.35) which, for the solitary-wave solution, leads to Q=dx(A2+B2)=1 κβ (κ+1)β2 g2(m+ω)1/κ Iκ(α2),(2.36) 046602-3 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) where Iκ(α2)=1 −1 dy 1+α2y2 (1 −y2)(κ−1)/κ [1 −α2y2](κ+1)/κ =B1 2,1 κ2F11+1 κ,1 2,1 2+1 κ;α2 +α2B3 2,1 κ2F11+1 κ,3 2,3 2+1 κ;α2, (2.37) and 2F1is a hypergeometric function. Here B(x,k) denotes the βfunction. We also have energy-momentum conservation Eq. (2.6) leading to conservation of energy and momentum, E=T00dx;P=T01dx. (2.38) Because of Lorentz invariance it is sufficient to calculate the energy-momentum tensor in the comoving frame v=0. The energy momentum tensor in an arbitrary frame is then given by Tμν =μ αν βTαβ ;μ α=cosh ηsinh η sinh ηcosh η.(2.39) In the rest frame of the solitary wave, for the unperturbed system, one has that T00 =h11−1 κ+h2,(2.40) where h1=R2(x)dθ dx =κβ2 m+ωcosh(2κβx) ×(κ+1)β2 g2(m+ωcosh(2κβx))1/κ ,(2.41) h2=m¯ ψψ =m(A2−B2) =m(κ+1)β2 g2(m+ωcosh(2κβx))1/κ .(2.42) Integrating in the rest frame, we get, for the rest-frame energy, E0=H11−1 κ+H2,(2.43) where H1=dxh1=β m+ω(κ+1)β2 g2(m+ω)1/κ ×B1 2,1+1 κ2F11+1 κ,1 2,3 2+1 κ;α2,(2.44) H2=dxh2=1 κβ (κ+1)β2 g2(m+ω)1/κ ×B1 2,1 κ2F11 κ,1 2,1 2+1 κ;α2.(2.45) Since in the rest frame for stationary solutions T11 =T01 = 0, the energy of the solitary wave in the moving frame is just E=E0cosh η=γE 0;P=E0sinh η, (2.46) so the norm E2−P2=E2 0=M2 0. In particular, for κ=1and m=1, we have that M0=2 g2Qsinh−1g2Q 2;Q=2√1−ω2 g2ω.(2.47) We also have H1=− 2√1−ω2−2 tanh−11−ω 1+ω g2, (2.48) H2= 4 tanh−11−ω 1+ω g2. III. STABILITY OF EXACT SOLUTIONS A. Stability to changes in the frequency at fixed charge Bogolubsky [17] suggested that the stability could be ascertained by looking at variations of the wave function, keeping the charge fixed and seeing if the solution was a minimum (stable to that variation) or maximum (unstable to that variation) of the Hamiltonian as a function of the parameter ω. This principle has been very useful in the past in determining the stability of scalar wave equations that are Hamiltonian dynamical systems. If the variation decreased the energy, it turned out that the solitary waves were unstable. Since in higher dimensions there are many degrees of freedom for perturbing the system, this criterion is a sufficient condition for instability. For the Dirac case, we have found from our numerical simulations that this criterion does not determine the critical ωexcept when κ=1[18], the case originally studied by Bogolubsky [17]. Assuming we know the wave function at the value of ωcorresponding to a fixed charge Q, if we change the parametric dependence on ω, this also changes the charge. This can be corrected by assuming that the new wave function has a new normalization that corrects for this. That is, if we parametrize a rest-frame solitary-wave solution of the NLDE which has a charge Q[ω]by ψs(x,t)=χs(x,ω)e−iωt,(3.1) then we choose our slightly changed wave function to be ˜ ψ[x,t,ω,ω]=√Q[ω] √Q[ω]χs(x,ω)e−iωt ≡f(ω,ω)χs(x,ω)e−iωt.(3.2) Thewavefunction ˜ ψ[x,t,ω,ω] then has the same charge as ψ[x,t,ω]. Inserting this wave function into the Hamiltonian, we get a new probe Hamiltonian Hpdepending on both ω,ω. As a function of ωthis new Hamiltonian is stationary at the value ω=ω. The criterion Bogolubsky proposed [17]is that the solitary wave is stable (unstable) with respect to this variation in ωaccording to whether this new Hamiltonian has a minimum (maximum) at ω=ω. What we will find for κ=1 is that there is a critical value of ω(determined by the coupling g and Q) below which the solitary wave is unstable, and this result is borne out by numerical simulations 046602-4 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) which we will present below. However, we will present in another paper numerical simulations at arbitrary κwhich suggest that this method provides only a sufficient condition for instability [18]. The probe Hamiltonian has the form Hp[ω,ω]=H1[ω]f(ω,ω)2−1 κf(ω,ω)2(κ+1) +H2[ω]f(ω,ω)2.(3.3) For κ=1wehavethatf(ω,ω)2=β[ω]ω β[ω]ω,where β[ω]=√1−ω2. We then find that the first derivative of Hpwith respect to ωevaluated at ω=ωis indeed zero. The second derivative evaluated at ω=ωleads to the following expression: ∂2Hp ∂ω2ω=ω=− 2√1−ω2(ω2−3) +4 tanh−11−ω 1+ω g2ω2(ω2−1)2. (3.4) This function is zero at ωc=0.697586 and the second derivative is negative below this value of ωshowing an instability. In our numerical simulations of the unforced NLDE [18], we find that, below this value, the solitary waves are metastable, with the time for the instability to set in increasing exponentially as a function of ωfor ω<ω c. IV. NLDE IN EXTERNAL ELECTROMAGNETIC FIELDS We add electromagnetic interactions through the gauge covariant derivative i∂μ→(i∂μ−eAμ), (4.1) and then, under the combined transformations →ei(x);Aμ→Aμ−1 e∂μ;¯ →¯ e−i(x),(4.2) the Lagrangian is invariant. Again, the conserved current is given by Eq. (2.34). The gauge-invariant Lagrangian for the external field problem is L=i 2[¯ γμ∂μ−∂μ¯ γμ]−m¯  +g2 κ+1(¯ )κ+1−e¯ γμAμ. (4.3) One again finds that energy-momentum is conserved, with E and Pgiven by Eq. (2.38). Because of Lorentz invariance, it is sufficient to calculate the energy-momentum tensor in the comoving frame v=0. The energy-momentum tensor in an arbitrary frame is then given by Eq. (2.39). Using the freedom of gauge invariance, one can choose the axial gauge A1=0, eA0=V(x). In the axial gauge the Dirac equation becomes iγμ∂μ−m +g2(¯ )κ−γ0V(x)=0.(4.4) Going into the rest frame and choosing for ψ0the representation of Eq. (2.4), we find that the Dirac equation becomes ∂xA+(m+ω)B−g2[A2−B2]κB−V(x)B=0, (4.5) ∂xB+(m−ω)A−g2[A2−B2]κA+V(x)A=0. As in the case of the standard NLDE equation, we can again look for static solutions of the form ψ0, where A=Rcos θ and B=Rsin θ. The conservation of T11 for static solutions that vanish at infinity yields the equation T11 =ωψ†ψ−m¯ ψψ +LI=0,(4.6) where LI=g2 κ+1(¯ ψψ)κ+1−V(x)ψ†ψ. (4.7) Multiplying the Dirac equation [Eq. (4.4)] on the left by ¯  and using Eq. (2.3), we obtain ωψ†ψ+i¯ ψγ1∂xψ−m¯ ψψ +g2(¯ ψψ)κ+1−V(x)ψ†ψ=0. (4.8) The energy density is given by h=i 2[¯ γ1∂x−∂x¯ γ1]+m¯  −g2 κ+1(¯ )κ+1+V(x)† ≡h1+h2−h3+h4.(4.9) From Eq. (4.6),(4.7), and Eq. (4.8) we again find that h3=1 κh1.(4.10) Multiplying the Dirac equation [Eq. (4.4)] on the left by ¯  and using Eq.(2.3), eliminating the self-interaction term using Eq. (4.6), we obtain i¯ ψγ1∂xψ+κ[m¯ ψψ +V(x)ψ†ψ−ωψ†ψ]=0.(4.11) Realizing that the equations for the standard NLDE are modified by replacing ωby ω−V(x), we obtain, for θ(x), dθ dx =−κω +κmcos(2θ)+κV(x),(4.12) and for R(x) we obtain R2=(κ+1)(mcos 2θ−ω+V(x)) g2(cos 2θ)κ+11/κ .(4.13) Notice, if we now choose V(x)=μcos 2θ,wearriveatthe solutions found earlier with m→m+μ. Note that for a bound state without the external potential we needed ω<m.Sonow we need ω<m+μfor βkto be real. The potential is then V1(x)=μ(m+μ)+ωcosh 2βκx ω+(m+μ) cosh 2βκx; (4.14) βκ=κ(m+μ)2−ω2. The eigenvalue ωis set by fixing the charge Qof the solitary wave. We can choose a negative μsatisfying ω<m+μ so the potential at small xlooks much like a harmonic trap. Namely, choosing κ=1, μ=−1/4, m=1, ω=1/2, then V1(x)=−2 cosh √5x 2+3 12 cosh √5x 2+8 .(4.15) This is plotted in Fig. 1. The charge density for the solitary wave corresponding to this external potential is plotted in Fig. 2. 046602-5 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) 2 1 1 2 x 0.24 0.23 0.22 0.21 0.20 V x FIG. 1. (Color online) Potential V1(x) corresponding to Eq. (4.14) with μ=−1/4, ω=1/2. At small x,V(x) has the expansion V1(x)≃−1 4+x2 32 −13x4 1536 +O(x6).(4.16) Another solution is found if we let V2(x)=μsin(2θ(x)).(4.17) 4 2 2 4 x 0.1 0.2 0.3 0.4 0.5 x FIG. 2. (Color online) Charge density R2(x) corresponding to the potential V1(x)givenbyEq.(4.14) with μ=−1/4, ω=1/2. The solution to the differential equation, θ(x)=−κω +κmcos(2θ(x)) +κμsin(2θ(x)), (4.18) then is found to be θ(x)=tan−1μ+μ2+m2−ω2tanh [κ(x−x0)μ2+m2−ω2] m+ω.(4.19) Letting x0=0, we have that tan θ=μ+μ2+m2−ω2tanh (κxμ2+m2−ω2) m+ω.(4.20) The potential is, thus, given by V2(x)=μsin 2θ=2μtan θ 1+tan2θ=2μ[μ+μ2+m2−ω2tanh (κxμ2+m2−ω2)] (m+ω)[μ+√μ2+m2−ω2tanh (κx√μ2+m2−ω2)]2 (m+ω)2+1 .(4.21) For the values κ=1, μ=−1/4, m=1, ω=1/2, we find that V2(x) has a kinklike form shown in Fig. 3. The charge density R2corresponding to these values of the parameters in V2(x) is shown in Fig. 4. For small xfor these parameters V2(x)≃3 37 −1365x 5476 −54249x2 810448 +O(x3),(4.22) and if instead we choose μ=1/4 and leave the other parameters the same, we get the opposite type of kink shown in Fig. 5. Now for small xwe have V2(x)≃3 37 +1365x 5476 −54249x2 810448 +O(x3).(4.23) We can also solve for the more general case V3(x)=μ1cos(2θ(x)) +μ2sin(2θ(x)),(4.24) and then the solution to the differential equation is θ(x)=−κω +κ(m+μ1) cos(2θ(x)) +κμ2sin(2θ(x)),(4.25) so we just need to change m→m+μ1and μ→μ2to obtain tan θ= μ2+μ2 2+(m+μ1)2−ω2tanh κxμ2 2+(m+μ1)2−ω2 m+μ1+ω.(4.26) 046602-6 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) The potential is then given by V3(x)=μ1cos 2θ+μ2sin 2θ=μ1 1−tan2θ 1+tan2θ+2μ2 tan θ 1+tan2θ=N(x) D(x),(4.27) where N(x)= 2μ2μ2+μ2 2+(μ1+m)2−ω2tanh κxμ2 2+(μ1+m)2−ω2 μ1+m+ω +μ11−μ2+μ2 2+(μ1+m)2−ω2tanh κxμ2 2+(μ1+m)2−ω22 (μ1+m+ω)2(4.28) and D(x)=μ2+μ2 2+(μ1+m)2−ω2tanh κxμ2 2+(μ1+m)2−ω22 (μ1+m+ω)2+1.(4.29) For ω=1/2, μ1=μ2=−1/4, κ=1,m =1, V3(x) has the form of Fig. 6. The charge density corresponding to this potential is given in Fig. 7. One can also consider the potential V4=μtan2θ. In that case, on substituting t=tan θ, we obtain the differential equation for t, dt dx =κ[t2(μ−m−ω)+m+μt4−ω].(4.30) Thus, we obtain κx = √2μtan−1√2μt μ−√(μ+ω)2+m2+2m(ω−3μ)−m−ω (μ+ω)2+m2+2m(ω−3μ)μ−(μ+ω)2+m2+2m(ω−3μ)−m−ω − √2μtan−1√2μt μ+√(μ+ω)2+m2+2m(ω−3μ)−m−ω (μ+ω)2+m2+2m(ω−3μ)μ+(μ+ω)2+m2+2m(ω−3μ)−m−ω .(4.31) A. “ω” stability of exact solutions with κ=1 in the presence of an external field Here we follow the method of Bogolubsky [17] and consider stability to changes in the frequency ω, keeping the charge Qfixed. That is, if we parametrize a rest-frame 2 1 1 2 x 0.1 0.1 0.2 V x FIG. 3. (Color online) Potential V2(x) corresponding to the exact solution Eq. (4.21) with ω=1/2, μ=−1/4. solitary-wave solution by ψs(x,t)=χs(x,ω)e−iωt,(4.32) 4 2 2 4 x 0.5 1.0 1.5 x FIG. 4. (Color online) The charge density for the solitary wave in the potential V2(x)givenEq.(4.21) with ω=1/2, μ=−1/4. 046602-7 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) 2 1 1 2 x 0.1 0.1 0.2 V x FIG. 5. (Color online) Potential V2(x) corresponding to the exact solution Eq. (4.21) with μ=1/4. then we choose our perturbed wave function to be ˜ ψ[x,t,ω,ω]=√Q[ω] √Q[ω]χs(x,ω)e−iωt ≡f(ω,ω)χs(x,ω)e−iωt.(4.33) Again inserting this wave function into the Hamiltonian, we get a new probe Hamiltonian Hpdepending on both ω,ω. We will consider the most general solution we have found where the potential has the form of V3(x) given by Eq. (4.24). The energy density is given by the four terms in Eq. (4.9), (h1+h2−h3+h4), where h1=R2dθ dx;h2=mR2cos 2θ, h3=g 2(R2cos 2θ)2;h4=R2(μ1cos 2θ+μ2sin 2θ). We also have that h3=h1. Following our previous discussion, the probe Hamiltonian has the form Hp[ω,ω]=H1[ω][f(ω,ω)2−f(ω,ω)4] +(H2+H4)[ω]f(ω,ω)2.(4.34) For our choice of potential, we have the following relationship: dθ dx =(m+μ1) cos 2θ+μ2sin 2θ−ω. (4.35) 3 2 1 1 2 3 x 0.3 0.2 0.1 0.1 V x FIG. 6. (Color online) Potential V3(x) corresponding to Eq. (4.27) with ω=1/2,μ1=μ2=−1/4. 10 5 5 10 x 0.2 0.4 0.6 0.8 x FIG. 7. (Color online) Charge density for the solitary wave in the potential V3(x) corresponding to Eq. (4.27) with ω=1/2, μ1= μ2=−1/4. Solving this equation, we obtain tan θ≡T=D+Etanh Bx, (4.36) where B=(m+μ1)2+μ2 2−ω2;E=B d1 ; (4.37) D=μ2 d1 ;d1=m+μ1+ω. We also find R2=2 g2(m+μ1) cos 2θ+μ2sin 2θ−ω cos22θ =2B2 d1g2sech2Bx [1 +(D+Etanh Bx)2] [1 −(D+Etanh Bx)2]2.(4.38) The charge is given by Q[ω]=∞ −∞ dxR2= 2ω(m+μ1)2+μ2 2−ω2 g2(ω2−μ2 2).(4.39) Also, H1=−ωQ +H2+H4;H2+H4=(m+μ1)T1+μ2T2, (4.40) where T1=dxR2cos 2θ=2[tanh−1(E−D)+tanh−1(E+D)] =2 g2tanh−1(μ1+m)2+μ2 2−ω μ1+m,(4.41) T2=dxR2sin 2θ =2μ2(m+μ1) g2ω2−μ2 2−2μ2 g2 tanh−1√(μ1+m)2+μ2 2−ω2 m+μ1 (μ1+m)2+μ2 2−ω2 . (4.42) 046602-8 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) Thus, Hp[ω]=−ωQ[ω]f(ω,ω)2[1 −f(ω,ω)2] +f(ω,ω)2(H2[ω]+H4[ω])[2 −f(ω,ω)2]. (4.43) If we take the first derivative of this probe Hamiltonian with respect to ωand set ω=ω, we find that this is zero only when μ2=0. Thus, the probe Hamiltonian is stationary to ω variations only when μ2=0.However, we can use this method for the external potential V1(x)=μ1cos 2θ[see Eq. (4.14)]. In that case, the first derivative is automatically zero when ω=ω. The second derivative changes sign at a particular value of ω=ω, which, for m=1, is the solution to the equation 2(μ1+1)3coth−1μ1+1 √μ1(μ1+2)−ω2+1 μ2 1+2μ1−ω2+1 −3μ2 1−6μ1+ω2−3=0.(4.44) Nowwehavethatfortheretobeanallowedrealsolution ω<m+μ1.Soifμ1is positive, ωis shifted upward compared to ωcand so is the region in ωspace where the solutions are real. If μ1is negative ωis shifted downward, but also the allowed regime of ωis decreased. The net result is that the possible regime of stability is approximately 30% of the allowed region for real solutions independent of the value of μ1. Preliminary simulations show that some of the exact solutions with the potential (4.14) are metastable. For the well potential, μ1<0, the solitary wave is metastable, probably due to the fact that the initial solitary wave extends beyond the inflection points of the potential well. Details of these simulations will be discussed in a future paper. V. VARIATIONAL ANSATZ FOR THE NLDE IN EXTERNAL FIELDS The gauge-invariant Lagrangian for the external field problem is given by Eq. (4.3). Using the freedom of gauge invariance, one can choose the axial gauge A1=0, eA0= V(x). Our ansatz for the trial variational wave function is to assume that, because of the smallness of the perturbation, the main modification to our exact solutions to the NLDE equation without an external field is that the parameters describing the position, momentum, boost, and phase become time dependent. That is, we replace vt →q(t); η→η(t); γωv →p(t); ωt=γω(t−vx)→φ(t)−p(t)[x−q(t)],(5.1) where φ(t)=ωγ t −p(t)q(t). Thus, our trial wave function in component form is given by 1(x,t)=cosh η 2A(x)+isinh η 2B(x)e−iφ+ip(x−q), (5.2) 2(x,t)=sinh η 2A(x)+icosh η 2B(x)e−iφ+ip(x−q), where x=cosh η(t)[x−q(t)]. Using this trial wave function we can determine the effective Lagrangian for the variational parameters. Writing the Lagrangian density as L=L1+L2+L3,(5.3) where L1=i 2(¯ γμ∂μ−∂μ¯ γμ), L2=−m¯  +g2 κ+1(¯ )κ+1,(5.4) L3=−eA0¯ γ0≡−V(x)†. Integrating over xand changing integration variables to z= (x−q) cosh η, one obtains L1=dxL1 =Q(p˙ q+˙ φ−ptanh η)−I0(cosh η−˙ qsinh η),(5.5) where Q=dz[A2(z)+B2(z)] (5.6) is as given by Eq. (2.36). Note that I0=dz(BA−AB)=H1,(5.7) where H1is the rest-frame kinetic energy and is given by Eq. (2.44).HereB(x)=dB(x) dx, and L2=L2dx =− m cosh ηI1+g2 (κ+1) cosh ηI2,(5.8) where I1=dz[A2(z)−B2(z)]; I2=dz[A2(z)−B2(z)](κ+1), (5.9) and L3=−dzρ(z)Vz cosh η+q(t)=−U[η(t),q(t)]. (5.10) Putting these terms together, we obtain L=Q(p˙ q+˙ φ−ptanh η)−I0(cosh η−˙ qsinh η) −m cosh ηI1+g2 (κ+1) cosh ηI2−U[η(t),q(t)].(5.11) We now get the following Lagrange’s equations: d dt δL δ˙ φ=0→dQ dt =0→Q=const,(5.12) i.e., the charge is canonically conjugated to the phase φ.The canonical solitary-wave momentum, which is conjugated to the solitary-wave position, is Pq=δL δ˙ q=Qp +I0sinh η, dPq dt =Q˙ p+I0cosh η˙η=δL δq =−∂U ∂q , 046602-9 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) -20 -10 0 10 20 30 x 0 0.5 1 1.5 2 2.5 3 |ψ1(x,t*)|2+|ψ2(x,t*)|2 -20 0 20 40 x 0 0.5 1 1.5 2 2.5 3 |ψ1(x,t*)|2+|ψ2(x,t*)|2 0 200 400 600 800 t -5 0 5 10 15 q(t) 0 200 400 600 800 t 0 2 4 6 P(t) FIG. 13. (Color online) Harmonic potential, V(x)=(V2/2)x2, with ωin the unstable regime. (Left upper panel) Charge density ρQat t∗=0; 133.3. (Right upper panel) Charge density at t∗=150. (Lower panels) Solitary-wave position q(t) and momentum P(t) from the numerical solutions of the CC equations (black solid lines) and numerical simulations (dashed lines) of the forced NLDE. The energy (red or middle curve) and charge (blue or upper curve) are also plotted. Parameters: g=1, m=1, ω=0.3, and V2=10−4. Initial condition: exact solitary wave of the unperturbed NLDE with initial velocity v(0) =0.1. and sn(x,l) is a Jacobi elliptic function. 1. Energy conservation From Eq. (5.33) we have that the solitary-wave energy is given by E=γM 0−cos kqI4[˙ q]+kγ cos kq ˙ q2I5[˙ q].(7.36) In the nonrelativistic limit we obtain E=M0−cos kqI0 4+M0 2+kI0 5cos kq˙ q2.(7.37) In the case of a weak potential (except for M0→0 when ω→1) E=1+˙ q2 2M0−cos kqI0 4.(7.38) 2. Solitary-wave momentum and dynamical stability The solitary-wave momentum is given by Eq. (5.36) and becomes P=γ(M0+kI5[˙ q] cos kq)˙ q. (7.39) In the nonrelativistic regime we obtain P=M0+kI0 5cos kq˙ q. (7.40) The necessary condition for stability, dP d˙ q=M0+kI0 5cos kq>0,(7.41) is satisfied except in the regime where M0→0, which is when ω→1. In that regime, the solitary wave is very broad and the condition 2π/k bis not fulfilled. 3. Numerical results for q(t)and P(t) For the pendulum equation there is a critical initial velocity at which the coordinate q(t) makes a transition from periodic motion to unbounded motion. This occurs when the modulus parameter l=1. This yields the condition vc=4I0 4 M0 .(7.42) Depending on our choice of parameters, for small-enough , vcwill be in the nonrelativistic regime. We choose ωto be in the stability region for the unforced problem (see Sec. III). For g=1, m=1, ωc=0.697586 and choosing ω=0.9, then the width of the solitary wave is 1/(2β)=1.15. If we choose k=0.1, then the characteristic wavelength 2π/k =62.8 1/(2β). From Eqs. (2.47) and (2.48) we have that Q=0.968644; H1=0.0625108 =I0,(7.43) M0=H1+ωQ =0.934291.(7.44) 046602-16 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) -20 -10 0 10 20 0 0.05 0.1 0.15 0.2 |ψ1(x,t*)|2+|ψ2(x,t*)|2 0 1000 2000 3000 4000 t -4 -2 0 2 4 q(t) 0 1000 2000 3000 4000 t -0.01 -0.005 0 0.005 0.01 P(t) FIG. 14. Periodic potential, V(x)=−cos(kx), very low initial velocity. (Upper panel) Charge density ρQat t∗=0; 2666.6. (Middle panel) Solitary-wave position q(t) from the numerical solution of Eq. (7.30) (solid line), approximate analytical expression (7.35) (dotted line), and numerical simulations (dashed line) of the forced NLDE. The three curves are superimposed. (Lower panel) Momentum P(t) from the numerical solutions of Eq. (7.30) (solid line) and numerical simulations (dashed line) of the forced NLDE. The curves are superimposed. Parameters: g=1, m=1, ω=0.9, =0.001, and k=0.1. Initial condition: exact solitary wave of the unperturbed NLDE with initial velocity v(0) =0.01. The other constants for this initial condition from Eq. (7.31) are I0 4=0.94632; I0 5=0.433477.(7.45) We have, first, compared the analytical solution Eq. (7.35) of the pendulum equation with the numerical solution of Eq. (7.30).For<0.1 the results are practically identical, for ⩾1 deviations occur. -20 0 20 40 60 x 0 0.05 0.1 0.15 0.2 |ψ1(x,t*)|2+|ψ2(x,t*)|2 0 1000 2000 3000 4000 5000 t -30 -20 -10 0 10 20 30 q(t) 0 1000 2000 3000 4000 5000 t -0.06 -0.04 -0.02 0 0.02 0.04 0.06 P(t) FIG. 15. Spatially periodic potential, V(x)=−cos(kx), initial velocity just below vc. (Upper panel) Charge density ρQat t∗=0; 5000. (Middle panel) Solitary-wave position q(t) from the numerical solution of Eq. (7.30) (solid line), approximate analytical expression (7.35) (dotted line), and numerical simulations (dashed line) of the forced NLDE. Solid and dotted lines are superimposed. (Lower panel) Momentum P(t) from the numerical solution of Eq. (7.30) (solid line) and numerical simulations (dashed line) of the forced NLDE. Parameters: g=1, m=1, ω=0.9, =0.001, and k=0.1. Initial condition: exact solitary wave of the unperturbed NLDE with initial velocity v(0) =0.0626619. Specifically we have chosen the initial condition q(0) = 0,˙ q=v0for the three cases (1) v0vc1 and then v0slightly below (2) and above (3) the critical value vc, namely v0=vc∓.001.(7.46) 046602-17 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) 0 30 60 90 120 150 0 0.05 0.1 0.15 0.2 |ψ1(x,t*)|2+|ψ2(x,t*)|2 0 1000 2000 3000 4000 t 0 50 100 150 q(t) 0 1000 2000 3000 4000 t 0 0.02 0.04 0.06 P(t) FIG. 16. Spatially periodic potential, V(x)=−cos(kx), initial velocity just above vc. (Upper panel) Charge density ρQat t∗=0; 4000. (Middle panel) Solitary-wave position q(t) from the numerical solution of Eq. (7.30) (solid line), approximate analytical expression (7.35) (dotted line), and numerical simulations (dashed line) of the forced NLDE. Solid and dotted lines are superimposed. (Lower panel) Momentum P(t) from the numerical solution of Eq. (7.30) (solid line) and from numerical simulations (dashed line) of the forced NLDE. Parameters: g=1, m=1, ω=0.9, =0.001, and k=0.1. Initial condition: exact solitary wave of the unperturbed NLDE with initial velocity v(0) =0.0646619. Choosing =0.001 yields vc=0.0636619, which is in the nonrelativistic regime, so we expect Eq. (7.35) to hold. In Fig. 14 we show that for v0=0.01 the analytic nonrelativistic result and the numerical solution of Eq. (7.30) give the same results as the solution of the NLDE. We also see that the shape of the charge density does not change in time. In Figs. 15 and 16 we show that just below and above the critical velocity, respectively, the analytical result (7.35) agrees with the numerical solution of Eq. (7.30), but both results differ very slightly from the simulation results. A summary of the result of our simulations of solitary waves in different external fields is presented in Table I. VIII. CONCLUSIONS In this study we have found exact solutions to the NLDE with scalar-scalar interactions of the form g2 κ+1(¯ )κ+1in certain external fields. We have discussed their stability with respect to “ω” variations. We have also introduced a collective coordinate method for studying the time evolution of solitary waves in external fields and determined simple equations for the collective coordinates that parallel those of a relativistic point particle. We found that unless (or until) the solitary waves displayed an instability, the collective coordinates describing the position and momentum in the CC equations gave remarkably good agreement with their counterparts from our simulations. We then presented a generalization of a dynamical stability criterion, based only on solving the CC equations, that was useful in studying the stability of solitary waves in the forced NLSE problem. For our simulations of the exact evolution, as well as the evolution of the collective coordinates, we concentrated on κ=1. For the forcing terms we used simple test potentials such as ramp, harmonic, and periodic potentials. In many instances, we found that the instability of the solitary wave solution was related to the metastability of the solitary wave in the absence of external forces and the critical time for breakup was quite close to the time found for the unforced problem. We had hoped that a generalization of the method used to map out domains of instability in the NLSE using the much simpler solutions of the collective coordinate problem would also work for the NLDE equation. Unfortunately for the problems we studied, we obtained dp d˙ q>0, which fulfills a necessary condition for stability, so this method did not give any information about instabilities. What we did find using the collective coordinate approximation was that, starting with exact solutions of the unforced problem, these solitary waves maintained shape in the CC approximation apart from the parameters becoming functions of time. The collective variables in the simulations, namely q(t) and P(t), were smooth functions for a reasonable period of time, even in the case when the solitary waves were only metastable. When these collective variables began to rapidly oscillate and/or diverge from their values found in the collective coordinate calculation, then that “defined” the onset of the instability. The criterion we used for the onset of instability using the collective coordinates is a sufficient condition and we found no cases where the condition for this dynamic instability was satisfied. Possibly this is a result of the fact that external fields differ from external sources. For the external source problem, we would expect in the nonrelativistic regime that we would recover the known results for the forced NLSE with source terms due to the arguments of Comech [15]. The simulations in this paper are confined to the κ=1 case. The numerical stability of solitary waves in the absence 046602-18 NONLINEAR DIRAC EQUATION SOLITARY WAVES IN ... PHYSICAL REVIEW E 86, 046602 (2012) TABLE I. Simulation results for three simple potentials using different parameter sets; g=1andm=1. Potential Cases Results V(x)=−V1xω=0.9>ω c=0.697586 Stable soliton, width Lorentz contracted, V1=10−2,10 −3,10 −4height increases ω=0.3<ω cAsymmetric shape, metastable for t⩽100, V1=0.01 unstable for t110 ω=0.3<ω cMetastable for t⩽100, splits into V1=0.0001 two solitons and radiation for t120 V(x)=1 2V2x2ω=0.9>ω cStable soliton, harmonic oscillations V2=0.0001, v0=0.1 ω=0.9>ω c, Metastable for t350, V2=0.0001, v0=0.9 unstable for t350 ω=0.3<ω c, Metastable for t120, splits into V2=0.0001, v0=0.1 two solitons and radiation for t120 V(x)=−cos(kx)ω=0.9>ω c,k=0.1, =0.0001 Stable soliton, harmonic oscillations v0=0.01 vc=0.0636619 ω=0.9>ω c,k=0.1, =0.0001 Stable soliton, very anharmonic v0=vc−0.001 oscillations ω=0.9>ω c,k=0.1, =0.0001 Stable soliton, translational motion v0=vc+0.001 plus oscillations of external potentials for general κwill be presented in a subsequent publication [18]. The semiclassical reduction of NLDE to NLSE and implications for solitary-wave stability have been recently discussed in a rigorous fashion by Comech [15]. Our numerical findings [18] agree with his analysis in the nonrelativistic regime. ACKNOWLEDGMENTS This work was supported in part by the US Department of Energy. F.G.M. acknowledges the hospitality of the Mathematical Institute of the University of Seville (IMUS) and of the Theoretical Division and Center for Nonlinear Studies at Los Alamos National Laboratory and financial support by the Plan Propio of the University of Seville and by Junta de Andalucia. N.R.Q. acknowledges financial support from the Humboldt Foundation through Research Fellowship for Experienced Researchers SPA 1146358 STP and by the MICINN through FIS2011-24540 and by Junta de Andalucia under Projects No. FQM207, No. FQM-00481, No. P06-FQM-01735, and No. P09-FQM-4643. APPENDIX: REST-FRAME IDENTITIES In the rest frame, energy-momentum conservation leads to identities among the various integrals that arise. The Lagrangian in the axial gauge is given by L=i 2[¯ γμ∂μ−∂μ¯ γμ]−m¯  +g2 κ+1(¯ )κ+1 −V(x)†. (A1) In the rest frame, the wave function is given by 0=ψe−iωt =A(x) iB(x)e−iωt.(A2) The energy-momentum conservation is given by Eq. (2.6) and leads to two independent equations. The first is ∂0T01 +∂1T11 =0.(A3) In the rest frame T0x=0, so Txx =const. If the solution goes to zero at infinity, then the constant is zero. We then have the relationship T11 =i 2[¯ γ1∂1−∂1¯ γ1]+L =ωψ†ψ−m¯ ψψ +g2 k+1(¯ ψψ)k+1−V(x)ψ†ψ=0. (A4) Integrating over space we get the relations ωQ −mI1+g2 κ+1I2−dxρ(x)V(x)=0.(A5) The second conservation law is ∂0T00 +∂1T10 =0,(A6) which leads to the conservation of energy. The energy of the solitary wave in the rest frame defines the rest mass M0, E=T00dx =M0.(A7) We have that T00 =i 2[¯ ψγ1∂1ψ−∂1¯ ψγ1ψ] +m¯ ψψ −g2 k+1(¯ ψψ)k+1+V(x)ψ†ψ =(AB1−BA1)+m(A2−B2) −g2 k+1(A2−B2)κ+1+(A2+B2)V(x).(A8) 046602-19 MERTENS, QUINTERO, COOPER, KHARE, AND SAXENA PHYSICAL REVIEW E 86, 046602 (2012) Integrating, we obtain M0=I0+mI1−g2 k+1I2+ρ(x)V(x).(A9) Using the identity of Eq. (A5), we then have, even in the presence of interactions, that M0=I0+ωQ. (A10) [1] R. J. Finkelstein, C. Fronsdal, and P. Kaus, Phys. Rev. 103, 1571 (1956). [2] U. Enz, Phys. Rev. 131, 1392 (1963). [3] M. Soler, Phys.Rev.D1, 2766 (1970). [4] W. A. Strauss and L. Vazquez, Phys. Rev. D 34, 641 (1986). [5] S. Y. Lee, T. K. Kuo, and A. Gavrielides, Phys. Rev. D 12, 2249 (1975). [6] Y. Nogami and F. M. Toyama, Phys.Rev.A45, 5258 (1992). [7]D.J.GrossandA.Neveu,Phys.Rev.D10, 3235 (1974). [8] W. Thirring, Ann. Phys. 3, 91 (1958). [9] A. Alvarez and B. Carreras, Phys. Lett. A 86, 327 (1981). [10] F. Cooper, A. Khare, B. Mihaila, and A. Saxena, Phys. Rev. E 82, 036604 (2010). [11] N. R. Quintero, F. G. Mertens, and A. R. Bishop, Phys. Rev. E 82, 016606 (2010). [12] F. G. Mertens, N. R. Quintero, and A. R. Bishop, Phys. Rev. E 81, 016608 (2010). [13] F. G. Mertens, N. R. Quintero, I. V. Barashenkov, and A. R. Bishop, Phys.Rev.E84, 026614 (2011). [14] F. Cooper, A. Khare, N. R. Quintero, F. G. Mertens, and A. Saxena, Phys. Rev. E 85, 046607 (2012). [15] A. Comech, arXiv:1203.3859, and references therein. [16] N. G. Vakhitov and A. A. Kolokolov, Radiophys. Quantum Electron. 16, 783 (1973). [17] I. L. Bogolubsky, Phys. Lett. A 73, 87 (1979). [18] N. R. Quintero, F. G. Mertens, F. Cooper, A. Khare, and A. Saxena (unpublished). 046602-20