Discrete Klein–Gordon Equation

Updated 5 April 2026

The discrete Klein–Gordon equation is a lattice-based model describing scalar fields with nonlinear interactions and dispersion via finite-difference operators.
Its analytical framework bridges discrete spatial formulations with continuum PDE properties using Fourier methods and dispersive decay estimates.
The model underpins practical applications in lattice field theory, numerical schemes, and the study of discrete breathers and localized oscillations.

The discrete Klein–Gordon equation is a fundamental model in mathematical physics describing scalar fields on spatial lattices or in semi-discretized domains. Its analysis bridges discrete and continuum wave propagation, nonlinear dynamics, spectral theory, and numerical methods for partial differential equations. Discrete realizations appear in lattice field theory, nonlinear dynamical lattices (such as coupled mechanical oscillators and quantum crystals), and provide testbeds for bridging analysis between finite-difference numerical schemes and the functional analysis of PDEs.

1. Discretizations and Lattice Formulations

The classical Klein–Gordon equation on Euclidean space is

$\partial_t^2 u - \Delta u + m^2 u + \mathbb N(u) = f(x,t),$

where $m > 0$ is the mass, $\mathbb N(u)$ denotes nonlinearity, and $f$ is an external source. The discrete Klein–Gordon equation (DKG) arises from (i) intrinsic lattice models and (ii) finite-difference/semi-discretizations.

Lattice dKG (nearest-neighbor, regular grid)

For $n \in \mathbb{Z}^d$ ,

$\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$

Here, the discrete Laplacian is

$(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$

and the model is supplemented with suitable initial/boundary data.

Functional-Discrete Methods (FD-method, Goursat context)

For a domain discretized to mesh $x_i = i h_1$ , $y_j = j h_2$ , consider the problem

$\frac{\partial^2 u}{\partial x\,\partial y} + \mathbb N(u(x, y)) = f(x, y),\quad u(x,0) = \psi(x),\ u(0,y) = \phi(y)$

and approximate $m > 0$ 0 by “frozen” piecewise-constant values per cell ( $m > 0$ 1), yielding a set of locally linear problems for cellwise auxiliary functions $m > 0$ 2 that are then combined to form an approximate global solution (Makarov et al., 2012).

Fractional dKG

On $m > 0$ 3, the fractional Laplacian $m > 0$ 4 (for $m > 0$ 5) is constructed by convolution with fractional-centered difference kernels, leading to discrete analogues of fractional wave propagation with dispersion relations reflecting the underlying difference operator (Dasgupta et al., 2022).

Discrete Phase-Space Variant

Discrete phase-space, continuous-time models introduce “quantum” difference operators $m > 0$ 6 (weighted finite differences) and lead to Klein–Gordon-type equations where the spatial Laplacian is replaced by sums of such operators, particularly in the context of discrete quantum field theory (Das et al., 2022).

2. Linear and Nonlinear Solution Theory

Linear Discrete Propagation and Dispersive Estimates

Discrete Klein–Gordon equations exhibit mode decomposition in Fourier space,

$m > 0$ 7

with dispersion $m > 0$ 8, $m > 0$ 9 for nearest-neighbor Laplacians.

Sharp $\mathbb N(u)$ 0 dispersive decay estimates hold:

$\mathbb N(u)$ 1 in $\mathbb N(u)$ 2 (Borovyk et al., 2013)
$\mathbb N(u)$ 3 in $\mathbb N(u)$ 4, $\mathbb N(u)$ 5 in $\mathbb N(u)$ 6 (Cuenin et al., 2020)

On $\mathbb N(u)$ 7, with small quasi-periodic potential, the decay persists at $\mathbb N(u)$ 8, with explicit KAM-based reducibility arguments and Van der Corput-type stationary phase estimates (Wan et al., 27 Dec 2025).

Nonlinear Well-posedness and Scattering

For nonlinearities of the form $\mathbb N(u)$ 9, small initial data in $f$ 0 yield global solutions and scattering for thresholds depending on space dimension and $f$ 1. For instance, Strichartz estimates yield global small-data well-posedness for $f$ 2 in $f$ 3 and for $f$ 4 in $f$ 5 (Cuenin et al., 2020). Analogous results extend to the quasi-periodic potential setting (Wan et al., 27 Dec 2025).

Continuum Limit and Interpolation

For the rescaled lattice $f$ 6 and nonlinear dKG

$f$ 7

the continuum limit $f$ 8 is addressed using the Shannon interpolation operator $f$ 9 mapping $n \in \mathbb{Z}^d$ 0 to a band-limited function in $n \in \mathbb{Z}^d$ 1. For sufficiently smooth initial data, solutions $n \in \mathbb{Z}^d$ 2 converge strongly in $n \in \mathbb{Z}^d$ 3 to solutions of the continuum NLKG at rate $n \in \mathbb{Z}^d$ 4, uniformly on bounded time intervals (Chauleur, 2024).

3. Nonlinear Localized States and Breathers

Existence and Stability of Discrete Breathers

Contrary to the continuous NLKG, where the Fermi–Golden–Rule and resonance with the continuous spectrum preclude small-amplitude periodic solutions, the presence of a bounded lattice spectral band (e.g., $n \in \mathbb{Z}^d$ 5 for $n \in \mathbb{Z}^d$ 6) enables the non-resonant persistence of time-periodic, spatially localized breathers in the fully discrete system (Maeda, 2016). Under a spectral gap hypothesis and a non-resonance condition for all harmonics of the “internal mode,” analytic bifurcation and Lyapunov–Schmidt-type decompositions yield small-amplitude breathing solutions exponentially localized in space: $n \in \mathbb{Z}^d$ 7 Here, $n \in \mathbb{Z}^d$ 8 with $n \in \mathbb{Z}^d$ 9 outside $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 0, and $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 1 the ground state of the discrete Schrödinger operator.

Multiscale Reductions to Discrete NLS

For small-amplitude, weakly coupled, driven/damped lattices, the rotating-wave and multiple-scale ansatz reduce dKG to modulation equations of discrete nonlinear Schrödinger (dNLS) type. Rigorous justification of this approximation centers on energy estimates for the error dynamics and produces explicit error bounds (e.g., $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 2 for the damped, driven dKG over times $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 3 (Muda et al., 2019); similar constructions for the parametrically driven case yield $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 4 errors for $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 5 (Muda et al., 2019)). The dNLS family captures the main amplitude and phase dynamical regimes and provides insight into breather stability and instability mechanisms.

4. Geometric and Spectral Properties of Propagation

Propagation Geometry and Dispersive Zones

In 2D lattices, the shape of the “light cone” is determined by the group-velocity image,

$\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 6

which for generic parameters forms a convex compact region in $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 7 with a distinctive “astroid” inner boundary and four cusp points. Fundamental solutions decay exponentially outside this region. The polynomial decay rates for the fundamental solution, correlated to stationary phase degeneracies (nondegenerate, fold, cusp), are invariant under variation of system parameters (Borovyk et al., 2013).

Quantum Harmonic Lattice Correspondence

The quantum harmonic lattice is the canonical quantization of the classical dKG system. The decay exponents for classical propagators yield parallel decay rates for commutators of time-shifted Weyl operators, implementing Lieb–Robinson-type finite-propagation-speed bounds in the quantum setting (Borovyk et al., 2013).

5. Discrete-to-Continuum Limits and Fractional Variants

Fractional and Semiclassical Extensions

The fractional discrete Klein–Gordon equation on $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 8 replaces the nearest-neighbor Laplacian with a fractional difference operator parameterized by $\ddot u_n(t) + m^2 u_n(t) + \sum_{j=1}^d [u_{n+e_j}(t) - 2u_n(t) + u_{n - e_j}(t)] + \mathbb N(u_n(t)) = f_n(t).$ 9. The resulting models are globally well-posed in $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 0 and admit strong convergence to continuum fractional KG equations as $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 1, provided appropriate Sobolev regularity for the initial data (Dasgupta et al., 2022).

Notably, the symbol $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 2 remains bounded on the Brillouin zone and yields a maximal generator norm, in contrast to the unbounded continuum fractional Laplacian.

Discrete Phase-Space KG and Quantum Field Models

With spatial coordinates replaced by discrete quantum numbers $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 3, and evolution given by weighted difference operators $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 4, the Klein–Gordon equation is adapted to discrete phase-space representation. The Hermite function basis enables explicit (quasi-)momentum-plane-wave solutions matching continuum dispersion relations. The resulting Fock quantization, propagators, and commutators are softened in the ultraviolet and yield divergence-free Yukawa and Coulomb potentials through beta function regularization (Das et al., 2022).

6. Numerical Methods and Algorithmic Implementations

The FD-method provides a convergent algorithm for the Goursat problem for nonlinear KG equations, working via cellwise freezing and Adomian polynomial expansion of nonlinearities. Cellwise Riemann function solutions are patched to global solutions, and convergence is superexponential in the correction rank $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 5 when the step size $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 6 (Makarov et al., 2012). Numerical experiments show errors decreasing exponentially in $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 7 and proportionally to $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 8 in the basic solution, matching theoretical predictions.

7. Model Variants: Sine-Gordon, $(\Delta_{\rm disc} u)_n = \sum_{j=1}^d [u_{n+e_j} - 2u_n + u_{n-e_j}]$ 9, and $x_i = i h_1$ 0-Symmetric Systems

Discrete sine–Gordon and $x_i = i h_1$ 1 chains are archetypal cases of discrete Klein–Gordon systems. Their regimes interpolate between continuum and anti-continuum (decoupled oscillators), yielding kinks, breathers, and multibreather solutions with explicit constructions and precise bifurcation diagrams. The stability of these solutions, tracked via Gerschgorin disk arguments and Floquet multipliers, is determined by location (onsite/intersite), coupling strength, and symmetry considerations, including Peierls–Nabarro bifurcations and mode-exchange mechanisms. $x_i = i h_1$ 2-symmetric extensions demonstrate gain/loss-induced instability/stability via eigenvalue shifts in spectral analysis (Chirilus-Bruckner et al., 2014).