Short-Time Drift Propagator for Fokker-Planck

• Short-Time Drift Propagator is a framework that decouples the drift term from the diffusion kernel, yielding precise short-time solutions for the Fokker-Planck equation.
• It employs a pure-diffusion Gaussian convolution with an external drift shift, enabling efficient error control and preservation of density positivity.
• Systematic perturbative improvements and operator norm techniques in STDP provide enhanced accuracy for simulating stochastic trajectories in complex systems.

The Short-Time Drift Propagator (STDP) is an analytical and numerical framework for the solution of the Fokker-Planck equation, specifically optimized for high accuracy over short time increments. Its key innovation is a reformulation of the short-time transition kernel such that the drift term is separated as an external shift, while the convolution kernel is reduced to the pure-diffusion component. This approach enables efficient and rigorously positive-preserving time propagation, robust error control, and superior accuracy compared to standard Gaussian (Euler–Maruyama based) methods, while allowing for systematic improvement to arbitrary order in the time increment. The STDP has been validated for both constant and state-dependent drift and diffusion, including rigorous perturbative expansions, and is connected at a fundamental level with hypocoercivity theory, decay estimates, and operator norm techniques for Fokker-Planck type evolution.

1. Mathematical Formulation of the Fokker-Planck Problem

The Fokker-Planck equation describes the time evolution of a probability density $p(x,t)$ associated with an $n$-dimensional stochastic differential equation (SDE)

$$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$

with the associated Fokker-Planck operator

$$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$

where $D = \frac12 G Q G^\top$ is the diffusion tensor. For small time increments $\Delta t$, the conventional short-time propagator is a displaced multivariate Gaussian whose mean is directly affected by the state-dependent drift, making direct evaluation in general infeasible when $v(x)$ varies with $x$ (Mangthas et al., 2023).

2. Construction of the Short-Time Drift Propagator

The central idea of the STDP is to decouple the drift term from the covariance structure of the short-time transition kernel. The propagator is re-expressed such that the drift enters solely as an overall shift in state space,

$$p(x+v\,\Delta t, t+\Delta t) = \int K(x-x'; \Delta t)\; p(x',t) \, dx',$$

with

$$K(y; \Delta t) = (4\pi\Delta t)^{-\frac{n}{2}} |D|^{-\frac12} \exp\left[-\frac{1}{4\Delta t} y^\top D^{-1} y\right].$$

Hence, the convolution at each step is through a pure-diffusion Gaussian, while the drift is realized by a translation $n$0 (Mangthas et al., 2023). The Kramers-Moyal expansion demonstrates that, for small $n$1, only the second-order (diffusion) term contributes non-trivially inside the kernel.

3. Iterative Numerical Scheme and Extensions

The STDP iterative update consists of:

1. Quadrature evaluation: At each grid point $n$2, compute nodes $n$3, where $n$4 is the matrix square root of $n$5 and $n$6 are Gaussian–Hermite roots.
2. Interpolation: Interpolate $n$7 to the off-grid locations $n$8.
3. Weighted sum: Aggregate with quadrature weights to obtain $n$9.
4. Shift: Update the density by shifting $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$0.

For nonconstant drift or diffusion, local linearization and Trotter product splitting can be utilized, and higher-order accuracy is obtainable through Strang splitting (Mangthas et al., 2023). Adaptive quadrature and local error control are also directly integrable.

4. Systematic Perturbative Improvements

Recent work formalizes the STDP as a perturbative expansion about the frozen-coefficient Gaussian kernel. If the Fokker-Planck equation is

$$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$1

the STDP expansion has the form

$$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$2

where coefficients $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$3 are explicit functions of the derivatives of $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$4 and $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$5 at $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$6 (Kappler, 2024). The expansion is constructed to each order in $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$7, yielding controllable accuracy, and each truncation preserves exact normalization.

This enables closed-form expressions for finite-time Kramers–Moyal coefficients and for the entropy production rate. The error with respect to the exact propagator in $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$8 norm scales as $$\dot{x}(t) = v(x,t) + G(x,t)\, \eta(t), \quad \mathbb{E}[\eta(t)\eta^\top(s)] = Q\delta(t-s),$$9 for truncation order $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$0, which is substantially better than the Gaussian propagator's sublinear scaling (Kappler, 2024).

5. Short-Time Drift Propagator for Linear Drift: Operator Norm and Hypocoercivity

For Fokker-Planck equations with linear drift, i.e., $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$1 with normalized $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$2, the $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$3-propagator norm coincides precisely with that of the drift ODE, $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$4. The short-time expansion of the propagator norm is determined by the hypocoercivity index $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$5 of the drift (the minimal number of commutators needed to detect all directions under the Kalman condition): $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$6 for small $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$7, with $$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$8. For coercive cases ($$\frac{\partial p}{\partial t} = -\sum_{i=1}^n \frac{\partial}{\partial x_i}[v_i(x,t)\,p] + \sum_{i,j=1}^n \frac{\partial^2}{\partial x_i \partial x_j}[D_{ij}(x,t)\,p],$$9) the leading term is $D = \frac12 G Q G^\top$0, while for hypocoercive ($D = \frac12 G Q G^\top$1) the leading correction shifts to higher powers. The operator-theoretic underpinning is a block-diagonalization in Hermite function bases, establishing the Fokker-Planck semigroup as the second quantization of the drift ODE (Arnold et al., 2020).

6. Benchmarks, Advantages, and Limitations

Benchmarks: For 1D and 2D Wiener-drift processes with constant coefficients, the STDP scheme numerically replicates the exact solution within mean absolute error $D = \frac12 G Q G^\top$2 (1D) and $D = \frac12 G Q G^\top$3 (2D) under typical discretization, with convergence rate scaling linearly in $D = \frac12 G Q G^\top$4. In complex multiplicative-noise systems, higher-order STDP outperforms the standard Gaussian propagator by 5–9 orders of magnitude in $D = \frac12 G Q G^\top$5 error for small $D = \frac12 G Q G^\top$6 (Mangthas et al., 2023, Kappler, 2024).

Advantages: STDP factorizes the drift term, leaving a simple convolution kernel, reducing computational complexity, and avoiding spurious negative oscillations or large stencils. The method guarantees positivity of the propagated density and supports trivial parallelization.

Limitations: The core formulation assumes constant or divergence-free drifts; for nonzero divergence ($D = \frac12 G Q G^\top$7), a Jacobian-weight correction is required. For strongly state-dependent coefficients, accurate local linearization or splitting strategies become critical (Mangthas et al., 2023).

7. Applications and Outlook

STDP delivers highly accurate finite-time estimates of propagation, the moments $D = \frac12 G Q G^\top$8, Kramers–Moyal coefficients, and entropy production rates directly, supporting applications in inference from time series and simulation of stochastic trajectories (Kappler, 2024). Extensions to multidimensional, nonconstant, and degenerate-diffusion Fokker-Planck systems are enabled via operator-splitting and block-diagonal techniques. The approach provides a theoretically and numerically robust mechanism for high-precision short time stepping and for understanding decay and regularization mechanisms via operator norm analyses rooted in hypocoercivity theory (Arnold et al., 2020, Mangthas et al., 2023).