Stochastic Keller–Segel Systems

• Stochastic Keller–Segel systems are chemotaxis models that merge classical drift–diffusion dynamics with random forcing to capture variability in aggregation and diffusion.
• They encompass a range of formulations—SPDEs, particle systems, and random media approaches—each addressing noise effects and logistic damping in unique ways.
• Analytical and numerical studies reveal that noise acts as a model-specific regularizer, influencing well-posedness, mean-field limits, and solution behaviors in these complex systems.

Searching arXiv for recent and foundational papers on stochastic Keller–Segel systems. I’m unable to invoke the arXiv search tool in this interface, so I will rely on the arXiv records provided in the data block and cite them directly. A stochastic Keller–Segel system is a chemotaxis model in which the classical Keller–Segel drift–diffusion mechanism is coupled to randomness, either at the SPDE level, through random coefficients or heterogeneous media, or through stochastic particle systems whose empirical measures approximate chemotactic densities. In the literature represented here, the term covers parabolic–elliptic, parabolic–parabolic, and elliptic–parabolic chemotaxis equations; systems coupled to fluid dynamics; McKean–Vlasov and branching particle models; and stochastic homogenization in random media (Huang et al., 2019, Matzavinos et al., 2013). Across these formulations, the central analytical issue remains the same as in deterministic chemotaxis: diffusion and transport compete with aggregation, while stochastic forcing, logistic damping, nonlocal source terms, and fluid coupling alter both the well-posedness mechanism and the relevant notion of solution (Zhang et al., 2023, Chen et al., 2024).

1. Model classes and deterministic core

The deterministic core is the Keller–Segel mechanism

$$\partial_t \rho = \Delta \rho - \chi \nabla\cdot(\rho \nabla c),$$

with $c$ determined either elliptically or parabolically. One stochastic extension studied in dimensions $d=2,3$ is the parabolic–elliptic SPDE

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$

which incorporates both idiosyncratic and common noise (Huang et al., 2019). A distinct 1D formulation on $[0,1]$ uses multiplicative spatial Wiener processes and, in Itô form,

$$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$

with Neumann boundary conditions (Hausenblas et al., 2022).

Other model classes differ by how the chemoattractant is represented and where the noise acts. In bounded 2D domains, an elliptic–parabolic stochastic chemotaxis system is written as

$$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$

with either linear-growth multiplicative noise or nonlinear noise of the form $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$ (Chen et al., 2024). In fluid-coupled settings, one obtains stochastic Keller–Segel–Navier–Stokes systems. A representative 2D fractional-dissipation model on $\mathbb{R}^2$ is

$$\begin{aligned} &\mathrm{d} n +u\cdot \nabla n \,\mathrm{d}t = \Delta n \,\mathrm{d}t- \mathrm{div}(n \nabla c)\,\mathrm{d}t+(n-n^2)\,\mathrm{d}t,\ &\mathrm{d} c+ u\cdot \nabla c\,\mathrm{d}t = \Delta c\,\mathrm{d}t-n c\,\mathrm{d}t,\ &\mathrm{d} u+ (u\cdot \nabla) u\,\mathrm{d}t+\nabla P \,\mathrm{d}t = - (-\Delta)^\alpha u\,\mathrm{d}t + n\nabla \phi\,\mathrm{d}t+ f(t,u)\,\mathrm{d} W_t, \end{aligned}$$

with $c$0 and logistic source $c$1 (Zhang et al., 2023). A self-consistent 2D parabolic–elliptic fluid system instead couples

$c$2

on $c$3 (Kong et al., 16 May 2026).

The phrase also covers stochastic representations that are not SPDE perturbations in the narrow sense. One line of work derives Keller–Segel equations from interacting particle systems with Brownian noise, common noise, or branching and competition (García et al., 2017, Cavallazzi et al., 23 Dec 2025). Another studies Keller–Segel equations in random heterogeneous media, where the randomness is encoded in stationary ergodic coefficients $c$4 and $c$5, and the limit is obtained by stochastic homogenization rather than by a temporal noise term (Matzavinos et al., 2013).

2. Sources of stochasticity and notions of solution

The stochastic forcing mechanisms vary sharply across the literature. In the parabolic–elliptic SPDE with common noise, idiosyncratic noise contributes only to the diffusion matrix, while common environmental noise survives in the mean-field limit as the transport-type term

$c$6

and the limiting density is the conditional law of a McKean–Vlasov SDE given the common noise filtration (Huang et al., 2019). In the 1D Patlak–Keller–Segel model of Hausenblas–Mukherjee–Tran, randomness enters through time-homogeneous spatial Wiener processes $c$7, first in Stratonovich form and then in Itô form with correction $c$8 (Hausenblas et al., 2020, Hausenblas et al., 2022). In the elliptic–parabolic bounded-domain problem, the noise may be globally Lipschitz with linear growth or genuinely nonlinear through $c$9 (Chen et al., 2024). In the fractional KS–SNS model, the noise acts only in the fluid equation and is multiplicative, globally Lipschitz, of linear growth, with additional vorticity control assumptions (Zhang et al., 2023). In the recent bounded-domain model with nonlocal Fisher–KPP source, the noise is super-linear in the sense

$d=2,3$0

which forces a different compactness framework (Li et al., 10 Jun 2026).

The corresponding solution concepts are equally heterogeneous. The common-noise parabolic–elliptic equation is studied as an SPDE whose weak solution satisfies a duality formula against $d=2,3$1 test functions, and whose McKean–Vlasov counterpart is interpreted via conditional densities (Huang et al., 2019). The 1D spatial-noise model first yields martingale solutions and later, after pathwise uniqueness, probabilistically strong solutions in

$d=2,3$2

(Hausenblas et al., 2022). The fractional KS–SNS paper distinguishes analytically weak solutions, martingale weak solutions, and probabilistically strong solutions, and proves pathwise uniqueness for martingale solutions (Zhang et al., 2023). The bounded 2D elliptic–parabolic chemotaxis paper works with local and global mild solutions in $d=2,3$3 and uses stopping times to encode maximal existence (Chen et al., 2024). In the parabolic–parabolic 2D particle interpretation, the limiting object is formulated as a nonlinear martingale problem on path space, because the drift depends on the whole past of the law (Fournier et al., 2022).

A recurrent point of terminology is that “stochastic Keller–Segel system” may mean at least three different things: an SPDE for the density, a stochastic particle system approximating a deterministic or random PDE, or a random-media problem with deterministic time evolution for each realization. The literature here uses all three meanings, and the distinction is structural rather than merely linguistic (Matzavinos et al., 2013, Talay et al., 2017).

3. Well-posedness, thresholds, and global regimes

Global well-posedness is known in several regimes, but the stabilizing mechanisms are model-dependent. In one dimension, the stochastic Patlak–Keller–Segel system with multiplicative spatial Wiener processes admits martingale solutions on finite intervals, and subsequent work proves additional regularity, pathwise uniqueness, and the existence of a unique strong solution (Hausenblas et al., 2020, Hausenblas et al., 2022). For the parabolic–elliptic common-noise SPDE in $d=2,3$4, there is existence and uniqueness for small $d=2,3$5-data in spaces such as

$d=2,3$6

together with a quantitative mean-field limit from the particle system (Huang et al., 2019).

In bounded 2D domains, global existence can be enforced either by logistic damping or by specific nonlinear noise structures. For the elliptic–parabolic system with linear-growth multiplicative noise, the hypotheses $d=2,3$7, $d=2,3$8, and sufficiently strong logistic damping yield a unique global nonnegative mild solution (Chen et al., 2024). The same paper shows that for the nonlinear-noise model

$d=2,3$9

one also has unique global nonnegative mild solutions under explicit exponent constraints. The paper interprets this as a regularization-by-noise mechanism for a specific superlinear multiplicative structure rather than as a generic feature of stochastic chemotaxis (Chen et al., 2024).

For fluid-coupled problems, two distinct global regimes are currently available. In the 2D stochastic KS–SNS system with fractional fluid dissipation and logistic source, one has a unique global pathwise weak solution for $$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$0, with

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$1

almost surely (Zhang et al., 2023). In the self-consistent 2D parabolic–elliptic KS–Navier–Stokes system on $$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$2, global mild solutions exist under the subcritical mass condition

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$3

and satisfy a stochastic free-energy estimate based on

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$4

(Kong et al., 16 May 2026).

Recent work also shows global nonnegative martingale solutions for a bounded-domain stochastic Keller–Segel equation with nonlocal Fisher–KPP source and super-linear multiplicative noise,

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$5

without any smallness assumption on the initial data, provided the nonlocal source is dominant and

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$6

(Li et al., 10 Jun 2026). On the deterministic side of the logistic Keller–Segel equation,

$$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$7

global bounded mild solutions exist when $$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$8; in dimension $$d \rho_t =\frac{1}{2}\sum_{i,j=1}^d D_{ij}\Big(\rho_t \sum_{k=1}^{d'}(\nu^{ik}_t\nu^{jk}_t+ \sigma^{ik}_t\sigma^{jk}_t)\Big)\,dt -\chi\nabla\cdot(\nabla c_t\rho_t)\,dt -\sum_{i=1}^{d}D_i\Big( \rho_t \sum_{k=1}^{d'}\sigma^{ik}_t\,d W_t^k\Big), \quad -\Delta c_t + c_t = \rho_t,$$9, any positive logistic damping suppresses blow-up (Cavallazzi et al., 23 Dec 2025).

A common misconception is that stochastic forcing itself is the universal regularizer in chemotaxis. The available results do not support that uniform claim. One line of work explicitly states that “noise alone” is not the key stabilizer in the stochastic KS–SNS setting, where fractional dissipation and logistic damping are decisive (Zhang et al., 2023). Another proves only small-data well-posedness for common-noise parabolic–elliptic KS and does not claim any shift of the deterministic blow-up threshold (Huang et al., 2019). By contrast, the bounded 2D elliptic–parabolic model shows that a carefully chosen nonlinear multiplicative noise can play a role similar to a logistic term (Chen et al., 2024). The stabilizing mechanism is therefore model-specific.

4. Particle systems, mean-field limits, and propagation of chaos

A major branch of the theory studies stochastic Keller–Segel systems as interacting particle models. For the 2D parabolic–elliptic equation, a singular $[0,1]$0-particle diffusion

$[0,1]$1

is constructed directly below the critical mass. For $[0,1]$2 and

$[0,1]$3

there exists a unique non-explosive solution in distribution starting from any configuration in the set $[0,1]$4 of configurations with at most one collision pair, and the process admits the symmetric $[0,1]$5-finite invariant measure

$[0,1]$6

(Cattiaux et al., 2016). A quantitative microscopic derivation with cutoff $[0,1]$7, $[0,1]$8, proves that the maximal particle-wise distance between the interacting system and the mean-field system stays below $[0,1]$9 with probability at least $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$0 on fixed finite intervals, and deduces propagation of chaos toward the 2D Keller–Segel PDE in the full subcritical regime $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$1 (García et al., 2017).

For the 1D parabolic–parabolic problem, a new McKean–Vlasov interpretation uses a time-memory drift,

$$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$2

where $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$3 is the density of $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$4, $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$5 encodes the initial chemical field, and $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$6 is the gradient of the heat kernel. This formulation yields global well-posedness in one dimension for all parameters and connects the time marginals of the McKean–Vlasov process to the mild solution of the parabolic–parabolic Keller–Segel system (Talay et al., 2017). In two dimensions, the corresponding doubly parabolic particle system is non-Markovian: each particle interacts with the entire past of every other particle through

$$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$7

For sufficiently small $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$8, the paper proves existence of the $$du = \bigl(r_u A u - \chi \operatorname{div}(u\nabla v) + \gamma\,u\bigr)\,dt + u\,dW_1,\qquad dv = \bigl(r_v A v - a v + B u\bigr)\,dt + v\,dW_2,$$9-particle system, tightness of empirical measures, and that every weak limit solves a nonlinear martingale problem whose time marginals satisfy the parabolic–parabolic Keller–Segel equation in weak sense (Fournier et al., 2022).

Mean-field limits have also been established for stochastic KS equations with common noise. Starting from the particle system

$$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$0

one obtains a McKean–Vlasov limit with common noise and a propagation-of-chaos estimate that is conditional on the common noise filtration (Huang et al., 2019). For degenerate diffusion, a moderately interacting stochastic particle system with extra Brownian regularization and logarithmic cutoffs rigorously converges, in the propagation-of-chaos sense, to a degenerate parabolic–elliptic Keller–Segel equation with porous-medium diffusion in the subcritical regime (Chen et al., 2023).

A particularly recent development adds demographic branching to the microscopic model. The branching moderately interacting particle system with attractive Coulomb-type interaction and local competition converges to the logistic Keller–Segel PDE with rate $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$1, and the birth–death mechanism is explicitly interpreted as the microscopic origin of the macroscopic term $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$2 (Cavallazzi et al., 23 Dec 2025).

5. Analytical structures and proof techniques

The modern theory relies on several distinct analytical templates. For SPDE well-posedness in bounded domains, localization and $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$3-Itô formulas are fundamental. In the elliptic–parabolic 2D chemotaxis system, the drift is cut off by a smooth function $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$4, local mild solutions are produced by a Banach fixed point argument for the truncated equation, and global existence follows from Yosida-regularized $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$5-Itô formulas, Gagliardo–Nirenberg estimates, and generalized Gronwall–Bellman inequalities (Chen et al., 2024). In the recent bounded-domain model with nonlocal Fisher–KPP source and super-linear noise, the Galerkin approximation must itself be further cut off in finite $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$6 norm, positivity of the approximants must be proved to exploit the Fisher–KPP structure, and compactness is obtained through a generalized tightness criterion combined with Jakubowski’s version of the Skorokhod theorem (Li et al., 10 Jun 2026).

For the stochastic KS–SNS system with fractional dissipation, the main technical innovation is a three-layer approximation: mollification $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$7, Fourier truncation $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$8, and high-norm cut-off $$du = \bigl( A u - \nabla\cdot(\chi u \nabla v) + g(u) \bigr)\,dt + \sum_{k=1}^\infty \sigma_k(u)\,dW^k_t,\qquad 0 = A v + u - v,$$9. This produces globally Lipschitz coefficients in $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$0 and allows the use of infinite-dimensional SDE theory. The convergence to the original system then proceeds through entropy–energy inequalities, stochastic compactness, and a final Littlewood–Paley/Besov uniqueness argument. The uniqueness proof uses bilinear estimates such as

$\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$1

together with vorticity estimates and BDG control of the noise term (Zhang et al., 2023).

For stochastic particle and McKean–Vlasov formulations, singular-kernel control requires a different toolkit. The 2D parabolic–elliptic particle system is analyzed through Dirichlet forms, $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$2-symmetric diffusions, and comparison with squared Bessel processes that detect multi-particle collisions and total collapse (Cattiaux et al., 2016). The 1D parabolic–parabolic McKean–Vlasov theory uses a sharp density estimate for one-dimensional diffusions with bounded drift,

$\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$3

a Picard scheme for the time-memory drift, and a singular Gronwall lemma adapted to kernels of order $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$4 (Talay et al., 2017). The 2D doubly parabolic particle system introduces a “Markovianization” argument: the path-dependent interaction is controlled by a Coulomb-type interaction depending only on current pairwise distances, which allows the non-Markovian singular drift to be estimated through pair-separation functionals (Fournier et al., 2022).

A separate analytical tradition concerns random media. In stochastic homogenization of Keller–Segel, the random coefficients $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$5 and $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$6 are treated by stochastic two-scale convergence, stochastic gradients $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$7, and corrector problems in the ergodic probability space. The limit is a deterministic homogenized chemotaxis system with effective tensors $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$8 and $\sum_i b_i\|u\|_r^\gamma u\,dW_t^i$9 defined through random cell problems rather than through stochastic integration (Matzavinos et al., 2013).

6. Computation, random media, and current directions

Numerical work on stochastic Keller–Segel systems has recently moved toward particle–field hybrids. In three dimensions, an efficient stochastic interacting particle-field algorithm approximates the fully parabolic Keller–Segel system by empirical particle measures for $\mathbb{R}^2$0 and a spectral representation for $\mathbb{R}^2$1. The crucial numerical idea is to replace the parabolic heat-kernel history formula by an implicit Euler step for the chemoattractant equation,

$\mathbb{R}^2$2

which leads to a one-step recursion through the Yukawa kernel

$\mathbb{R}^2$3

This avoids history dependence and allows the method to study aggregation, multimodal cluster merging, and finite-time singularity indicators in 3D (Wang et al., 2023). The same paper reports approximately first-order time convergence, $\mathbb{R}^2$4, in an early-time test and uses the scaling of $\mathbb{R}^2$5 with the number of Fourier modes to distinguish near-singular from smooth regimes (Wang et al., 2023).

Random-media formulations point in a different direction. In stochastic homogenization, the microscopic coefficients are stationary ergodic random fields, and the large-scale limit is again a deterministic Keller–Segel system but with effective diffusion and chemotaxis tensors computed by stochastic cell problems. The periodization procedure over large cubes provides an approximation of these homogenized coefficients and converges almost surely as the period size tends to infinity (Matzavinos et al., 2013). This line of work shows that “stochastic Keller–Segel” does not always mean temporal noise; it may also mean chemotaxis in a spatially random environment.

Several directions are now explicit in the literature. For stochastic KS–SNS, natural extensions include noise acting directly on the chemotaxis equations, exponents $\mathbb{R}^2$6, higher dimensions, and more realistic geometries (Zhang et al., 2023). For particle approximations, the main open problems are removing cutoffs, understanding critical and supercritical mass regimes, and extending propagation-of-chaos results to singular interactions with branching or memory (Cavallazzi et al., 23 Dec 2025, Fournier et al., 2022). For self-consistent fluid coupling, the subcritical mass theory on $\mathbb{R}^2$7 leaves open the stochastic behavior at and above $\mathbb{R}^2$8, as well as long-time statistical questions such as invariant measures and ergodicity (Kong et al., 16 May 2026). The current picture is therefore technically rich but not unified: stochastic Keller–Segel systems form a family of related models whose qualitative behavior depends decisively on how randomness enters, on whether aggregation is counterbalanced by damping or fluid dissipation, and on whether the stochastic description is macroscopic, microscopic, or homogenized.