A fast stochastic interacting particle-field method for 3D parabolic parabolic Chemotaxis systems: numerical algorithms and error analysis

TL;DR

This work addresses the computational bottlenecks of simulating the 3D parabolic-parabolic Keller-Segel system by introducing the SIPF-PIC method, a hybrid stochastic particle-field framework accelerated with FFT-based spectral solvers. The authors derive a per-step complexity of $\mathcal{O}(P + H^3 \log H)$ and establish a rigorous error bound $\mathcal{O}(H^{-16/13} + P^{-1/2}H^{4/13})$ for the discretization, leveraging a Brownian-coupling analysis. Numerical experiments validate convergence rates, demonstrate substantial speedups over the original SIPF approach, and reveal the method’s ability to resolve complex blow-up phenomena, including ring-like and tetrahedral patterns, in 3D. The results provide a scalable tool for exploring multiscale chemotaxis dynamics with rigorous theoretical guarantees and practical computational efficiency.

Abstract

In this paper, we develop a novel numerical framework, the stochastic interacting particle-field method with particle-in-cell acceleration (SIPF-PIC), for the efficient simulation of the three-dimensional (3D) parabolic-parabolic Keller-Segel (KS) systems. The SIPF-PIC method integrates Lagrangian particle dynamics with spectral field solvers, by leveraging localized particle-grid interpolations and fast Fourier transform (FFT) techniques. For $P$ particles and $H$ Fourier modes per spatial dimension, the SIPF-PIC method achieves a computational complexity of $\mathcal{O}(P + H^3 \log H)$ per time step, a significant improvement over the original SIPF method (proposed in \cite{SIPF1}), which has a complexity of $\mathcal{O}(PH^3)$, while preserving numerical accuracy. Moreover, we establish a rigorous error analysis, proving that the discretization errors are of order $\mathcal{O}(H^{-16/13}+P^{-1/2}H^{4/13})$. Finally, we present numerical experiments to validate the theoretical convergence rates and demonstrate the computational efficiency of our new method. Notably, these experiments also show that the method captures complex blowup dynamics beyond single-point collapse, including ring-type singularities, where mass dynamically concentrates into evolving annular structures.

Figures (9)

• Figure 1: Extended vertices for fourth-order interpolation
• Figure 1: Convergence characteristics with dual parameter analysis.
• Figure 2: Convergence with increasing grid resolution $H$ under a fixed particle count ($P = 10^4$) and truncation frequency $H_0 = 16$.
• Figure 3: Two-dimensional particle distributions at $T = 2.0$ for different initial masses. (a)-(d): Subcritical to near-critical regimes ($M_0 = 22.0$--$25.0$); (e)-(f): Supercritical regimes ($M_0 = 26.0$--$27.0$).
• Figure 4: Three-dimensional particle distributions at $T = 1.0$ for different initial masses. (a)-(d): Subcritical to near-critical regimes ($M_0 = 45.0$--$48.0$); (e)-(f): Supercritical regimes ($M_0 = 49.0$--$50.0$).

Theorems & Definitions (33)

• Definition 3.1
• Definition 3.2
• Definition 3.3
• Lemma 3.4
• Proof 1
• Lemma 3.5
• Proof 2
• Theorem 3.6
• Proof 3
• Lemma 3.7