CFL Criterion: Stability in Numerical PDEs

• CFL Criterion is a fundamental stability condition that limits the time step in explicit PDE solvers by ensuring the numerical domain of dependence encapsulates the physical domain.
• Its derivation via von Neumann stability analysis highlights the requirement that numerical stencils must adhere to physical wave propagation constraints.
• Enforcing the CFL condition is crucial in diverse applications like CFD and electromagnetics to achieve stable, convergent, and physically accurate simulations.

The Courant-Friedrichs-Lewy (CFL) criterion is a fundamental stability condition that constrains the maximum allowable time step in explicit time integration schemes for hyperbolic and parabolic partial differential equations (PDEs). It arises from the requirement that the numerical domain of dependence must contain the physical domain of dependence—a prerequisite for both numerical stability and physical fidelity in wave propagation, transport, and diffusion equations. The CFL condition is central to the design and analysis of finite-difference, finite-volume, finite-element, and kinetic solvers across computational fluid dynamics, electromagnetics, plasma physics, and related fields.

1. Mathematical Foundation and Derivation

The CFL condition is obtained via von Neumann stability analysis applied to discretizations of the prototypical scalar wave equation or linear advection equation. For the 1D wave equation $u_{tt} = c^2 u_{xx}$, central finite differences in space and time yield the three-level leapfrog scheme: $$u^{n+1}_k = 2u^n_k - u^{n-1}_k + r^2 (u^n_{k+1} - 2u^n_k + u^n_{k-1}), \quad r = c \frac{\Delta t}{\Delta x}.$$ Stability analysis prescribes $|r| \leq 1$, i.e.

$\Delta t \leq \frac{\Delta x}{c}$

for 1D, and in 3D Cartesian meshes,

$$\Delta t \leq \left[ c \sqrt{\Delta x^{-2} + \Delta y^{-2} + \Delta z^{-2}} \right]^{-1}.$$

Thus, the maximal admissible time step is inversely determined by the maximum signal speed and the finest grid spacing.

This constraint generalizes to any explicit scheme for hyperbolic PDEs: $$CFL \equiv \max_i \left( \frac{|\lambda_i| \Delta t}{\Delta x} \right) \leq C_\text{stab},$$ where $\lambda_i$ are the maximal local characteristic wave speeds and $C_\text{stab}$ is a scheme-dependent constant (e.g., $1$ for classical leapfrog, $<1$ for higher-order TVD/SSP Runge-Kutta methods) (Li et al., 2014, Lilly et al., 2024).

2. Physical Interpretation and Practical Implications

The CFL number quantifies the distance traversed by the fastest physical signal during a single time step, relative to the local grid spacing. Enforcing the CFL condition ensures that the numerical stencil does not update with information outside the computational domain of dependence, which would otherwise render the discrete problem ill-posed and subject to instability and spurious amplification of unresolved modes (Li et al., 2014, Ikeya et al., 2014).

Fine geometric features (leading to small $\Delta x$) or spatially adaptive meshes further reduce $\Delta t$. In explicit schemes, this can make the number of time steps required for a global simulation prohibitive, especially in multiscale or multiphysics settings (Peterseim et al., 2016, Zhang et al., 2016).

3. Extensions, Generalizations, and Context-Specific Forms

The precise structure of the CFL condition depends on the PDE class, discretization, and update kinetics:

• Advection-Diffusion-Reaction Models: For explicit discretizations of Burgers' equation or the Navier-Stokes system, distinct "convective" and "diffusive" CFL numbers appear:

$$\text{advective: } \Delta t \leq \frac{\Delta x}{U_\text{max}}, \quad \text{diffusive: } \Delta t \leq \frac{\Delta x^2}{2 \nu},$$

where $U_\text{max}$ is the maximum local velocity and $\nu$ the viscosity (Potluri et al., 18 Jan 2026, Kaya et al., 2024).

• Moment Closure and Kinetic Models: Moment and kinetic solvers impose feasibility and stability via a CFL constraint tied to the maximal discrete velocity,

$\Delta t \leq \frac{\Delta x}{2 \max_i |\xi_i| },$

where $\xi_i$ are discrete velocity nodes. This condition ensures positivity and $L^2$ stability of the constrained minimization dynamics (Sarna, 2020).

• Nonlinear Conservation Laws with IRP: For systems of conservation laws (e.g., multispecies traffic, sedimentation), explicit WENO-SSP-RK methods require

$$\Delta t \leq \mu \hat w_1 \frac{\Delta x}{\alpha},$$

where $\alpha$ is the maximum local signal speed, $\hat w_1$ a Gauss–Lobatto quadrature weight, and $\mu$ a flux-dependent factor (Barajas-Calonge et al., 4 Jun 2025).

A summary table for several explicit CFL constraints:

Problem Type	Explicit CFL Limit	Reference
3D FDTD electromagnetics	$$\Delta t \leq [c\sqrt{\Delta x^{-2}+\Delta y^{-2}+\Delta z^{-2}}]^{-1}$$	(Li et al., 2014)
1D Advection	$\Delta t \leq \Delta x/U_\text{max}$	(Potluri et al., 18 Jan 2026)
1D Diffusion	$\Delta t \leq \Delta x^2/(2\nu)$	(Potluri et al., 18 Jan 2026)
IRP WENO schemes	$\Delta t \leq \mu \hat w_1 \Delta x/\alpha$	(Barajas-Calonge et al., 4 Jun 2025)
Moment closure (Boltzmann)	$\Delta t \leq \Delta x/(2 \max |\xi_i|)$	(Sarna, 2020)

4. Strict Enforcement in Local Time Stepping and Multiscale Schemes

In spatially heterogeneous or adaptively refined meshes, strict CFL enforcement presents a challenge. Local time-stepping (LTS) schemes permit different regions to be updated with locally maximal stable time steps, but can violate the global CFL condition if not properly synchronized. Weak CFL violations often induce non-physical delays, phase errors, or interface artifacts, even if the simulation remains linearly stable (Gnedin et al., 2018). To guarantee physical wave speeds and convergence, full global synchrony—requiring that each cell step finishes simultaneously with all immediate neighbors—is mandatory. Algorithmically, this involves neighbor communication (e.g., "pull-back" and "step-drop" messaging) and time-step quantization (Gnedin et al., 2018).

Locally optimized explicit time integration (e.g., multi-level LTS with CFL-optimized Runge-Kutta schemes such as FB-LTS) raises the effective $C_\text{stab}$ (to $\sqrt{8} \approx 2.828$ for FB-RK(3,2)) and couples coarse–fine regions via conservative predictor-corrector interfaces, delivering significant speedup while preserving both conservation and stability (Lilly et al., 2024).

5. Methods for Relaxing or Bypassing the CFL Constraint

Several advanced approaches aim to mitigate the restrictive nature of the classical CFL bound:

• Model Order Reduction (MOR): Structure-preserving reduction (e.g., block-diagonal SPRIM) projects the full FDTD system onto a much lower-dimensional explicit model. For time steps beyond the standard CFL bound, singular-value perturbation is applied to enforce Schur complement-based stability ($\sigma_i < 2/\Delta t$) in the reduced system. This explicit MOR+perturbation approach enables $\Delta t$ up to 9x the classical limit with minimal error and robust stability (Li et al., 2014, Zhang et al., 2016).
• Spectral/S-fraction ROMs: Spectrally accurate reduced-order models constrain the operator spectrum to physically relevant bands, allowing stable time steps close to the Nyquist limit ($\Delta t \sim 1/(2f_\text{max})$), vastly outperforming traditional FDTD or FV schemes (Druskin et al., 2014).
• Pseudo-spectral Analytical Time-Domain (PSATD): Analytical field updates in Fourier space inherently relax the CFL restriction; when coupled with time-averaging to suppress temporal aliasing (NCI), these schemes enable effectively CFL-free operation even in highly anisotropic, Lorentz-boosted PIC environments (Shapoval et al., 2021).
• Implicit and Semi-Implicit Solvers: Fully implicit (e.g., BTCS for Burgers' equation) and θ-implicit FDTD (optimal at CFL=1.0) can eliminate or greatly relax the explicit CFL restriction, albeit at the expense of solving global linear or nonlinear systems at every step (Ikeya et al., 2014, Potluri et al., 18 Jan 2026).

6. CFL-Conditioned Stability in Kinetics, Two-Phase Flows, and Advanced Moment Methods

• Unified Gas-Kinetic Schemes (UGKS): Unlike traditional schemes, the UGKS decomposes the update into strong-stability-preserving (SSP) convex combinations of hyperbolic sub-methods, each with explicit CFL constraints but without coupling to the collision time or diffusive limit. The resulting stability condition retains only hyperbolic-type bounds:

$$\max \left\{ \max_{k}|c_k|, \frac{v a^2 + 2\theta}{\mathrm{erf}(a/\sqrt{v\theta})} \right\} \frac{\Delta t}{\Delta x} \leq 1.$$

No parabolic or subcollision restrictions are present (Chen et al., 1 May 2025).

• IMPES for Two-Phase Flow: Generalized CFL timestep criteria employing numerically differentiated wave speeds deliver robust, regime-adaptive time stepping for compressible and incompressible flows, outperforming classical Coats-type criteria, especially in capillary-driven cases (Burr et al., 1 Oct 2025).

7. Practical Performance and Trade-Offs

Empirical tests across electromagnetic, fluid, and kinetic domains confirm that systematically relaxing or optimizing the CFL condition—when done with guaranteed stability—enables substantial speedups (3–10×), minimal loss of physical accuracy (sub-1% error in resonances or S-parameters), and improved efficiency in multiscale and multi-physics problems (Li et al., 2014, Zhang et al., 2016, Lilly et al., 2024). In contrast, weak enforcement or neglect of global CFL synchronization may yield only modest speedups while introducing phase artifacts and loss of physical fidelity (Gnedin et al., 2018).

Moreover, higher-order and invariant-region-preserving discretizations impose correspondingly stricter (smaller) CFL constants, dictated by quadrature weights and nonlinear flux constraints (Barajas-Calonge et al., 4 Jun 2025).

• “Structure-Preserving Reduction of Finite-Difference Time-Domain Equations with Controllable Stability Beyond the CFL Limit” (Li et al., 2014)
• “A Stable FDTD Method with Embedded Reduced-Order Models” (Zhang et al., 2016)
• “Enforcing the Courant-Friedrichs-Lewy Condition in Explicitly Conservative Local Time Stepping Schemes” (Gnedin et al., 2018)
• “Stability Property of Numerical Cherenkov Radiation and its Application to Relativistic Shock Simulations” (Ikeya et al., 2014)
• “Explicit and Implicit Finite Difference Solvers Implemented in JAX for Shock Wave Physics” (Potluri et al., 18 Jan 2026)
• “Overcoming timestep limitations in boosted-frame Particle-In-Cell simulations of plasma-based acceleration” (Shapoval et al., 2021)
• “S-fraction multiscale finite-volume method for spectrally accurate wave propagation” (Druskin et al., 2014)
• “Invariant-region-preserving WENO schemes for one-dimensional multispecies kinematic flow models” (Barajas-Calonge et al., 4 Jun 2025)
• “A Positive and Stable L2-minimization Based Moment Method for the Boltzmann Equation of Gas dynamics” (Sarna, 2020)
• “Local Time-Stepping for the Shallow Water Equations using CFL Optimized Forward-Backward Runge-Kutta Schemes” (Lilly et al., 2024)
• “Error analysis of a pressure correction method with explicit time stepping” (Kaya et al., 2024)
• “A new Timestep Criterion for the Simulation of Immiscible Two-Phase Flow with IMPES Solvers” (Burr et al., 1 Oct 2025)
• “Stability of the first-order unified gas-kinetic scheme based on a linear kinetic model” (Chen et al., 1 May 2025)
• “Relaxing the CFL condition for the wave equation on adaptive meshes” (Peterseim et al., 2016)