Dirac Equation Numerics

A numerical Dirac solver evolves spinor components with a discretized first-order spatial operator. Preserving its discrete inner product is essential, but not sufficient: a norm-conserving grid can still contain extra low-energy modes or misrepresent the boundary. Starting from the covariant Dirac equation, this page develops a one-dimensional Hamiltonian calculation in a prescribed external field.

Required background. Dirac current defines density and flux, and representations as arrays fixes component ordering.

A one-dimensional spinor Hamiltonian

For a two-component $1+1$ dimensional model, choose $\alpha=\sigma_x$, $\beta=\sigma_z$ and a real electrostatic potential energy $V(x)$:

$$H=-i\hbar c\,\sigma_x\partial_x +mc^2\sigma_z+V(x)I_2.$$

This is a declared reduced model, not the full four-component three-dimensional spinor written on a line. Store each grid site’s two components in a fixed order and construct derivative and component matrices consistently with that ordering.

On a periodic uniform grid, the centered derivative $D_0$ is skew-Hermitian, so $H_h=-i\hbar c(D_0\otimes\sigma_x)+I\otimes mc^2\sigma_z+V_h$ is Hermitian. With nonuniform quadrature weights $W$, the correct condition is

$$H_h^\dagger W=WH_h, \qquad \|\psi_h\|_W^2=\psi_h^\dagger W\psi_h.$$

Test that identity before propagation. A small ordinary Euclidean residual is not enough when the physical metric is weighted.

Boundaries must control the current

The one-dimensional current is $j=c\psi^\dagger\sigma_x\psi$. Integration by parts produces a boundary bilinear form. A self-adjoint closed boundary must make that form vanish for every pair of allowed states, not merely for one numerical eigenvector.

For example, imposing $\psi_2=i\lambda\psi_1$ with fixed real $\lambda$ at each endpoint gives zero normal current and a consistent separated boundary subspace. Imposing both components to vanish at both ends can overconstrain a first-order problem. Open scattering boundaries require a separate incoming/outgoing construction or a calibrated absorber.

Free dispersion and branch preparation

For a continuum plane wave, $H(k)=\hbar ck\sigma_x+mc^2\sigma_z$ has energies $E_\pm=\pm\sqrt{m^2c^4+\hbar^2c^2k^2}$. The centered grid replaces $k$ by $\sin(ka)/a$:

$$E_{\pm,h}(k)= \pm\sqrt{m^2c^4+\hbar^2c^2\frac{\sin^2(ka)}{a^2}}.$$

Check the numerical matrix eigenvalues against this lattice dispersion first, then against the continuum at fixed small physical $k$ as $a\to0$. Agreement with the lattice formula tests implementation; convergence to the continuum tests discretization. The additional low-energy region near $ka=\pi$ is fermion doubling.

To prepare a free positive-energy packet, project each resolved momentum spinor with $P_+(k)=[I+H(k)/E_+(k)]/2$ and renormalize. Use the discrete $H_h(k)$ if testing the discrete solver. At an exactly zero-energy mode, specify the degenerate subspace explicitly instead of dividing by zero. Arbitrary spinor components generally mix both branches.

Time evolution and spectral limits

A midpoint Crank–Nicolson update is unitary in the correct metric for a time-independent self-adjoint discrete Hamiltonian, up to solve errors. Converge its phase error and spatial dispersion separately; exact norm conservation does not establish correct group velocity.

Ordinary imaginary-time ground-state filtering is inappropriate for the unrestricted continuum Dirac Hamiltonian, which is unbounded below. On a finite grid it selects the most negative cutoff state. Bound-state calculations in a spectral gap need sector-aware or shift-targeted methods and checks for spectral pollution.

Strong external fields can mix positive and negative sectors. A one-particle solver does not by itself compute pair creation or a many-body vacuum response; those questions require a field-theoretic interpretation and additional state information. See the Klein paradox.

Exercises

For a massless centered-grid model, compare the energies at $k=0$ and $k=\pi/a$.

Solution

Both are zero because $\sin0=\sin\pi=0$. The second zero persists at the grid edge even though the continuum energy there would be $\hbar c\pi/a$. Norm conservation cannot distinguish this extra low-energy species from the intended one; the discretization itself must be examined.

References

• Hammer, René, and Walter Pötz, “Staggered grid leap-frog scheme for the (2+1)D Dirac equation,” arXiv (2013), arXiv:1306.5895 — structure-preserving grid design and dispersion analysis.
• Thaller, Bernd, The Dirac Equation, Springer (1992), author’s book page — Hamiltonian domains, spectral structure, and physical interpretation.