Speed-of-light pulses in the massless nonlinear Dirac equation with a potential

PHYSICAL REVIEW E 96, 052219 (2017) Speed-of-light pulses in the massless nonlinear Dirac equation with a potential Niurka R. Quintero* Department of Applied Physics, EPS, Virgen de África 7, 41011, Sevilla, Spain and Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, 18071 Granada, Spain Franz G. Mertens† Physikalisches Institut, Universität Bayreuth, D-95440 Bayreuth, Germany Fred Cooper‡ Santa Fe Institute, Santa Fe, New Mexico 87501, USA and Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA Avadh Saxena§ Theoretical Division and Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA A. R. Bishop Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Received 23 May 2017; revised manuscript received 27 September 2017; published 22 November 2017) We consider the massless nonlinear Dirac (NLD) equation in 1 +1 dimension with scalar-scalar self-interaction g2 2(¯ )2in the presence of three external electromagnetic real potentials V(x), a potential barrier, a constant potential, and a potential well. By solving numerically the NLD equation, we find different scenarios depending on initial conditions, namely, propagation of the initial pulse along one direction, splitting of the initial pulse into two pulses traveling in opposite directions, and focusing of two initial pulses followed by a splitting. For all considered cases, the final waves travel with the speed of light and are solutions of the massless linear Dirac equation. During these processes the charge and the energy are conserved, whereas the momentum is conserved when the solutions possess specific symmetries. For the case of the constant potential, we derive exact analytical solutions of the massless NLD equation that are also solutions of the massless linearized Dirac equation. Decay or growth of the initial pulse is also predicted from the evolution of the charge for the case of a non-zero imaginary part of the potential. DOI: 10.1103/PhysRevE.96.052219 I. INTRODUCTION The relativistic generalization of the nonlinear Schrödinger (NLS) equation, namely the nonlinear Dirac (NLD) equation has emerged as a natural model in many physical systems, such as extended particles [1,2], light solitons in waveguide arrays and experimental realization of an optical analog for relativistic quantum mechanics [3], Bose-Einstein condensates in honeycomb optical lattices [4–6], and phenomenological models of quantum chromodynamics [7], among many others. Solitary waves in the 1 +1 dimensional nonlinear Dirac (NLD) equation have been studied [8,9] in the past in the case of massive Gross-Neveu [10] (with N=1, i.e., just one localized fermion) and massive Thirring [11] models. Soler [12] proposed in 1970 that the self-interacting four- Fermi theory was a useful model to study extended fermions. Subsequently, Strauss and Vázquez [13] were able to study the stability of this model under dilatation and found the domain of stability for the Soler solutions. These solutions are solitary *[email protected] †[email protected] ‡[email protected] §a[email protected]v [email protected] waves which can have either one or two humps, depending on the value of the frequency ω∈(0,1). Recent studies using a split operator method suggested that all stable NLD solitary waves have a one-hump profile, but not all one-hump waves are stable, while all waves with two humps are unstable [14]. In particular for the scalar-scalar self-interaction (g2/2)( ¯ )2 the solitary waves were stable in simulations only if ω∈ [0.56,1) [14]. The interaction between solitary waves of different initial charge was studied in detail for the scalar-scalar case in the work of Alvarez and Carreras [15] by Lorentz boosting the static solutions and allowing them to scatter. More accurate simulations have been performed by Shao and Tang in [16], where a new quasi-stable long-lived oscillating bound state from the binary collisions of a single-humped soliton and a two-humped soliton was observed. The dynamics of single solitons also has been studied for the NLD equation with external electromagnetic fields as well as under forcing conditions [9,17,18]. The functional shape of the initial soliton does not change, however the soliton is accelerated in the ramp potential, or its center oscillates around the initial position in the cases of harmonic and periodic potentials. For smaller values of ωthe soliton, after some transient time, is again unstable. In this paper we study the very interesting behavior of the time evolution of solitary wave solutions of the Soler model when we add a potential barrier, a constant potential or a 2470-0045/2017/96(5)/052219(10) 052219-1 ©2017 American Physical Society QUINTERO, MERTENS, COOPER, SAXENA, AND BISHOP PHYSICAL REVIEW E 96, 052219 (2017) potential well, where the evolution is governed by the massless NLD equation in that external potential and is controlled by the initial conditions (ICs). For that problem we will numerically show that the initial pulse can: (i) propagate unidirectionally, (ii) rapidly split into two pulses moving in opposite directions, or (iii) focus before splitting. Moreover, in all these cases, the resulting pulses travel with the speed of light and the charge Qas well as the energy Eare conserved. We find that the self-interaction goes to zero after a critical time and one is left with two solutions of the massless linear Dirac equation having solitary wave shape, traveling at the speed of light in opposite directions. For the case of a constant external potential we find exact solutions to the full massless NLD equation which are also exact solutions of the massless linear Dirac equation. These analytic solutions provide insight into our numerical results as they represent a single pulse at time zero which becomes two pulses moving in opposite directions at later times. We remark that the linearized massless NLD equation with a constant potential can be transformed to the (1 +1)-dimensional wave equation (i.e., the d’Alembert equation) using light cone coordinates. The general solution consists of two pulses of arbitrary shapes which move in opposite directions with the speed of light [19]. The massless Dirac equation is of interest in its own right, e.g., graphene harbors massless Dirac fermions [20]. Similarly, the general solution of the massless Maxwell-Dirac equations in 1 +1 dimensions corresponds to two wave packets propagating in opposite directions with the speed of light without changing shape, which is useful in the study of concentrated propagating waves [21]. This paper is organized as follows. In Sec. II, we present the numerical solutions of the massless NLD equation using different initial conditions for: (a) a potential barrier, (b) a constant potential, and (c) a potential well. In Sec. III we find the exact analytical solutions of the massless NLD equation with a constant external potential. In Sec. IV we discuss our main findings and conclusions. The conservation of the charge and the energy is discussed in the Appendix, where we also show that if the external potential is symmetric, V(x)=V(−x), for certain symmetries of the massless NLD equation the momentum is also conserved. II. NUMERICAL SOLUTIONS FOR THE MASSLESS NLD EQUATION The NLD equation in 1 +1 dimensions with scalar-scalar self-interaction is given by iγμ∂μ−m +g2(¯ )=0,(1) where (x,t)=ψ(x,t) χ(x,t)(2) is a two-component spinor field, ¯ =†γ0is the adjoint spinor, γμare the Dirac matrices. We choose the representation γ0=σ3and γ1=iσ2, where σjare the Pauli matrices, mis the mass and gis the coupling constant. We add electromagnetic interactions through the gauge covariant derivative i∂μ→(i∂μ−eAμ). (3) Using the freedom of gauge invariance, we choose the axial gauge A1=0, eA0=V(x). In this gauge the nonlinear Dirac equation becomes iγμ∂μ−m +g2(¯ )=γ0V(x). (4) As initial condition (IC) for the numerical solution of Eq. (4) we take the exact static solitary wave solution of Eq. (1)[17,22] ψ(x,0) =A(x),χ(x,0) =iB(x),(5) with A(x)=√2β1√m+ω g cosh β1x m+ωcosh 2β1x,(6) B(x)=√2β1√m−ω g sinh β1x m+ωcosh 2β1x,(7) where β1=√m2−ω2. We cannot set m=0, because β1 would be imaginary. Therefore we replace min Eqs. (6), (7) by a new parameter μ. That is, we use as IC Eq. (5) with A(x)=√2β√μ+ω g cosh βx μ+ωcosh 2βx,(8) B(x)=√2β√μ−ω g sinh βx μ+ωcosh 2βx,(9) where β=μ2−ω2. This IC has essentially the same properties as Eqs. (6), (7): it is localized, A(x) is symmetric and B(x) is antisymmetric. As we want to investigate the massless NLD equation, we now set m=0inEq.(4). Since the fields can be scaled, without loss of generality we can consider g=1. Moreover, we choose μ=1 which means that ωis in the range 0 <ω<1. For the potential in Eq. (4) we consider three cases: a potential barrier V1(x)=ω+μsech 2βx, a constant potential V2(x)=ω, and a potential well V3(x)=ω−μsech 2βx.For all three cases we have included the constant term ω, because it will turn out that the energy, E=dx T00, in the second case is equal to ωQ. The energy and the charge Qare conserved (see the Appendix). Here, T00 denotes the density of the energy (in the energy-momentum tensor notation). Our numerical simulations have been performed by means of a fourth-order Runge-Kutta method. We choose N+1 points starting at n=0 and vanishing boundary conditions (±L,t)=0. The other parameters related with the discretization of the system are x∈[−L,L], x =0.02, L= 100, and t =0.0001. We first choose the parameter ω=0.9 for which the initial charge density ρ(x,0) =|ψ(x,0)|2+|χ(x,0)|2=A(x)2+B(x)2(10) is a single pulse. The simulation shows that this initial pulse splits into two pulses which move in opposite directions with the speed of light (Fig. 1). Interestingly, this scenario is the same for the three potentials used for Fig. 1. The splittings 052219-2 SPEED-OF-LIGHT PULSES IN THE MASSLESS . . . PHYSICAL REVIEW E 96, 052219 (2017) -15 -10 -5 0510 15 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (a) -60 -40 -20 0 20 40 60 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (b) -15 -10 -5 0510 15 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (c) -60 -40 -20 0 20 40 60 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (d) -15 -10 -5 0510 15 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (e) -60 -40 -20 0 20 40 60 x 0 0.05 0.1 0.15 0.2 |ψ(x,t*)|2+|χ(x,t*)|2 (f) FIG. 1. Snapshots of the charge density. (a), (c), (e): t=0 (black solid lines), t=2 (red dashed lines), and t=4 (blue dotted lines). (b), (d), (f): t=0 (black solid lines), t=20 (red dashed lines), and t=40 (blue dotted lines). (a), (b): potential barrier V1(x)=ω+μsech 2βx. (c), (d): constant potential V2(x)=ω. (e), (f): potential well V3(x)=ω−μsech 2βx. Parameters: g=1, ω=0.9, and μ=1. IC: Eq. (5) with Eqs. (8), (9). take place for short times 0 <t<t s≈10 and differ for the different potentials (Fig. 1). Notice that the NLD Eq. (4) with m=0 is invariant under the transformation ψ(x,t)→ψ(−x,t), χ(x,t)→−χ(−x,t), and x→−x. Therefore, ψ(x,t)=ψ(−x,t),χ(x,t)=−χ(−x,t).(11) These symmetries are fulfilled by the numerical solutions in Fig. 1, since the initial conditions (5) with Eqs. (8), (9)also satisfy Eq. (11). This is observed for short and long times in Figs. 2and 3. Figure 3shows an additional feature: for tts (tsis a transient time), i.e., for the red pulses at t=20 and the blue pulses at t=40, the following holds: ψ(x,t)=sign(x)χ(x,t).(12) As a consequence, ¯  =|ψ(x,t)|2−|χ(x,t)|2→0.(13) This means that the nonlinear term in Eq. (4) vanishes. This fact has also been observed in the simulations of Fig. 1for both the potential well and the potential barrier. As the nonlinear term of the NLD equation approaches zero for tts, the pulses which travel to the right and to the left with speed of light are solutions of the linearized NLD equation for ttsand will be given in the next section. When we choose ωμ=1, e.g., ω=0.1, the initial charge density Eq. (10) exhibits two humps. This facilitates the splitting into two pulses. Therefore the transient time ts=5is much shorter here than in the case ω=0.9. Defining f(x,t)=ψ(x,t)+χ(x,t),(14) h(x,t)=ψ(x,t)−χ(x,t),(15) 052219-3 QUINTERO, MERTENS, COOPER, SAXENA, AND BISHOP PHYSICAL REVIEW E 96, 052219 (2017) -20 -15 -10 -5 0510 15 20 x -0.2 0 0.2 0.4 ψr(x,t*) (a) -20 -10 0 10 20 x -0.4 -0.3 -0.2 -0.1 0 0.1 ψi(x,t*) (b) -20 -15 -10 -5 0510 15 20 x -0.2 -0.1 0 0.1 0.2 χr(x,t*) (c) -20 -10 0 10 20 x -0.2 -0.1 0 0.1 0.2 χi(x,t*) (d) FIG. 2. Snapshots of the real [(a), (c)] and imaginary [(b), (d)] parts of the spinor components ψand χ. Black solid lines: t=0, red dashed lines: t=2, blue dotted lines: t=4. Constant potential V2(x)=ω. Parameters and ICs: same as in Fig. 1. another interesting implication of Eq. (12) is that f(x,t)=0,x<0,(16) h(x,t)=0,x>0.(17) As an example we present a snapshot of the real and imaginary parts of f(x,40) for all xin Fig. 4. Interestingly, the pulse in Fig. 4resembles a semi-compacton [23], because it practically vanishes for x<0 and on the right-hand side it vanishes exponentially. So far we have used different potentials, but always the same initial condition, namely Eq. (5) with Eqs. (8), (9), and for ttswe obtained always the same scenario for the dynamics of the pulses. The splitting of the initial pulse into two pulses, which again move in opposite directions, can also be obtained using other ICs, provided that the charge density of the ICs is a localized function. For instance, one can choose ψ(x,0) =2a1sech βx, χ(x,0) =0,(18) where a1=√(μ−ω)/2 instead of Eq. (5) with Eqs. (8), (9). We numerically verified (not shown here) that the charge densities |ψ(x,t)|2and |χ(x,t)|2approach each other which means that the nonlinear term in the NLD equation vanishes for tts. In this case, ψ(x,t) and χ(x,t) also possess the symmetries in Eq. (11). However, for other localized ICs, a collision of two traveling waves can be observed before they split. Indeed, in Fig. 5, for a constant potential, starting with the initial conditions ψ(x,0) =√2a 2eiπ/4(sech[b(x+x0)] −isech[b(x−x0)]),(19) χ(x,0) =√2a 2eiπ/4(sech[b(x+x0)] +isech[b(x−x0)] ),(20) the evolution of the charge density is shown for a=1, b=2, and x0=10. The wave, initially centered at x=−x0travels to the right, whereas the one centered at x=x0travels to the left. They meet at x=0 when t=10 and then separate without changing their shapes. III. EXACT ANALYTICAL SOLUTIONS OF THE MASSLESS NLD EQUATION In the previous section we have used three different potentials and also different ICs and obtained very similar results for tts. Therefore, we now concentrate on the simplest case, namely a constant potential V2(x)=ω. The two components of the NLD Eq. (4) with m=0 satisfy i∂tψ+i∂xχ+g2{|ψ|2−|χ|2}ψ=ωψ, (21) −i∂tχ−i∂xψ+g2{|ψ|2−|χ|2}χ=−ωχ. (22) 052219-4 SPEED-OF-LIGHT PULSES IN THE MASSLESS . . . PHYSICAL REVIEW E 96, 052219 (2017) -60 -40 -20 0 20 40 60 x 0 0.1 0.2 0.3 0.4 ψr(x,t*) (a) -60 -40 -20 0 20 40 60 x -0.05 0 0.05 0.1 0.15 0.2 ψi(x,t*) (b) -60 -40 -20 0 20 40 60 x -0.1 0 0.1 χr(x,t*) (c) -60 -40 -20 0 20 40 60 x -0.2 -0.1 0 0.1 0.2 χi(x,t*) (d) FIG. 3. Same as Fig. 2, but for longer times: black solid lines: t=0, red dashed lines: t=20, blue dotted lines: t=40. Using the transformations (14) and (15) we obtain ∂tf+∂xf−ig2 2{fh +fh}h=−iωf, (23) ∂th−∂xh−ig2 2{fh +fh}f=−iωh. (24) First, we seek localized solutions of the linearized Eqs. (23), (24): ∂tf+∂xf+iωf =0,(25) ∂th−∂xh+iωh =0 (26) -100 -50 0 50 100 x -0.2 0 0.2 0.4 f(x,t*) FIG. 4. Snapshots of fr(x,t) (solid line) and fi(x,t) (dashed line) at t=t=40. Constant potential: V(x)=ω=0.9, same parameters andICasinFig.1. with ICs f(x,0) =ρ1(x)eiθ1(x),h(x,0) =ρ2(x)eiθ2(x),(27) where the amplitudes ρ1and ρ2are real functions, which go to zero for x→±∞, and the phases θ1and θ2are also real functions. Equations (25), (26) are decoupled and have the solutions f(x,t)=e−iωtρ1(x−t)eiθ1(x−t), h(x,t)=e−iωtρ2(x+t)eiθ2(x+t).(28) Substituting Eqs. (28)inEqs.(14), (15) we obtain ψ(x,t)=1 2e−iωt[ρ1(x−t)eiθ1(x−t) +ρ2(x+t)eiθ2(x+t)],(29) χ(x,t)=1 2e−iωt[ρ1(x−t)eiθ1(x−t) −ρ2(x+t)eiθ2(x+t)].(30) This is the solution of the linearized massless NLD equation. Notice that the linearized massless NLD equation with a constant potential can be transformed to the (1 +1)- dimensional d’Alembert equation by invoking light cone coordinates [19]. The solution of the latter equation comprises two arbitrarily shaped pulses moving in opposite directions. Our solution [Eqs. (29), (30)] is composed of two localized waves with the shapes ρ1(x) and ρ2(x), which travel in opposite directions with the speed of light, namely unity. For the particular cases ρ1=0orρ2=0, the initial pulse propagates unidirectionally with the speed of light. For complex constant potentials, i.e., ωis a complex constant, 052219-5 QUINTERO, MERTENS, COOPER, SAXENA, AND BISHOP PHYSICAL REVIEW E 96, 052219 (2017) (b) (a) -20 -15 -10 -5 0510 15 20 x 0 0.5 1 1.5 2 |ψ(x,t*)|2+|χ(x,t*)|2 FIG. 5. Parameters: m=0, g=1, and V(x)=0.1. IC: using Eqs. (19), (20) with a=1, b=2, and x0=10. (a): density of the charge versus x. Simulations: solid line (t=0), dashed line (t=10), and dotted line (t=20) (solid and dotted lines are superimposed). Blue and red lines, exact solutions at t=10,20, respectively; are superimposed with dashed and dotted lines. (b): contour plot of the density of charge evolution. the solution [Eqs. (29), (30)] predicts a decaying or growing function depending on the sign of the complex part of the potential (see Appendix). The solution, Eqs. (29), (30), is very similar to the numerical solutions of the full NLD equation for tts, which we obtained in Sec. II. For this time regime the nonlinear term in Eq. (4) vanished because ¯  =|ψ(x,t)|2−|χ(x,t)|2→0. As we want to find exact analytical solutions of the full NLD equation we require that the nonlinear terms in Eqs. (21), (22) and (23), (24) vanish for all times |ψ|2−|χ|2=1 2[fh +fh]=0.(31) When this condition is fulfilled, Eqs. (29), (30) are exact solutions of the full NLD equation for all times. The condition, Eq. (31), explicitly reads ρ1(x−t)ρ2(x+t) cos (x,t)=0 (32) with (x,t)=θ2(x+t)−θ1(x−t), and is fulfilled in three cases: ρ1(x−t)=0 (33) or ρ2(x+t)=0 (34) or θ2(x+t)−θ1(x−t)=±π 2.(35) For ρ1≡0, here, f(x,t)=0 and ψ(x,t)=−χ(x,t)= 1 2ei[θ2(x+t)−ωt]ρ2(x+t). In this case a pulse with the shape ρ2(x) travels to the left with the speed of light. Similarly, for ρ2≡0 a pulse travels to the right with the speed of light. For instance, assuming ψ(x,0) =χ(x,0) = exp(−x2), i.e., θ1(x)=0, ρ1(x)=2exp(−x2), and ρ2(x)=0, the solution [Eqs. (29), (30)] becomes ψ(x,t)=χ(x,t)= e−iωt exp[−(x−t)2]. Let us consider Eq. (35). As it must be fulfilled for all x and t, the only solution is that θ1and θ2are constants. Thus, (b)(a) -20 -10 0 10 20 x 0 0.5 1 1.5 2 2.5 |ψ(x,t*)|2+|χ(x,t*)|2 FIG. 6. Parameters: m=0, g=1, and V(x)=0.1. IC: using Eqs. (42), (43) with a=1, b=2, and x0=10. (a): density of the charge versus x. Simulations: solid line (t=0), dashed line (t=10), and dotted line (t=25) (solid and dotted lines are superimposed). Blue and red lines, exact solutions at t=10,25, respectively, are superimposed with dashed and dotted lines. (b): contour plot of the density of charge. 052219-6 SPEED-OF-LIGHT PULSES IN THE MASSLESS . . . PHYSICAL REVIEW E 96, 052219 (2017) -30 -20 -10 0 10 20 30 x -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 ψr(x,t*) (a) -30 -20 -10 0 10 20 30 x -1 -0.5 0 0.5 ψi(x,t*) (b) -30 -20 -10 0 10 20 30 x -0.5 0 0.5 χr(x,t*) (c) -30 -20 -10 0 10 20 30 x -0.4 -0.2 0 0.2 0.4 0.6 0.8 χi(x,t*) (d) FIG. 7. Parameters: m=0, g=1, and V(x)=0.1. IC: using Eqs. (42), (43) with a=1, b=2, and x0=10. Real [(a), (c)] and imaginary [(b), (d)] parts of ψ(x,t∗)andχ(x,t∗)versusx. Simulations: solid line (t=0), dashed line (t=10), and dotted line (t=30). Blue and red lines, exact solutions at t=10,30, respectively, are superimposed with dashed and dotted lines. choosing the minus sign in the condition (35), the other exact solution of the massless NLD equation finally reads ψ(x,t)=1 2ei(θ1−ωt)[ρ1(x−t)−iρ2(x+t)],(36) χ(x,t)=1 2ei(θ1−ωt)[ρ1(x−t)+iρ2(x+t)].(37) This solution contains an arbitrary constant θ1and arbitrary, localized functions ρ1(x,t) and ρ2(x,t), whereas the solution (29), (30) of the linearized massless NLD equation contained four arbitrary functions θ1(x,t), θ2(x,t), ρ1(x,t), and ρ2(x,t). Considering instead the plus sign in Eq. (35) is equivalent to interchanging ψ(x,t) with χ(x,t). Clearly, this interchange leaves invariant Eqs. (21), (22). To proceed it is shown that if the IC together with the massless NLD equation satisfy further symmetries, the solution given by Eqs. (36), (37) is simplified. Indeed, Eqs. (21), (22) are invariant under the transformation ψ(x,t)→−iχ(−x,t), and x→−x. This means that ψ(x,t)=−iχ(−x,t).(38) The exact solutions Eqs. (36), (37) satisfy the symmetry (38) only when ρ1(x−t)=ρ2(−x+t),ρ 1(x+t)=ρ2(−x−t).(39) In this case, the exact solution of the full NLD equation reads ψ(x,t)=1 2ei(θ1−ωt)[ρ1(x−t)−iρ1(−x−t)],(40) χ(x,t)=1 2ei(θ1−ωt)[ρ1(x−t)+iρ1(−x−t)].(41) This has an important consequence: In Eqs. (36), (37)the right and the left running pulses generally have different shapes ρ1(x) and ρ2(x), as shown in the evolution of the density of charge in Fig. 6, where the initial conditions satisfy f(x,0) =√2asech[b(x+x0)] eiπ/4,(42) h(x,0) =√2a1+tanh[b(x−x0)] cosh[b(x−x0)] e−iπ/4.(43) Although θ1is constant, ρ1(x)= ρ2(−x), therefore the symmetry (38) is not expected (see Fig. 7). However, choosing ψ(x,0) =asech[bx],χ(x,0) =ia sech[bx],(44) where aand bare constants, ρ1(x) and ρ2(x) satisfy Eq. (39), and the symmetry in Eq. (38) is fulfilled. Notice that in this case, ρ1(x,0) =ρ2(x,0) =ρ(x,0) =√2asech[bx], 052219-7 QUINTERO, MERTENS, COOPER, SAXENA, AND BISHOP PHYSICAL REVIEW E 96, 052219 (2017) -30 -20 -10 0 10 20 30 x -0.5 0 0.5 1 ψr(x,t*) (a) -30 -20 -10 0 10 20 30 x -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 ψi(x,t*) (b) -30 -20 -10 0 10 20 30 x 0 0.1 0.2 0.3 0.4 0.5 0.6 χ r (x,t*) (c) -30 -20 -10 0 10 20 30 x -0.5 0 0.5 1 χi(x,t*) (d) FIG. 8. Real [(a), (c)] and imaginary [(b), (d)] parts of the spinor components ψand χ, using the ICs given by Eq. (44). Simulation: solid lines: t=0, dashed lines: t=10, and dotted lines: t=20. Blue and red lines, exact solutions at t=10,20, respectively, are superimposed with the dashed and dotted lines. Parameters: m=0, g=1, ω=0.1, a=1, and b=2. θ1=π/4, and θ2=−π/4. The exact solution (36), (37)now becomes ψ(x,t)=1+i 2ae −iωt{sech[b(x−t)] −isech[b(x+t)]},(45) χ(x,t)=1+i 2ae −iωt{sech[b(x−t)] +isech[b(x+t)]}.(46) Figures 8and 9show the solutions represented by Eqs. (45), (46). These solutions serve as a test of our simulations, i.e., the numerical solution of the massless NLD equation, using the ICs, Eqs. (44). The simulation also confirms that the charge and the energy are the constants Q=4a2/|b|=2 and E=ωQ =0.2, respectively, predicted by the theory. Moreover, the momentum Pis zero. Other symmetries, as for instance the ones represented by Eq. (11) drastically reduce Eqs. (36), (37) to the trivial solution ρ1(x,t)=ρ2(x,t)=0, i.e., ψ(x,t)=χ(x,t)=0. IV. CONCLUSIONS In this paper we have investigated a notable property of the massless NLD equation in the presence of an external field—the fact that initial pulses centered at the origin of the external field break into two pulses traveling in opposite directions (with the speed of light), which after a short time become solutions of the massless linear Dirac equation. We have presented numerical simulations starting with various pulse initial conditions for the case of three different external potentials: a potential barrier, a potential well and a constant potential. By considering exact solutions of the massless NLD equation in the presence of a constant potential which are also solutions of the linear Dirac equation we were able to gain insight into why an initial pulse at time zero becomes two solutions of the linear Dirac equation moving in opposite directions at later times. The latter exact solution also served as a test for our numerical procedure. We note that not only does the initial pulse break into two pulses traveling with the speed of light in opposite directions but also that the initial pulse can propagate along one direction, or that two initial pulses can first focus (i.e., merge) and then split. Moreover, the pulse can either decay or grow when the potential is complex. ACKNOWLEDGMENTS We acknowledge financial support from the Ministerio de Economía y Competitividad of Spain through Grant No. FIS2014-54497-P (N.R.Q., F.G.M., and A.S.), from the Junta de Andalucía through Grant No. P11-FQM-7276 (N.R.Q.), and from the University of Seville through the Plan Propio (F.G.M.). F.G.M. acknowledges financial support and hospitality of the Theoretical Division and Center for Nonlinear Studies at Los Alamos national Laboratory. This work was supported in part by the US Department of Energy. 052219-8 SPEED-OF-LIGHT PULSES IN THE MASSLESS . . . PHYSICAL REVIEW E 96, 052219 (2017) -30 -20 -10 0 10 20 30 x 0 0.2 0.4 0.6 0.8 1 |ψ(x,t*)|2 (a) -30 -20 -10 0 10 20 30 x 0 0.2 0.4 0.6 0.8 1 |χ(x,t*)|2 (b) -30 -20 -10 0 10 20 30 x 0 0.5 1 1.5 2 |ψ(x,t*)|2+|χ(x,t*)|2 (c) (d) FIG. 9. (a), (b): Charge densities of spinor components versus x. Simulation: solid lines: t=0, dashed lines: t=10, and dotted lines: t=20. Blue and red lines, exact solutions at t=10,20, respectively, are superimposed with the dashed and dotted lines. (c): Charge density versus x. (d): contour plot of the charge density. Same parameters and ICs as in Fig. 8. APPENDIX: CONSERVATION LAWS The NLD Eq. (4) can be derived in a standard fashion from the Lagrangian density L=i 2[¯ γμ∂μ−∂μ¯ γμ]−m¯  +g2 2(¯ )2−¯ γ0V(x). (A1) Multiplying Eq. (4) on the left by ¯ , multiplying the adjoint NLD equation on the right by and summing up, we obtain the continuity equation ∂ ∂t(¯ γ0)+∂ ∂x(¯ γ1)=2VI(x)¯ γ0, (A2) where ¯ γ0=|ψ|2+|χ|2is the density of charge and VI(x) is the imaginary part of the potential. Integrating over xand assuming that vanishes at ±∞ we obtain that, for a real potential, the charge is conserved, i.e., dQ dt =∂ ∂t dx(¯ γ0)=0.(A3) Moreover, if VI(x) is constant the evolution of the charge is governed by dQ dt =2VIQ, (A4) which predicts a decay and growth of the initial pulse for VI<0 and VI>0, respectively. Multiplying Eq. (4) on the left by ¯ x, multiplying the adjoint NLD equation on the right by xand taking the difference between these two expressions, we obtain the continuity equation ∂T01 ∂t +∂T11 ∂x =¯ xγ0V(x)+¯ γ0V(x)x,(A5) where the density of momentum is T01 =i 2(¯ xγ0−¯ γ0x),(A6) T11 =−m¯  +g2 2(¯ )2+i 2[¯ γ0t−¯ tγ0].(A7) Integrating over xand assuming that T11(−∞)−T11(∞)=0 we obtain dP dt =dx[¯ xγ0V(x)+¯ γ0V(x)x],(A8) where the momentum is P=dxT01.(A9) In the case of a real constant potential V(x)=V2(x)=ω, the momentum is conserved. Moreover, if the real potential 052219-9