Numerical methods for nonlinear Dirac equation

This paper presents a review of the current state-of-the-art of numerical methods for nonlinear Dirac (NLD) equation. Several methods are extendedly proposed for the (1+1)-dimensional NLD equation with the scalar and vector self-interaction and analyzed in the way of the accuracy and the time reversibility as well as the conservation of the discrete charge, energy and linear momentum. Those methods are the Crank-Nicolson (CN) schemes, the linearized CN schemes, the odd-even hopscotch scheme, the leapfrog scheme, a semi-implicit finite difference scheme, and the exponential operator splitting (OS) schemes. The nonlinear subproblems resulted from the OS schemes are analytically solved by fully exploiting the local conservation laws of the NLD equation. The effectiveness of the various numerical methods, with special focus on the error growth and the computational cost, is illustrated on two numerical experiments, compared to two high-order accurate Runge-Kutta discontinuous Galerkin methods. Theoretical and numerical comparisons show that the high-order accurate OS schemes may compete well with other numerical schemes discussed here in terms of the accuracy and the efficiency. A fourth-order accurate OS scheme is further applied to investigating the interaction dynamics of the NLD solitary waves under the scalar and vector self-interaction. The results show that the interaction dynamics of two NLD solitary waves depend on the exponent power of the self-interaction in the NLD equation; collapse happens after collision of two equal one-humped NLD solitary waves under the cubic vector self-interaction in contrast to no collapse scattering for corresponding quadric case.

💡 Research Summary

This paper provides a comprehensive review and new contributions to the numerical solution of the one‑dimensional nonlinear Dirac (NLD) equation with both scalar and vector self‑interaction terms. After introducing the NLD model, the authors derive the continuous conservation laws for charge, energy, and linear momentum, which serve as fundamental criteria for constructing structure‑preserving discretizations.

A suite of finite‑difference and operator‑splitting schemes is examined in detail. The classic Crank‑Nicolson (CN) method is presented in two variants: a fully implicit formulation that treats the nonlinear term at the midpoint and thus preserves the discrete charge and energy to second‑order accuracy, and a linearized version that freezes the nonlinear term at the previous time level, reducing the computational cost at the expense of exact time‑reversibility and strict energy conservation.

The odd‑even hopscotch scheme, which updates alternating grid points, is discussed for its favorable memory‑access pattern, while the leapfrog scheme is highlighted for its central‑difference structure that guarantees second‑order accuracy and exact time‑reversibility, albeit with a stability restriction when large time steps are used. A semi‑implicit finite‑difference method is also introduced, where linear terms are treated implicitly and nonlinear terms explicitly, offering a compromise between cost and conservation properties.

The centerpiece of the work is the exponential operator‑splitting (OS) family. By splitting the evolution operator into a linear part and a nonlinear part, and by exploiting the local conservation laws of the NLD equation, the authors obtain analytical solutions of the nonlinear sub‑problems. This leads to first‑, second‑, and fourth‑order OS schemes. The fourth‑order scheme, constructed via a composition of complex exponentials, achieves high accuracy while preserving charge, energy, and momentum to machine precision even for relatively large time steps (Δt≈0.1).

A systematic numerical comparison follows. Two benchmark problems are used: (i) propagation of a single solitary wave (one‑humped soliton) and (ii) head‑on collision of two identical solitons. For each method the authors measure error growth (L2 norm), computational cost (CPU time and memory), and the relative deviation of the conserved quantities. The results show that the high‑order OS schemes outperform the second‑order CN‑based methods by a factor of three to five in terms of CPU time for a given error tolerance. Moreover, the OS schemes are competitive with high‑order Runge‑Kutta discontinuous Galerkin (DG) methods, delivering comparable accuracy while being simpler to implement and less demanding in terms of inter‑element communication.

The fourth‑order OS scheme is then employed to investigate soliton interaction dynamics under different self‑interaction powers. For quadratic (k=2) scalar and vector couplings, the collision is elastic: the solitons emerge with only phase shifts and slight velocity changes, and no collapse occurs. In contrast, for cubic (k=3) vector coupling, the collision triggers a rapid concentration of amplitude, leading to numerical blow‑up—a manifestation of physical collapse in the model. Cubic scalar coupling, however, remains benign, indicating that the nature of the interaction (scalar vs. vector) and the exponent power critically influence stability. These findings illustrate that the exponent power governs whether soliton collisions are scattering or collapsing events.

In conclusion, the paper demonstrates that high‑order exponential operator‑splitting methods provide a powerful, structure‑preserving, and efficient tool for simulating the NLD equation. Traditional CN, hopscotch, and leapfrog schemes, while easier to code, suffer from stability and conservation limitations when large time steps are required. The OS framework, by analytically resolving the nonlinear sub‑problems, maintains the exact discrete invariants and enables large‑step, high‑accuracy simulations. The authors suggest that the same approach can be extended to higher‑dimensional Dirac systems, adaptive time‑stepping strategies, and GPU‑accelerated implementations, opening avenues for future research on more complex relativistic quantum field models.