Numerical stationary states for nonlocal Fokker-Planck equations via fixed points of consistency maps

Technical Summary: Numerical Stationary States for Nonlocal Fokker-Planck Equations via Fixed Points of Consistency Maps

Problem Statement
The paper addresses the numerical computation of stationary states for a class of nonlocal Fokker-Planck-type equations of the form:
$$\partial_t u = \nabla \cdot (u\nabla(H'[u] + \Phi[u](x) + V(x))),$$
where $u(x,t)$ is a density, $H(u)$ represents internal energy, $V(x)$ is an external potential, and $\Phi[u]$ encodes nonlocal, typically nonlinear, interactions. These equations arise in diverse fields such as fluid dynamics, phase separation, cell adhesion, and biological aggregation.

The primary challenge lies in the long-time dynamics of these systems. Standard time-dependent solvers (e.g., finite-volume schemes) are limited to detecting dynamically stable states. They often fail to identify unstable stationary states or complex bifurcation branches, particularly in systems exhibiting metastability or multiple time scales. Furthermore, existing numerical continuation software often struggles with the specific structural requirements of general nonlinear diffusion and nonlocal terms found in this class of equations.

Methodology
The authors propose a fixed-point-based numerical framework that avoids direct discretization of the differential operators. Instead, the stationary problem is reformulated as a nonlinear fixed-point equation $T u - u = 0$, where $T$ is a map derived from the original PDE and its nonlocal interaction terms.

1. Fixed-Point Reformulation:
   By setting $\partial_t u = 0$ and utilizing mass conservation and positivity, the stationary state $u^*$ is shown to satisfy a consistency condition on its support: $H'[u] + \Phi[u] + V = c_0$. Assuming $H'$ is invertible, the stationary solution is expressed as a fixed point of the map:
   $$T u := (H')^{-1}[c_0(u) - \Phi[u](x) - V(x)].$$
   For linear diffusion ($H[u] = \sigma u(\log u - 1)$), this inverse is an exponential function.

2. Numerical Solution Strategy:
   The authors solve the discretized system $F_n(u) = T_n(u) - u = 0$ using a Jacobian-free Newton-Krylov (JFNK) method.

   • Outer Loop: A Newton iteration updates the solution.
   • Inner Loop: A Krylov-subspace method (specifically GMRES) solves the linear system $J \delta u = -F(u)$.
   • Matrix-Free Approach: The method does not explicitly construct the Jacobian matrix. Instead, it computes the action of the Fréchet derivative $DT[u](\phi)$ on a vector $\phi$.
   • Derivative Approximation: The authors compare two approaches for the Fréchet derivative:
     • Analytic: Deriving the exact closed-form expression for the derivative of the map $T$.
     • Finite-Difference: Approximating the derivative using central differences.

3. Discretization:
   The nonlocal terms (convolutions) are handled via Fast Fourier Transforms (FFT) on periodic domains or appropriate quadrature rules. This shifts the accuracy burden from differentiation stencils to the precision of the convolution and integration rules.

Key Contributions

• Framework Development: The authors develop and test a robust fixed-point framework for three representative models: the McKean-Vlasov equation, a Cucker-Smale flocking model, and a neural Fokker-Planck equation.
• Analytic vs. Discrete Derivatives: The paper provides a direct comparison between using the analytic Fréchet derivative and a finite-difference approximation. It demonstrates that the analytic approach offers guaranteed accuracy and computational efficiency, though the finite-difference approach remains viable when analytic forms are intractable.
• Bifurcation Analysis: The method successfully reproduces known bifurcation diagrams and detects both stable and unstable stationary states, including supercritical and subcritical branches.
• Discovery of New Behavior: The framework reveals complex bifurcation structures not previously observed in these specific problems, particularly in cases involving multi-mode kernels and two-dimensional domains.

Results
The authors validate their approach against existing analytical results and known bifurcation theories:

• Accuracy: For the McKean-Vlasov equation with a Kuramoto-type kernel, the numerical solutions converge to the analytical reference solution with an error decay rate of $O(N_n^{-4})$, where $N_n$ is the number of grid points. The critical bifurcation thresholds are recovered with high precision ($\sim 10^{-5}$).
• Multi-Mode Kernels: Using bimodal interaction kernels, the method identifies multiple solution branches, including subcritical bifurcations and secondary branches emerging from the primary bifurcation points. In cases where the linearized operator has a two-dimensional kernel (degenerate bifurcation), the method detects transcritical-like behavior that standard local theory struggles to characterize.
• Stability: By combining the fixed-point solver with a time-dependent solver for perturbation analysis, the authors map the stability of various branches. They observe that while supercritical branches are generally stable, subcritical branches are unstable, and complex interactions (e.g., intersections of branches with distinct profiles) occur.
• Generalizability: The approach is successfully extended to:
  • 2D Domains: Demonstrating applicability on the torus $\mathbb{T}^2$.
  • Different Domains: Handling unbounded domains (Cucker-Smale) and semi-infinite domains with no-flux boundaries (Neural Fokker-Planck).
  • Various Kernels: Testing top-hat, triangular, and attractive-repulsive kernels, revealing rich solution structures including secondary bifurcations.

Significance
The paper claims that this fixed-point approach offers an excellent tool for exploring the stationary solutions of nonlocal Fokker-Planck equations. Its primary significance lies in its agnosticism to dynamical stability, allowing for the detection of unstable states and complex bifurcation branches that time-evolution solvers miss. By leveraging the structural properties of the PDE to recast the problem as a fixed-point map, the method avoids the difficulties associated with discretizing differential operators directly. The authors note that while the current work focuses on linear diffusion, the framework is designed to be extensible to nonlinear diffusion (e.g., porous medium) and systems of interacting populations, provided the necessary structural conditions for the fixed-point formulation are met.