Convergence Analysis of a Stochastic Interacting Particle-Field Algorithm for 3D Parabolic-Parabolic Keller-Segel Systems

Technical Summary: Convergence Analysis of a Stochastic Interacting Particle-Field Algorithm for 3D Parabolic-Parabolic Keller-Segel Systems

Problem Statement
The paper addresses the numerical simulation of the three-dimensional (3D) fully parabolic Keller-Segel (KS) system, a model describing chemotaxis where organisms (density $\rho$) move in response to chemical gradients (concentration $c$). The system is characterized by the equations:
$$\rho_t = \nabla \cdot (\mu \nabla \rho - \chi \rho \nabla c), \quad \epsilon c_t = \Delta c - \lambda^2 c + \rho$$
The primary challenges in solving this system numerically include the potential for finite-time blow-up (singularity formation) and the high computational cost associated with resolving the nonlinear coupling between the particle density and the chemical field, particularly in 3D. While mesh-based methods (finite difference, finite element) and previous particle-based approaches exist, they often struggle with stability, convergence guarantees, or the computational complexity of pairwise particle interactions ($O(P^2)$).

Methodology
The authors propose and analyze the Stochastic Interacting Particle-Field method with Random Batch approximation (SIPF-r). This hybrid algorithm combines particle-based representations for the density with a spectral field method for the chemical concentration:

1. Particle-Field Coupling: The density $\rho$ is approximated by an empirical measure of $P$ particles, while the chemical concentration $c$ is represented by a truncated Fourier series with $H$ modes. The chemical field is updated using an implicit Euler scheme and spectral methods, while particle positions are updated via an Euler-Maruyama scheme driven by the gradient of the chemical field.
2. Random Batch Method (RBM): To reduce the computational complexity of particle interactions from $O(P^2)$ to $O(PR)$ (where $R \ll P$ is the batch size), the algorithm employs the RBM. At each time step, particles are grouped into small random batches. Interactions are computed only within these batches, rendering the particle trajectories independent and identically distributed (i.i.d.) at each step. This bypasses the need to explicitly handle the propagation of chaos.
3. Discretization: The method utilizes a large spatial domain with periodic boundary conditions to emulate an infinite domain, assuming the singularity remains localized away from boundaries. The chemical field is computed using Green's functions derived in the Fourier domain.

Key Contributions and Theoretical Results
The paper provides a rigorous convergence analysis for the SIPF-r method, establishing that the empirical measure of the particle system converges to the exact solution of the limiting McKean-Vlasov process.

• Convergence Rates: Under mild regularity assumptions on the exact solution and boundedness of numerical approximations, the authors prove that the 1-Wasserstein distance between the numerical density $\tilde{\rho}$ and the exact density $\rho$ scales as:
  $$W_1(\tilde{\rho}_t, \rho_t) = O\left( \frac{1}{H^2} + \sqrt{\frac{\delta t}{P}} + \delta t \right)$$
  where $\delta t$ is the time step, $H$ is the number of Fourier modes, and $P$ is the number of particles.
• Error Decomposition: The analysis decomposes the error into components arising from:
  • Fourier Truncation: Decays as $O(1/H^2)$ for density and $O(1/H)$ for concentration gradients.
  • Time Discretization: First-order convergence $O(\delta t)$.
  • Statistical Sampling: The error due to the random batch approximation scales as $O(\delta t / \sqrt{R})$, which is shown to be subordinate to the leading-order discretization errors.
• Stability: The authors demonstrate that the SIPF-r method possesses inherent numerical stability. The spectral discretization and implicit time-stepping ensure that the numerical chemical concentration and its gradient remain uniformly bounded for any fixed $H$, even when the exact solution approaches blow-up. This allows the method to capture the growth trend of singularities without numerical breakdown.

Numerical Results
The theoretical findings are validated through extensive numerical experiments:

• Accuracy: Comparisons with a high-resolution radial finite difference method (FDM) benchmark show a mean relative error of approximately 0.055, confirming high accuracy.
• Convergence Verification: Numerical tests confirm the theoretical convergence rates. The error with respect to the time step $\delta t$ exhibits first-order convergence ($O(\delta t^{1.023})$), and the error with respect to the batch size $R$ follows the predicted $O(R^{-0.5})$ scaling.
• Blow-up Detection: The method successfully detects finite-time blow-up and critical mass thresholds in 3D. Experiments show that the algorithm can identify intense focusing and potential singularities for initial masses near critical thresholds (e.g., $M_0 \approx 73$) using moderate parameters ($P=10^4, H=24$). Notably, the method remains robust in detecting blow-up even under discretization parameters significantly more relaxed than those required for the convergence proof.

Significance
The paper claims that the SIPF-r method offers a computationally efficient and theoretically grounded approach for simulating 3D fully parabolic Keller-Segel systems. By integrating the Random Batch Method with a particle-field framework, the algorithm achieves a substantial reduction in computational cost while maintaining rigorous convergence guarantees. The work is significant for its ability to handle the complex dynamics of chemotaxis, including the detection of finite-time singularities, without requiring the stringent constraints often necessary for traditional mesh-based or full-interaction particle methods. The authors note that the method serves as a robust diagnostic tool for studying intense focusing phenomena in three dimensions.