Klein–Gordon Equation Numerics

The Klein–Gordon equation requires both a field and its time derivative as initial data. A solver must resolve its positive- and negative-frequency branches and preserve the appropriate discrete energy or charge. The squared field amplitude is not a general conserved probability density; the Klein–Gordon equation explains its physical interpretation.

Required background. Spatial discretization supplies grid operators and numerical stability supplies amplification-factor analysis.

A semidiscrete second-order system

For the free scalar field, write

$$\ddot\phi=-A_h\phi, \qquad A_h=-c^2\Delta_h+\omega_m^2I, \qquad \omega_m=\frac{mc^2}{\hbar}.$$

With periodic or suitable closed boundaries, choose $-\Delta_h$ nonnegative and self-adjoint in the quadrature metric $W$. Then $A_h$ is positive definite for $m>0$; a massless problem can have a zero mode. Introducing $\pi=\dot\phi$ gives the first-order system

$$\frac d{dt}\begin{pmatrix}\phi\\\pi\end{pmatrix} = \begin{pmatrix}0&I\\-A_h&0\end{pmatrix} \begin{pmatrix}\phi\\\pi\end{pmatrix}.$$

Its generator is not ordinarily skew-Hermitian in the naive Euclidean norm of $(\phi,\pi)$, whose components even have different units. It does preserve the quadratic energy

$$\mathcal E_h= \frac12\left(\pi^\dagger W\pi+ \phi^\dagger WA_h\phi\right)$$

for the continuous semidiscrete evolution. For positive-definite $A_h$, $\sqrt{2\mathcal E_h}$ is a norm on the pair $(\phi,\pi)$. If $A_hv=0$, the component along $v$ instead obeys $\ddot\phi_v=0$ and can grow linearly, $\phi_v(t)=\phi_v(0)+t\pi_v(0)$. The energy controls its velocity but not its displacement. This is physical zero-mode motion, not by itself a numerical instability.

Frequency branches and charge

For an eigenmode $A_hv=\omega_h^2v$ with $\omega_h>0$, the two time dependences are $e^{\mp i\omega_ht}$. Pure positive-frequency initial data satisfy $\pi_0=-iA_h^{1/2}\phi_0$, computed by a spectral or other justified operator action. Setting $\pi_0=0$ generally creates an equal mixture of the nonzero-frequency branches, not a stationary positive-energy particle. Treat any zero-frequency component separately.

Up to a chosen overall normalization, the conserved complex-field charge is

$$Q_h=i(\phi^\dagger W\pi-\pi^\dagger W\phi).$$

A positive-frequency mode has $Q_h=2\omega_h\|\phi\|_W^2>0$ and a negative-frequency mode has the opposite sign. A real field has zero value of this complex-field charge while retaining positive energy. Monitoring $\sum_j|\phi_j|^2$ alone therefore tests the wrong invariant.

Leapfrog stability includes the mass term

The centered time update is

$$\phi^{n+1}-2\phi^n+\phi^{n-1} =-\Delta t^2A_h\phi^n.$$

For each nonzero-frequency spatial mode, stability requires $\Delta t\,\omega_h<2$ with a margin; equality gives a repeated endpoint amplification root that can produce spurious linear growth. At $\omega_h=0$, the update instead reproduces the physical linear solution described above. On a periodic $d$-dimensional uniform grid of spacing $a$,

$$\omega_{h,\max}^2 \leq\omega_m^2+\frac{4dc^2}{a^2},$$

so a sufficient bound is $d(c\Delta t/a)^2+(\omega_m\Delta t/2)^2<1$. The familiar massless wave-equation CFL limit is incomplete for a massive field.

Start the scheme consistently with $\phi^1=\phi^0+\Delta t\,\pi^0-\frac12\Delta t^2A_h\phi^0+O(\Delta t^3)$. A first-order starting step can spoil the intended accuracy.

Plane-wave convergence test

In one dimension,

$$\omega_h(k)^2=\omega_m^2+ \frac{4c^2}{a^2}\sin^2(ka/2),$$

and the leapfrog time frequency obeys $\sin^2(\Omega\Delta t/2)=\Delta t^2\omega_h^2/4$. Compare measured phase accumulation with this fully discrete relation, then refine $a$ and $\Delta t$ independently toward $\sqrt{\omega_m^2+c^2k^2}$.

Use a near-grid-edge mode to test the actual stability bound; a smooth packet may fail to excite an unstable high-frequency mode until roundoff eventually does. Report energy and charge diagnostics with the discrete-in-time invariant appropriate to the chosen scheme, since the displayed semidiscrete energy need not be exactly preserved by leapfrog.

Interaction and interpretation boundaries

An electrostatic potential enters through covariant time derivatives and generally introduces additional time-derivative terms. Adding an arbitrary $V(x)$ to $A_h$ does not reproduce minimal electromagnetic coupling. Particle creation, antiparticle occupation, and vacuum observables require a field-theoretic model beyond this free initial-value solver.

Exercises

For a spatially uniform massive mode, what timestep restriction remains even as spatial gradients vanish?

Solution

The mode still oscillates with $\omega_m=mc^2/\hbar$, so leapfrog requires $\Delta t\,\omega_m<2$. A criterion based solely on $c\Delta t/a$ misses this rest-energy oscillation. Removing a fast carrier by a controlled reformulation can help, but it changes the evolution variables and must be derived.

References

• LeVeque, Randall J., Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM (2007), author’s book resources — modal stability, wave equations, and convergence.
• Thaller, Bernd, The Dirac Equation, Springer (1992), author’s book page — relativistic spectral branches and fixed-particle interpretation.