Stability of nonlinear Dirac solitons under the action of external potential

TL;DR

This work analyzes the stability of nonlinear Dirac solitons in the Gross-Neveu and ABS models when subjected to spatial potentials, using Lakoba's method of characteristics and a two-collective-coordinate ansatz for center-of-mass motion. By deriving an effective energy–momentum framework and Newton-like equations for the collective coordinates, the authors show that solitons are numerically stable across most regimes, with instabilities arising mainly for low-frequency modes in harmonic traps; these can be mitigated by absorbing boundary conditions in GN, though ABS instabilities at low frequencies persist. The combination of full PDE simulations and CC theory yields a consistent description of soliton dynamics, including exact linear-potential momentum evolution and harmonic-potential oscillatory behavior, aligning well with predicted trajectories and density profiles. The results reinforce the robustness of solitons under external potentials and provide a practical numerical approach for studying nonlinear Dirac solitons, with implications for stability analyses in related relativistic field theories and their perturbations.

Abstract

The instabilities observed in direct numerical simulations of the Gross-Neveu equation under linear and harmonic potentials are studied. The Lakoba algorithm, based on the method of characteristics, is performed to numerically obtain the two spinor components. We identify non-conservation of energy and charge in simulations with instabilities and we find that all studied solitons are numerically stable, except the low-frequency solitons oscillating in the harmonic potential over long periods of time. These instabilities, as in the case of Gross-Neveu equation without potential, can be removed by imposing absorbing boundary conditions. The dynamics of the soliton is in perfect agreement with the prediction obtained using an ansatz with only two collective coordinates, namely the position and momentum of the center of mass. We use the temporal variation of both field energy and momentum to determine the evolution equations satisfied by the collective coordinates. By applying the same methodology, we also demonstrate the spurious character of the reported instabilities in the Alexeeva-Barashenkov-Saxena model under external potentials.

Figures (14)

• Figure 1: Sketch of the method of characteristics applied to solve the nonlinear Dirac equation. The numerical solutions $u_{k}^{n}$ and $v_{k}^{n}$ fulfill the Eqs. \ref{['eq:uf']}-\ref{['eq:vf']}. Starting from $u_{k}^{n}$ and $v_{k}^{n}$, the equations are integrated along the characteristic lines $\bar{\xi}$ and $\bar{\eta}$, respectively. This procedure enables the solution $u_{k+1}^{n+1}$ and $v_{k-1}^{n+1}$ to be obtained.
• Figure 2: Simulations of GN model ($m=1$ and $\kappa=1$) \ref{['eq:sys']} with ramp potential, $V(x)=-V_1 x$, with $V_1=10^{-2}$. Time evolution of (a) position, (b) momentum, and the relative errors of (c) charge and (d) energy of the soliton from numerical simulations with step-size fixed at $h=10^{-1}$ (solid dark-blue line), $h=10^{-2}$ (dash-dotted red line), $h=10^{-3}$ (dotted orange line), and $h=10^{-4}$ (dashed dark-green line). (e) Comparison of soliton charge density, $\rho(x,t)=|\psi(x,t)|^2+|\chi(x,t)|^2$, at fixed times $t=50$ and $t=100$, obtained from numerical simulations with step-size fixed at $h=10^{-3}$ (solid dark-blue line) with those of the collective coordinate theory (dotted orange line). The simulated length is $L=100$ and the initial condition uses the ansatz \ref{['eq:psia']}-\ref{['eq:chia']} with $q(0)=0$, $\dot{q}(0)=0$, $p(0)=0$, and $\omega=0.9$.
• Figure 3: Simulations of GN model ($m=1$ and $\kappa=1$) \ref{['eq:sys']} with ramp potential, $V(x)=-V_1 x$, with $V_1=10^{-2}$. Comparison of soliton charge density, $\rho(x,t)=|\psi(x,t)|^2+|\chi(x,t)|^2$, at fixed times (a) $t=50$ and (b) $t=100$, obtained from numerical simulations (solid dark-blue line) with those of the collective coordinate theory (dotted orange line). (c) Time evolution of the relative errors of the soliton charge (dash-dotted red line) and energy (dashed dark-green line). Comparison of (d) position and (e) momentum of the soliton obtained from numerical simulations (solid dark-blue line) and from the collective coordinate theory (dotted orange line). The numerical step-size is fixed at $h=10^{-3}$, the simulated length is $L=100$, and the initial condition uses the ansatz \ref{['eq:psia']}-\ref{['eq:chia']} with $q(0)=0$, $\dot{q}(0)=0$, $p(0)=0$, and $\omega=0.3$.
• Figure 4: Simulations of GN model ($m=1$ and $\kappa=1$) \ref{['eq:sys']} with ramp potential, $V(x)=-V_1 x$, with $V_1=10^{-4}$. Comparison of soliton charge density, $\rho(x,t)=|\psi(x,t)|^2+|\chi(x,t)|^2$, at fixed times (a) $t=500$ and (b) $t=1000$, obtained from numerical simulations (solid dark-blue line) with those of the collective coordinate theory (dotted orange line). (c) Time evolution of the relative errors of the soliton charge (dash-dotted red line) and energy (dashed dark-green line). Comparison of (d) position and (e) momentum of the soliton obtained from numerical simulations (solid dark-blue line) and from the collective coordinate theory (dotted orange line). The numerical step-size is fixed at $h=10^{-3}$, the simulated length is $L=150$, and the initial condition uses the ansatz \ref{['eq:psia']}-\ref{['eq:chia']} with $q(0)=0$, $\dot{q}(0)=0$, $p(0)=0$, and $\omega=0.3$.
• Figure 5: Simulations of GN model ($m=1$ and $\kappa=1$) \ref{['eq:sys']} with ramp potential, $V(x)=-V_1 x$, with $V_1=10^{-5}$. Comparison of soliton charge density, $\rho(x,t)=|\psi(x,t)|^2+|\chi(x,t)|^2$, at fixed times (a) $t=1000$ and (b) $t=3000$, obtained from numerical simulations (solid dark-blue line) with those of the collective coordinate theory (dotted orange line). (c) Time evolution of the relative errors of the soliton charge (dash-dotted red line) and energy (dashed dark-green line). The numerical step-size is fixed at $h=10^{-3}$, the simulated length is $L=200$, and the initial condition uses the ansatz \ref{['eq:psia']}-\ref{['eq:chia']} with $q(0)=0$, $\dot{q}(0)=0$, $p(0)=0$, and $\omega=0.1$.