1 Introduction

Chemotaxis partial differential equations (PDEs) were introduced by Keller and Segel (KS [16]) to describe the aggregation of the slime mold amoeba Dictyostelium discoideum due to an attractive chemical substance. A related random walk model was known earlier, developed by Patlak [26]. For an analysis of basic taxis behaviors (such as aggregation, blowup, and collapse) based on reinforced random walks, we refer to [31]. In this work, we consider the parabolic-parabolic (fully parabolic) KS system of the following form:

$$\begin{aligned} \rho _{t}&= \nabla \cdot (\mu \, \nabla \rho - \chi \, \rho \, \nabla c ), \nonumber \\ \epsilon c_t&= \varDelta \, c - k^2 \, c + \rho , \end{aligned}$$
(1)

where \(\chi , \mu \) (\(\epsilon , k\)) are positive (non-negative) constants. The model is called elliptic if \(\epsilon =0\) (when c evolves rapidly to a local equilibrium), and parabolic if \(\epsilon >0\). Here, \(\rho \) is the density of active particles (such as bacteria), and c is the concentration of chemo-attractant (e.g. food). For a more detailed discussion, please refer to Sect. 2.

The KS systems (1) have been extensively studied for several decades, with various cases and dimensions explored. It is well known that the system (1) in \({\mathbb {R}}^2\) converges to Dirac-delta type function in finite time given sufficiently large initial mass in the parabolic-elliptic case (\(k=0\) and \(\epsilon =0\)) [27] and the fully parabolic case [13]. In addition to the \(\delta \) type singularity, \(C(T-t+|x|^2)^{-1}\) is shown to be an alternative profile for self-similar blowup [9, 29] in 3D parabolic elliptic case. Given the smallness assumption of the initial data, [1, 18, 25, 33] established global well-posedness results for 3D fully parabolic KS in various function spaces. To date, the complete characterization (criteria and profile) of the finite time 3D singularity of (1) with \(\epsilon > 0\) remains largely open. See also [14] for a related model that does not blow up and still exhibits spiky solutions by assuming a saturation concentration for the bacteria.

Many notable numerical methods have been developed for KS systems. Chertock et al. [5] developed a finite-volume method on a class of chemotaxis and haptotaxis models for accurate and efficient simulations. Shen et al. [28] proposed an energy dissipation and bound preserving scheme that is not restricted to specific spatial discretization methods. The bound preserving property is achieved through modifications of the system. Chen et al. [3] developed a fully-discrete finite element method (FEM) scheme for the 2D parabolic-elliptic KS, following the approach of Shen et al. [28]. They showed that the proposed scheme will blow up in a finite time, under assumptions similar to those in the continuous blow-up scenarios. In the classic setting, Liu and Wang [20] reformulated the equation using the Le Châtelier Principle to attain a positive-preserving scheme. See [19, 38] among others on discontinuous Galerkin methods. Besides the aforementioned numerical methods for 2D KS, we refer to [7, 32] on mesh-based (finite element, finite volume/difference with difference of potential domain decomposition) methods for 3D fully parabolic KS models.

As an alternative to the Eulerian discretization methods mentioned above, there have been steady advancements in the Lagrangian formulations for the KS system (1) and related equations. The Lagrangian framework approximates \(\rho \) in (1) as the density of an N-interacting particle system as N tends to infinity. Stevens [30] derived a convergent N-particle system for the fully parabolic KS. Craig and Bertozzi [6] proved the convergence of a blob method for a related aggregation equation. Liu et al. [21, 22] developed a random particle blob method with a mollified kernel for the parabolic-elliptic KS, and proved its convergence when the limiting (macroscopic mean field) equation admits a global weak solution. See also [2] for an extension to the fully parabolic KS. As noted in [24], the success of this line of analysis relies on the limiting nonlinear mean field equation, rather than the underlying many-particle Markov process.

On numerically implemented particle methods for KS systems, Havskovéč- Ševčovič [10, 11] developed a convergent regularized particle system for the 2D parabolic-elliptic KS and in the presence of a passive flow [17]. The authors of this paper conducted a deep learning study on chemotaxis aggregation in 3D laminar and chaotic flows, based on a kernel regularized particle method for parabolic-elliptic KS systems [35]. For 2D fully parabolic KS, [8] proposed a particle-field method and [34] investigated a fully particle approach. Though both methods may be generalized to higher dimensions, the former suffers from numerical instability when capturing blowup behavior and the latter is inefficient for long-time simulation due to memory costs. See Sect. 2.2 for an in-depth discussion and comparison with our proposed method.

To the best of our knowledge, an efficient particle-based (mesh-free) method capable of characterizing blow-up behaviors of 3D fully parabolic KS, especially the critical mass, is unknown. In this paper, we propose a novel stochastic interacting particle-field (SIPF) algorithm to serve this purpose.

Our method takes into account the coupled stochastic particle evolution (density \(\rho \)) and the accompanying field (concentration c) in the system and allows for a self-adaptive simulation of focusing and potentially singular behavior. In the SIPF algorithm, we represent the active particle density \(\rho \) using empirical particles, while the concentration field c is discretized using a spectral method instead of a finite difference method [8]. This is possible since the field c is smoother than density \(\rho \) and does not require interpolation of \(\rho \) to grid points as [8]. We demonstrate the effectiveness of our method through numerical experiments in 3D space, which, to the best of our knowledge, have not been systematically computed and benchmarked before.

It is worth noting that pseudo-spectral methods were employed to compute the nearly singular solutions of the 3D Euler equations [15]. Subsequently, the finite-time blowup of 3D axisymmetric Euler equations was computed using the adaptive moving mesh method [23]. These methods represent the cutting edge in the computation of nearly singular solutions of 3D Euler equations. Nevertheless, we also point out that the implementation of pseudo-spectral methods for 3D problems demands substantial computational resources, while the adaptive moving mesh method requires sophisticated design and advanced programming skills.

It is important to recognize that the Lagrangian algorithms used in the computation of parabolic-elliptic KS systems, such as the one proposed in [10], cannot be directly generalized to the fully parabolic KS. These algorithms rely on the assumption that the field c at time t can be accessed through particle density \(\rho \) at the same instant. Hence only the local update of the particle density is required. However, a direct generalization to the fully parabolic KS (e.g. [34]) will require historical particle density \(\rho \) from the starting time of the algorithm. An example and related convergence analyses can be found in [2]. Nonetheless, from a computational perspective, the volume of such historical data grows over time, posing a costly burden on memory and computational resources. In contrast, our SIPF algorithm computes particle and field values only once per time step, without relying on a long history. Therefore, our computational cost does not increase over time.

The main objective of this paper is to propose a novel stochastic interacting particle-field algorithm for the fully parabolic KS system. While we provide stability analysis and numerical verification of the convergence of the SIPF algorithm, a detailed theoretical analysis will be left as a future work.

The rest of the paper is organized as follows. In Sect. 2, we briefly review the blow-up behavior in the fully parabolic KS models under critical mass conditions and the Lagrangian formulations for computation. Section 3 outlines our proposed SIPF algorithms for solving the fully parabolic KS system. We simplify a theoretically equivalent method with history-dependent parabolic kernel functions, which is computationally undesirable, into efficient recursions. In Sect. 4, we present numerical results to demonstrate the performance of our method for both radial and multi-modal initial data. The SIFP results are in agreement with those from a fully resolved finite difference method in the radially symmetric case. Concluding remarks are in Sect. 5.

2 Parabolic-Parabolic KS System

In this section, we provide some theoretical analyses of singular behaviors and related computational methods for KS models in both parabolic-elliptic cases and parabolic-parabolic (fully parabolic) cases. We begin by revisiting the KS model:

$$\begin{aligned} \rho _{t}&= \nabla \cdot (\mu \, \nabla \rho - \chi \, \rho \, \nabla c ), \end{aligned}$$
(2)
$$\begin{aligned} \epsilon c_t&= \varDelta \, c - k^2 \, c + \rho , \end{aligned}$$
(3)
$$\begin{aligned} x\in&\varOmega \subseteq {\mathbb {R}}^d,\quad t\in [0,T]. \end{aligned}$$
(4)

The first Eq. (2) of \(\rho \) models the evolution of the density of active particles, such as bacteria. The bacteria diffuse with mobility \(\mu \) and drift in the direction of \(\nabla c\) with a velocity \(\chi \nabla c\), where \(\chi \) is called chemo-sensitivity. The second Eq. (3) of c models the evolution of the concentration of chemo-attractant, such as food. The increment in c is proportion to \(\rho \), which indicates the aggregation or attraction between active particles. An additional important physical parameter in the model is \(\epsilon \) in Eq.(3), which represents the timescale of the chemotaxis. When \(\epsilon \not =0\), the system is referred to as a parabolic-parabolic KS system. For \(\epsilon =0\) the system is reduced to the parabolic-elliptic case, which assumes that the chemical attractant released by the active particle instantaneously reaches equilibrium.

2.1 From Critical Collapse to Coexistence of Blow-up and Global Smooth Solutions

Well-known KS dichotomy (critical collapse) states that \(8\pi \) is the critical mass for the simplest two-dimensional parabolic-elliptic KS system in \(\varOmega ={\mathbb {R}}^2\), namely (1) with \(\epsilon =k=0\),

The well-known KS dichotomy, also known as critical collapse, states that in the simplest 2D parabolic-elliptic KS system in the domain \(\varOmega ={\mathbb {R}}^2\), the critical mass is \(8\pi \). This system, denoted by (1) with \(\epsilon =k=0\), i.e.,

$$\begin{aligned} \rho _{t}&= \nabla \cdot ( \nabla \rho - \rho \, \nabla c ), \nonumber \\ \varDelta \, c&=- \rho , \end{aligned}$$
(5)

exhibits a critical behavior where the solution either remains global and bounded or blows up in finite time depending on the initial mass of the density \(\rho \). Specifically, we have

  1. 1.

    If \(M_0<8\pi \), the system has a global smooth solution.

  2. 2.

    If \(M_0>8\pi \), the system blows up in finite time in the sense of \(|\cdot |_\infty \) norm.

It can be seen from the classical variance identity for system (5), [27], that,

$$\begin{aligned} \frac{d}{dt}\int _{x\in {\mathbb {R}}^2} |x|^2 \, \rho (x)\, dx=\frac{M}{2\pi }(8\pi -M). \end{aligned}$$
(6)

Then, the solution of (5) exhibits a quantized concentration of mass at the origin, which is a type of blow-up known as a \(\delta \)-blow up.

In the case of the system (5) on \({\mathbb {R}}^d\) with \(d\ge 3\), the identity (6) is no longer applicable, and the evolution of the KS system is not as straightforward. However, the coexistence of blow-up and global smooth solutions still depends on the size of the initial data. For fully parabolic KS, [36, 37] show that the system may blow up in finite time over a large set of radial initial data. On the other side of the dichotomy, it has been shown [33] that global strong solutions exist for small initial data in the fully parabolic system (1).

In addition, the blow-up profile can differ from the \(\delta \)-type blow-up observed in 2D cases. For instance, it has been shown in [12] that in 3D parabolic elliptic systems, there exist radial, positive, backward self-similar solutions of the form,

$$\begin{aligned} \rho (x,t)=\frac{V(x/\sqrt{T-t})}{T-t}, \quad \quad 0<t<T, \end{aligned}$$
(7)

where the radially decreasing profile function V satisfies \(\lim _{y\rightarrow \infty }y^2V(y)=L\in {\mathbb {R}}^{+}\). Later in a more refined result by [29], the blowup is said to be type I if

$$\begin{aligned} 0<\limsup _{t\rightarrow T}\, (T-t)\, \Vert \rho \Vert _{\infty }\, <\infty . \end{aligned}$$
(8)

Then for radial initial data in \(L^1({\mathbb {R}}^3)\), if a blowup is type I, \(\exists \, C>0\) such that

$$\begin{aligned} \rho (x,t)\le C(T-t+|x|^2)^{-1}, \quad 0<|x|\le R, \quad 0<t<T. \end{aligned}$$
(9)

To the best of our knowledge, the complete characterization of blow-up criteria with non-radial initial data and the profile of the fully parabolic KS system in \({\mathbb {R}}^3\) has not been analyzed, also as discussed in [37]. Therefore, numerical computations are necessary to investigate the potential singular behavior, which we will discuss in Sect. 4.3.

2.2 Lagrangian Formulations

As a fundamental step of deriving the algorithms, we begin by introducing the Lagrangian formulation of the active particle density \(\rho \) in the KS system (1). We focus on the elliptic system with \(\epsilon =k=0\), specifically (5), which can be generalized to any dimension d. By considering the equation \(\varDelta c = -\rho \) and utilizing the Green’s function of the Laplacian operator in \({\mathbb {R}}^d\), we can deduce the following:

$$\begin{aligned} c(x,t)=\left\{ \begin{array}{l} -\frac{1}{2\pi }\int \ln |x-y| \, \rho (y,t), \quad d=2\\ C_d\int \frac{1}{|x-y|^{d-2}}\, \rho (y,t)\, dy, \quad d\ge 3 \end{array}\right. , \end{aligned}$$
(10)

where \(C_d=\frac{\varGamma (d/2+1)}{d(d-2)\pi ^{d/2}}\). So the convection term in (2) can be expressed as follows:

$$\begin{aligned} \nabla c(x) = -\frac{\varGamma (d/2)}{2\pi ^{d/2}}\int \frac{x-y}{|x-y|^d}\, \rho (y,t)\, dy. \end{aligned}$$
(11)

Now we arrive at the interactive stochastic differential equation system of P particles, \(\{X^p_t\}_{p=1:P}\),

$$\begin{aligned} dX^p_t=-\chi \frac{M}{P}\sum _{q\not =p} \frac{\varGamma (d/2)}{2\pi ^{d/2}} \frac{X^p_t-X^q_t}{|X^p_t-X^q_t|^d}+\sqrt{2\mu } \, dW^p_t, \quad p=1,\cdots ,P, \end{aligned}$$
(12)

where \(W^p_t\) denotes independent identically distributed standard Brownian motions. In [24], it has been demonstrated, under mild regularity condition, that as \(P\rightarrow \infty \), the distribution of empirical particles \(\{X^p_t\}_{p=1:P}\) converges to \(\rho \) in the continuous PDE system (2). To study the singularity behavior in the parabolic elliptic KS systems, several novel numerical methods have been developed and implemented, as discussed in [10, 21, 35].

In the fully parabolic case (\(\epsilon \not =0\)), the solution of chemical concentration c is obtained by solving a parabolic equation, which is no longer Markovian as in (12). In this case, at time \(t>0\), the solution of \(\rho \) in the interval [0, t] needs to be incorporated in the representation of c. Specifically, we have to consider the following:

$$\begin{aligned} c(\cdot ,t)=e^{-k^2t}e^{t\varDelta } c(\cdot ,0) + \int _0^t e^{k^2(s-t)}e^{(t-s)\varDelta }\, \rho (\cdot ,s)\, ds, \end{aligned}$$
(13)

where the heat semigroup operator \(e^{t\varDelta }\) is defined by

$$\begin{aligned} (e^{t\varDelta }f)(x,t):=\int \frac{e^{-|x-y|^2/(4t)}}{(4\pi t)^{d/2}}\, f(y)\, dy. \end{aligned}$$
(14)

Similar to (12), the empirical particle system converging to density \(\rho \) reads:

$$\begin{aligned} d X^{p}_t = \, {\chi } \nabla _{X}\, c(X^{p}_{t},t) \, dt + \sqrt{2\,\mu } \, d W^p, \;\; p=1,\cdots , P, \end{aligned}$$
(15)

and \(W^p\)’s are independent Brownian motions in \({\mathbb {R}}^d\). However, due to the historic path dependence in the solution of c in (13), direct computation of the drift \(\nabla _{X}\, c(X^{p}_{t},t) \) in (15) as (12) would result in significant memory and computational cost in each step, which increases with computational time T. More precisely, by directly substituting (13) into (15), the empirical particle system is as follows,

$$\begin{aligned} d X^{p}_t =&\, {\chi } \,e^{-k^2t}e^{t\varDelta }\nabla _{X} c(X^{p}_t,0) dt + \sqrt{2\,\mu } \, d W^p, \nonumber \\&-\left( \chi \frac{M}{P}\int _0^t e^{k^2(s-t)}\sum _{q\not =p} \frac{1}{(4\pi (t-s))^{d/2}}\frac{X^p_t-X^q_s}{2(t-s)}\, ds \right) dt . \end{aligned}$$
(16)

In the discrete scheme of (16), at the n-th temporal step, one shall compute the interaction of \(X_{t_n}\) with other particle positions over time interval \([0,t_n]\). At the n-th step, the computational cost is \({\mathcal {O}}(n P^2)\) and memory cost is \({\mathcal {O}}(nP)\).

A purely probabilistic particle method based on (16) for 2D fully parabolic KS system is proposed in [34]. The method is efficient due to the interaction, as mentioned by [34]. In contrast, [8] proposed a memory-less approach for 2D systems similar to ours. However, the method in [8] computes c only on the spatial grid points, which leads to numerical inaccuracy in moving the particles. More details for computing 2D systems can be found in [34]. To the best of our knowledge, a memory-less algorithm (memory cost \({\mathcal {O}}(P)\) for saving particle positions) to compute fully parabolic KS systems in 3D, or higher dimensions, has not been developed yet. We will present such an algorithm in the following section.

3 SIPF Algorithms for Parabolic-Parabolic KS

In this section, we will present the SIPF algorithm for solving the fully parabolic KS models. Since we are interested in studying the spatially localized aggregation behavior as discussed in Sect. 2.1, it is reasonable to consider the system (2) and (3) in a large domain \(\varOmega =[-L/2,L/2]^d\) and assume Dirichlet boundary condition for the particle density \(\rho \) and Neumann boundary condition for the chemical concentration c.

As a discrete algorithm, we assume that the temporal domain [0, T] is partitioned into time steps by \(\{t_n\}_{n=0:n_T}\), where \(t_0=0\) and \(t_{n_T}=T\). We approximate the density \(\rho \) using particles, i.e.

$$\begin{aligned} \rho _t \approx {\frac{M_0}{P}}\, \sum _{j=1}^{P} \delta (x - X^{p}_{t}), \; \; P\gg 1, \end{aligned}$$
(17)

where \(M_0\) is the conserved total mass (integral of \(\rho \)). For the chemical concentration c, we approximate it using a set of Fourier basis. In the following derivation, we assume \(d=3\) for brevity, while the algorithms proposed work under any dimension. Specifically, \(c({{\textbf {x}}},t)\) can be represented as a series expansion:

$$\begin{aligned} \sum _{j,m,l \in {\mathcal {H}}}\, \alpha _{t;j,m,l}\, \exp (i 2\pi j \, x_1/L) \exp (i 2\pi m\, x_2/L)\exp (i2\pi l\, x_3/L), \end{aligned}$$
(18)

where \({\mathcal {H}}\) denotes index set

$$\begin{aligned} {\mathcal {H}}:= \left\{ (j,m,l)\in {\mathbb {N}}^3: |j|,|m|,|l|\le \frac{H}{2}\right\} , \end{aligned}$$
(19)

and \(i=\sqrt{-1}\).

Then at \(t_0=0\), we generate P empirical samples \(\{ X^{p}_{0}\}_{p=1:P}\) according to the initial condition of \(\rho _0\), and set up \(\alpha _{0;j,m,l}\) using the Fourier series representation of \(c_0\).

For ease of presenting our algorithm, we will use a slight abuse of notation. We will represent the density \(\rho \) at time \(t_n\) as \(\rho _n={\frac{M_0}{P}}\, \sum _{p=1}^{P} \delta (x - X^{p}_{n})\), and the chemical concentration c at time \(t_n\) as

$$\begin{aligned} c_n=\sum _{j,m,l \in {\mathcal {H}}}\, \alpha _{n;j,m,l}\, \exp (i 2\pi j \, x_1/L) \exp (i 2\pi m\, x_2/L)\exp (i2\pi l\, x_3/L).\end{aligned}$$

To discretize the time-stepping system (1) from \(t_n\) to \(t_{n+1}\), with \(\rho _n\) and \(c_{n-1}\) known, our algorithm, inspired by the operator splitting technique, consists of two sub-steps: updating the chemical concentration c and updating the organism density \(\rho \).

Updating chemical concentration c Let \(\delta t=t_{n+1}-t_n>0\) be the time step. We discretize the equation for c in (1) in time using an implicit Euler scheme:

$$\begin{aligned} \epsilon \, ( c_n - c_{n-1} )/\delta t = (\varDelta - k^2 )\, c_n + \rho _n. \end{aligned}$$
(20)

From (20), we obtain the explicit formula for \(c_n\) as:

$$\begin{aligned} (\varDelta - k^2 -\epsilon /\delta t )\, c_n = - \epsilon \, c_{n-1}/\delta t - \rho _n. \end{aligned}$$
(21)

It follows that:

$$\begin{aligned} c_n=c({{\textbf {x}}},t_n)&= - {\mathcal {K}}_{\epsilon ,\delta t} *(\epsilon \, c_{n-1}/\delta t + \rho _n) = - {\mathcal {K}}_{\epsilon ,\delta t} *(\epsilon \, c({{\textbf {x}}},t_{n-1})/\delta t +\rho (x,t_n)) \end{aligned}$$
(22)

where \(*\) is the spatial convolution operator, and \({\mathcal {K}}_{\epsilon ,\delta t}\) is the Green’s function of the operator \(\varDelta - k^2 -\epsilon /\delta t\). In the case of \({\mathbb {R}}^3\), the Green’s function \({\mathcal {K}}_{\epsilon ,\delta t}\) is given by:

$$\begin{aligned} {\mathcal {K}}_{\epsilon ,\delta t} = {\mathcal {K}}_{\epsilon ,\delta t}({{\textbf {x}}})= - \frac{\exp \{ - \beta |{{\textbf {x}}} | \}}{ 4 \pi |{{\textbf {x}}}|}. \quad \beta ^2 = k^2 + \epsilon /\delta t. \end{aligned}$$
(23)

The Green’s function admits a closed-form Fourier transform,

$$\begin{aligned} {\mathcal {F}}{\mathcal {K}}_{\epsilon ,\delta t} (\omega )=-\frac{1}{|\omega |^2+\beta ^2}. \end{aligned}$$
(24)

For the term \(-{\mathcal {K}}_{\epsilon ,\delta t} *c_{n-1}\) in (22), using Eq.(24), it is equivalent to modify the Fourier coefficients \(\alpha _{j,m,l}\) to \(\alpha _{j,m,l}/(4\pi ^2j^2/L^2+4\pi ^2m^2/L^2+4\pi ^2l^2/L^2+\beta ^2)\).

For the second term \({\mathcal {K}}_{\epsilon ,\delta t} *\rho \), we first approximate \({\mathcal {K}}_{\epsilon ,\delta t}\) with a \(\cos \) series expansion. Then, according to the particle representation of \(\rho \) in (17), we have

$$\begin{aligned} \begin{array}{ll} & ({\mathcal {K}}_{\epsilon ,\delta t} *\rho )_{j,m,l} \approx \\ & \frac{M_0}{P}\sum _{p=1}^P \frac{\exp (-{i}2\pi j X^p_{n,1}/L-{i}2\pi m X^p_{n,2}/L-{i}2\pi l X^p_{n,l}/L)(-1)^{j+m+l}}{4\pi ^2j^2/L^2+4\pi ^2m^2/L^2+4\pi ^2l^2/L^2+\beta ^2}. \end{array} \end{aligned}$$

Finally, we summarize the one-step update of the Fourier coefficients of the concentration field c in Algorithm 1.

figure a

Remark 1

Note that (20)–(22) are algebraic derivations. In the vanishing \(\epsilon \) regime, the related semi-discrete systems tend to an elliptic equation. In terms of numerical implementation, the updates of Fourier coefficients with Green’s kernel by (23) and (3) are also consistent. Hence our method is robust in the vanishing \(\epsilon \) regime.

Updating density of active particles \(\rho \) In the one-step update of the density \(\rho _n\) represented by particles \(\{X_n^p\}_{p=1:P}\), we apply the Euler-Maruyama scheme to solve the SDE (15):

$$\begin{aligned} X^{p}_{n+1} = X^{p}_{n} +\chi \nabla _{{\textbf {x}}} c(X^{p}_{n},t_n)\delta t + \sqrt{2\,\mu \, \delta t} \, N^p_n, \end{aligned}$$
(25)

where \(N^p_n\)’s are i.i.d. standard normal distributions with respect to the Brownian paths in the SDE formulation (15). For \(n>1\), substituting (22) in (25) gives:

$$\begin{aligned} X^{p}_{n+1} = X^{p}_{n} -\chi \nabla _{{\textbf {x}}} {\mathcal {K}}_{\epsilon ,\delta t} *(\epsilon \, c_{n-1}({{\textbf {x}}})/\delta t + \rho _n({{\textbf {x}}}))|_{{{\textbf {x}}}=X^{p}_{n}}\delta t + \sqrt{2\,\mu \, \delta t} \, N^p_n, \end{aligned}$$
(26)

from which \(\rho _{n+1}({{\textbf {x}}})\) is constructed via (17).

In such a particle formulation, the computation of spacial convolution is slightly different from the one in the update of c, namely (22).

For \(\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t} *\, c_{n-1}(X^p_n)\), to avoid the singular points of \(\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t}\), we evaluate the integral with quadrature points that are away from 0. Precisely, we denote the standard quadrature points in \(\varOmega \) as

$$\begin{aligned} x_{j,m,l}=( j\, L/H,m\, L/H, j\, L/H), \end{aligned}$$
(27)

where j, m, l are integers ranging from \(-H/2\) to \(H/2-1\). When computing the integral \(\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t} *\, c_{n-1}(X^p_n)\), we evaluate \(\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t}\) at \(\{X^p_n+{\bar{X}}^p_n-x_{j,m,l}\}_{j,m,l}\) where a small spatial shift is given by

$$\begin{aligned} {\bar{X}}^p_n=\frac{H}{2L}+\left\lfloor \frac{X^p_n}{H/L}\right\rfloor \frac{H}{L}-X^p, \end{aligned}$$
(28)

and c at \(\{x_{j,m,l}-{\bar{X}}^p_n\}_{j,m,l}\) correspondingly. The latter one is computed by the inverse Fourier transform of the shifted coefficients, with \(\alpha _{j,m,l}\) modified to \(\alpha _{j,m,l}\exp (-i2\pi j{\bar{X}}^p_{n;1}/L- i2\pi m{\bar{X}}^p_{n;2}/L-i2\pi l{\bar{X}}^p_{n;3}/L) \), where \(({\bar{X}}^p_{n;i})\) denotes the i-th component of \({\bar{X}}^p_{n}\).

The computation of the term \(\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t} *\, \rho (X^p_n,t_{n})\) is straightforward thanks to the particle representation of \(\rho (X^p_n,t_{n})\) in (17):

$$\begin{aligned} \nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon ,\delta t} *\, \rho _n(X^p_n)=\int \nabla _{{{\textbf {x}}}}\, {\mathcal {K}}_{\epsilon ,\delta t}(X^p_n-y)\rho (y)\approx \sum _{q=1,q\not =p}^P\frac{M}{P}\nabla _{{{\textbf {x}}}}\, {\mathcal {K}}_{\epsilon , \delta t}(X^p_n-X^q_n). \end{aligned}$$
(29)

We summarize the one-step update (for \(n>1\)) of density in SIPF as in Algorithm 2.

figure b

Combining (22) and (26), we conclude that the recursion from

\((\{X^{p}_{n}\}_{p=1:P}, \rho _n({{\textbf {x}}}),c_{n-1}({{\textbf {x}}}))\) to \((\{X^{p}_{n+1}\}_{p=1:P}, \rho _{n+1}({{\textbf {x}}}),c_n({{\textbf {x}}}))\) is complete. We summarize the SIPF method in the following Algorithm 3.

figure c

Numerical Stability Under discretization \({\mathcal {H}}\), the update in Alg.1 is equivalent to,

$$\begin{aligned} {\hat{c}}_n(\omega )=\frac{1}{1+(|\omega |^2+k^2)\cdot \frac{\delta t}{\epsilon }}{\hat{c}}_{n-1}(\omega )+\frac{1}{|\omega |^2+k^2+\frac{\epsilon }{\delta t}}{\hat{\rho }}_n(\omega ), \end{aligned}$$
(30)

where for brevity we used the notation, \(\omega =\frac{2\pi }{L}(n,m,l)\). Given \(|{{\hat{\rho }}}_n(\omega )|\) bounded by \(M_\rho \),

$$\begin{aligned} |{\hat{c}}_n(\omega )| \le M_c + A_c^n(|{\hat{c}}_0|-M_c), \end{aligned}$$
(31)

where \(A_c=\frac{1}{1+(|\omega |^2+k^2)\cdot \frac{\delta t}{\epsilon }}<1\) and \(M_c = \frac{M_\rho }{|\omega |^2+k^2}\). Furthermore, in the discrete setting,

$$\begin{aligned} |{{\hat{\rho }}}_n(\omega )| = |\frac{M}{P}\sum _{p=1}^P \exp (-i2\pi j X^p_{n;1}/L-i2\pi m X^p_{n;2}/L-i2\pi l X^p_{n;3}/L)| \le M. \end{aligned}$$
(32)

(32) implies that one can take \(M_\rho =M\) (total mass) and validate that the stability of the update of c in Alg.1 only requires \(k>0\), independent of \(\delta t\) and position of \(X_n\).

In the one-step update of \(\rho _n\) represented by \(X_n\) via Alg.2, the numerical stability relates to the increment from drift terms, namely \(\frac{\chi M \delta t}{P}\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon , \delta t}(X^p_n-X^q_n)\) and \(\epsilon \chi (F,{\check{G}}) \frac{L^3}{H^3}\). The upper boundedness of the latter is ensured by the spatial shifting (28). For the direct interaction term \(\frac{\chi M \delta t}{P}\nabla _{{{\textbf {x}}}}{\mathcal {K}}_{\epsilon , \delta t}(X^p_n-X^q_n)\), due to the random normal update step, namely \(X_{n+1}^p\leftarrow X_{n+1}^p+\sqrt{2\mu \delta t}N\), the probability that \(X^p_n\), \(X^q_n\) coincides is zero. Practically, we clip \(X_{n+1}\) to the range of domain \(\varOmega \).

Remark 2

The convergence of the distribution of many-particle solution to the continuous PDE (the mean-field limit) in similar systems can be found in [2, 30]. In the mean-field regime, our discretization for the c equation, namely (3), is implicit Euler in time and Fourier spectral discretization in space. We validate the first order convergence in time in Sect. 4.2. We will leave the full convergence analysis of the SIPF as a future work.

Remark 3

The memory usage throughout Alg.3 is \({\mathcal {O}}(P+\operatorname {Card} {\mathcal {H}})\), where \(\operatorname {Card} {\mathcal {H}}\) denotes the total number of Fourier terms in the discretization (19). The computational cost for each step of Alg.2 and Alg.1 are \({\mathcal {O}}(P\operatorname {Card} {\mathcal {H}})\) and \({\mathcal {O}}(P^2 + P\operatorname {Card} {\mathcal {H}} \log (\operatorname {Card}{\mathcal {H}}))\). So the one step update in Alg.3 costs \({\mathcal {O}}(P^2 + P\operatorname {Card} {\mathcal {H}} \log (\operatorname {Card}{\mathcal {H}}))\). To be noted, the complete algorithm Alg.3 admits direct parallelization over particles in each step.

4 Numerical Experiments

4.1 Aggregation Behaviors

To illustrate the performance of the algorithm, we start with two examples. In both cases, the initial distribution \(\rho _0\) is assumed to be a uniform distribution over a ball centered at \(\begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix}\) with radius 1, as shown in Fig. 1a. Also, in both cases, we assume the following model parameters,

$$\begin{aligned} \mu =\chi =1, \quad \epsilon =10^{-4} \text { and } k=10^{-1}. \end{aligned}$$
(33)

for the fully parabolic KS model (1). These parameter choices in (33) are made so that the model exhibits comparable behavior to the corresponding parabolic-elliptic KS system whose blow-up behavior is known. For the first example, the total mass is chosen to be \(M_0=20\), while for the second one, the mass is \(M_0=80\).

In the numerical computation of both examples, we use \(H=24\) Fourier basis in each spatial dimension to discretize the chemical concentration c, and use \(P=10000\) particles to represent the approximated distribution \(\rho \). The computational domain is \(\varOmega =[-L/2, L/2]^3\) with \(L=8\). We compute the evolution of c and \(\rho \) using Algorithm 3 with a time step of \(\delta t=10^{-4}\) up to \(T=0.1\). The total computation time is 284 seconds on a workstation with one RTX2070 Super GPU card.

Fig. 1
figure 1

Density \(\rho \) approximated by empirical distribution at \(T=0.1\): the mass effect on focusing

Full size image

In Fig. 1, we plot the distribution \(\rho \) by empirical samples, at the starting time \(T=0\) and final computation time \(T=0.1\). In Fig. 1b, we observe the diffusive behaviors compared to the initial distribution shown in Fig. 1a. While in Fig. 1c, we increase the total mass from \(M_0=20\) to \(M_0=50\), we can see particles become concentrated at the origin, which indicates the possible blow-up of the continuous system.

Fig. 2
figure 2

Chemical concentration c at final time \(T=0.1\), sliced at \(z=0\): the mass effect on focusing

Full size image

In Fig. 2, we present the chemical concentration c at the final time \(T=0.1\) and third component \(z=0\) for various initial total masses \(M_0\). By comparing the sub-figures, we can see that in the case of a large total mass, c exhibits a sharp profile at the origin. This behavior, together with the near singular behavior of \(\rho \), is shown in Fig. 1b, c, indicating a possible blow-up.

Fig. 3
figure 3

Maximum of chemical concentration c vs computation time T with different total mass \(M_0\): identifying blow up by refining the discretization

Full size image

Furthermore, if we assume, there exists a self-similar profile of \(\rho \) at origin as discussed in [29] and Sect. 2.1, namely, \(\rho (x,t)\sim \frac{1}{|x|^2}\), the Fourier coefficients of the chemical concentration c have the following asymptotics according to Eq. (1):

$$\begin{aligned} {\mathcal {F}}c(\omega )\sim \frac{1}{|\omega |^2+k^2} {\hat{\rho }}\sim \frac{1}{(|\omega |^2+k^2)|\omega |}. \end{aligned}$$
(34)

Then, the maximum of c in the computation shall vary with the discretization parameter H. Specifically, we observe that at the origin,

$$\begin{aligned} c(0)\sim \int \frac{1}{(|\omega |^2+k^2)|\omega |}e^{i\omega x}d\omega |_{x=0} = \int \frac{1}{(|\omega |^2+k^2)|\omega |}d\omega . \end{aligned}$$
(35)

In practical discretization, the range of the integral (35) is determined by the maximum frequency, which can be expressed as \([-\frac{\pi }{L}(\frac{H}{2}-1),\frac{\pi }{L}\cdot \frac{H}{2}]^3\). Then, for the type of blow-up with a profile proportional to \(\frac{1}{|x|^2}\),

$$\begin{aligned} \Vert c\Vert _{\infty }={\mathcal {O}}(\ln (H)). \end{aligned}$$
(36)

A similar analysis shows that for the type of blow-up with a profile proportional to \(\delta (x)\), we have the following asymptotic relation:

$$\begin{aligned} \Vert c\Vert _{\infty }={\mathcal {O}}(H). \end{aligned}$$
(37)

In Fig. 3, we show the maximum value of c as a function of computational time T for different numbers of Fourier modes H and total mass \(M_0\). In the case of a possible blow-up (Fig. 3b), we can see that the maximum of c varies dramatically for different values of H. This variation can be used as an indicator of a possible blow-up, which will be further investigated in the following experiments.

Additionally, in the case where \(M_0=80\) and \(T=1\) are chosen to achieve numerical stability for \(\Vert c\Vert _\infty \), we conducted tests for H ranging from 8 to 24. In Fig. 4, we plot \(\Vert c\Vert _\infty \) against H and observe that the maximum value of c grows nearly linearly with respect to H. This further supports our previous observation that the maximum of c depends on the discretization parameter H.

Fig. 4
figure 4

Maximum of c scales linearly with the number of Fourier modes H (in each dimension) under total mass \(M_0=80\) (super critical): a \(\delta \) type blow up

Full size image

Remark 4

Similar ideas for detecting blow-ups by comparing maximum values computed under different discretizations can be found in the literature on the finite volume approach to 2D KS systems. For example, in [4], the \(\delta \)-type singularities in the 2D system are identified when \(\Vert \rho \Vert _\infty ={\mathcal {O}}(\frac{1}{\varDelta x \varDelta y})\).

4.2 Numerical Convergence

In this subsection, we validate the convergence of the algorithms numerically. Throughout this subsection, we consider the same initial condition (\(\rho \) and c at \(t=0\)) and physical parameters as in the first example. For the parameters in the discretization, we will take uniform time step \(\delta _{ref}=2^{-11}T\), the number of Fourier modes in each dimension as \(H=24\), the number of particles as \(P=10000\), and the computational domain as \(\varOmega =[-L/2,L/2]^3\) with \(L=8\). Additionally, we set \(M_0=80\) and \(T=0.01\) when the system has not formed any singularities (as shown in Fig.  3b).

For the investigation for convergence of \(\delta t\), we consider the time step \(\delta t\) in the range between \(2^{-8}T\) to \(2^{-4}T\), and we take \(\delta t_{ref}=2^{-11}T\) as the reference solution.

Another identity that can be used to validate the accuracy of the computation for the system (2)–(3) is the total concentration of the chemical attractant c at a given time t, which is given by

$$\begin{aligned} c_\textbf{0}(t):=\int _\varOmega c(x,t)dx. \end{aligned}$$
(38)

By integrating both sides of (3) over physical domain \(\varOmega \), we have,

$$\begin{aligned} \epsilon \frac{d}{dt}c_\textbf{0}(t)=-k^2c_\textbf{0}(t) + M_0, \end{aligned}$$
(39)

which implies,

$$\begin{aligned} c_\textbf{0}(t)=\left( c_\textbf{0}(0)-\frac{M_0}{k^2}\right) e^{-\frac{k^2}{\epsilon }t}+\frac{M_0}{k^2}. \end{aligned}$$
(40)

In Fig. 5, we present both types of error to verify the convergence of the algorithm. The first type of error is the \(L_2\) relative error of the chemical concentration c at the final time T compared to the reference solution. The second type of error is the mean squared relative error of the total concentration (38) over the interval [0, T] compared to the ground truth (40). To calculate the mean squared relative error of the total concentration \(c_\textbf{0}\), we use the following formula:

$$\begin{aligned} c_\textbf{0}\ \text {error} = \sqrt{ \frac{1}{T}\int _0^T \Big (\frac{|{\tilde{c}}_\textbf{0}(t)-c_\textbf{0}(t)|}{|c_\textbf{0}(t)|}\Big )^2 dt} \end{aligned}$$
(41)

where \({\tilde{c}}_\textbf{0}\) is the approximated value of \({c}_\textbf{0}\) obtained using the SIPF algorithm with a time step \(\delta t\). In addition, we fit the slope of error versus \(\delta t\) in the logarithmic scale and find that the \(L_2\) relative error \(e(\delta t)\) follows an approximate first-order convergence rate, with \(e(\delta t)={\mathcal {O}}(\delta t^{1.011})\). This indicates that the algorithm is approximately first order in time (Table 1).

Fig. 5
figure 5

Relative errors of c vs. \(\delta t\), compared with \(\delta = 2^{-11}\times 0.01\). Fitted rate: \(e(\delta t)={\mathcal {O}}(\delta t^{1.011})\)

Full size image
Table 1 Relative errors of c vs. \(\delta t\), from data of Fig. 5
Full size table

In addition, we investigate the convergence of the algorithm with respect to other parameters of discretization. To this end, we keep the reference setting as Fig. 5, and alter Fourier mode number H or particle number P to compare with the reference solution. In Fig. 6 and Table 2 we present the \(L_2\) relative error of \(c(\cdot ,T)\) with varying \(P=100,\ 200,\ 400,\ 800,\ 1600\) comparing with reference solution by \(P=10000\). The fitted convergence rate is \(e(P)={\mathcal {O}}(P^{-0.287})\), significantly smaller than standard Monte Carlo. We conjecture that this is due to particle interaction and leave it for a future study.

In Fig. 7 and Table 3 we present the \(L_2\) relative error of \(c(\cdot ,T)\) with varying \(H=6,\ 8,\ 10,\ 12,\) comparing with reference solution by \(H=24\). The fitted convergence rate is \(e(H)={\mathcal {O}}(H^{-0.659})\).

We report that the convergent behavior of \(c_\textbf{0}\) while varying H or P is similar.

Fig. 6
figure 6

Relative errors of c vs. P, compared with \(P=10000\)

Full size image
Table 2 Relative errors of c vs. P, from data of Fig. 6
Full size table
Fig. 7
figure 7

Relative errors of c vs. H, compared with \(H=24\). Fitted rate: \(e(H)={\mathcal {O}}(H^{-0.659})\)

Full size image
Table 3 Relative errors of c vs. H, from data of Fig. 7
Full size table

4.3 Blow up behaviors

As mentioned in Sect. 2.1, it is well-known that \(8\pi \) is the critical mass for the simplest two-dimensional parabolic-elliptic KS system (5). Specifically, the system exhibits the following dichotomy:

  1. 1.

    If \(M_0<8\pi \), the system has a global smooth solution.

  2. 2.

    If \(M_0>8\pi \), the system has no global smooth solutions and can exhibit blow-up behavior.

In the case of a fully parabolic KS system or the specific parabolic-elliptic KS system (5) with passive advection, there is no known variance identity similar to (6). Therefore, to investigate the possible blow-up behaviors, numerical computation becomes necessary. By utilizing the asymptotics described in (36) and (37), we can conduct tests for two scenarios: \(H=24\) and \(H=12\). In these examples, we will compare the \(\Vert c\Vert _\infty \) values to detect any potential blow-up occurrences.

Mass dependence We begin by investigating the critical mass \(M_0\), which plays a central role in the dichotomy of the simple 2D parabolic elliptic system (5). To this end, we initialize the algorithm with a uniform distribution over the unit ball centered at the origin and \(c(0,x)=0\). We then apply the algorithm with two different values of H to compute the density and chemical concentration until \(T=1\). To identify the possible blow-up, we compute the ratio of \(|c|_\infty \) between the two cases. In Fig. 8a, we present the ratio, \(\frac{|c|_{\infty , H=24}}{|c|_{\infty , H=12}}\), over time for various values of \(M_0\). We observe a sharp increase in the ratio when a potential blow-up forms for \(M_0\ge 47.6\). Figure 8b presents the ratio at the final time \(T=1\), indicating that the critical mass of the aforementioned initial condition should fall between 47.6 and 47.8.

In addition to the SIPF algorithm, we also present the numerical results obtained using the finite difference method (FDM). We note that the KS system (2)–(3) admits radial solutions when given constant scalar physical parameters (33) and a radially symmetric initial condition \((\rho _0,c_0)\). Therefore, we re-write the system in the radial coordinate,

$$\begin{aligned} \rho _t&=\mu \left( \rho _{rr} + \frac{2}{r}\rho _r\right) -\chi \left( \rho _r c_r +\rho (c_{rr}+\frac{2}{r}c_r\right) , \end{aligned}$$
(42)
$$\begin{aligned} \epsilon c_t&=c_{rr}+\frac{2}{r}c_r-k^2c+\rho , \end{aligned}$$
(43)
$$\begin{aligned} r&\in {\mathbb {R}}^+, \quad t \in [0,T]. \end{aligned}$$
(44)

To formulate a finite difference scheme, we consider the system (42)–(43) on the domain [0, 20] with the Neumann boundary condition. We use a uniform partition with \(N_r=2\times 10^5\) intervals for spatial discretization. For the temporal domain, we employ a backward Euler scheme with a time step of \(\delta t=10^{-5}\). This discretization method requires a comparable computational time (approximately 150 seconds) to the proposed SIPF. In Fig. 8c, we present the maximum value of c over time for various initial mass \(M_0\), denoted as

$$\begin{aligned} |c|_{\infty , FDM} = \sup _{t,r}|c(r,t)|. \end{aligned}$$
(45)

We found that the FDM exhibits numerical instabilities for initial masses between 47.6 and 47.8, which matches the prediction made by the SIPF algorithm in Fig. 8b. This example further validates the accuracy of the proposed SIPF algorithm. It is worth noting that our SIPF algorithm applies readily to more general (non-radial) KS systems, whereas FDM in 3D with a fine uniform mesh will be computationally much more expensive.

Fig. 8
figure 8

Ratio of \(|c|_\infty \)’s from 2 SIPF runs with \(H=24\) and \(H=12\) and \(|c|_\infty \) from FDM run: in radial case, the proposed SIPF can identify the same critical mass as FDM simulation

Full size image

Geometry dependence In contrast to the simplest parabolic-elliptic KS system (5), where the total mass is the sole determinant of the aggregation behavior, we have experimentally observed that the critical mass can vary for different initial distributions of \(\rho \). For instance, in an experiment aimed at identifying the critical mass (as shown in Fig. 8), we replaced the initial distribution with a uniform distribution on a ball centered at the origin with a radius of \({\textbf {0.8}}\). With a more concentrated initial distribution, we found that the critical mass of the system decreases. To be more precise, in Fig. 9a, we present the ratio of \(|c|_{\infty }\) for various total masses \(M_0\) as a function of computational time T. We can see a significant change in the ratio when the total mass \(M_0\) is large enough (\(M_0\ge 39\)), indicating the formation of potential singularities. Conversely, for relatively small values of \(M_0\) (\(M_0\le 38.8\)), the ratio remains stable around 1 throughout the computational time. In Fig. 9b, we present the ratio at the final time \(T=0.1\) as a function of the total mass \(M_0\), which indicates the critical mass for this particular initial condition lies between 38.8 and 39.

Fig. 9
figure 9

Ratio of \(|c|_\infty \)’s from 2 runs with \(H=24,12\); particles stay within initial radius 0.8: the more concentrated initial, the smaller critical mass

Full size image

Dependence on physical and biological parameters. Here we investigate the dependence of critical mass on other physical coefficients in the KS systems. Here we take the base configuration as specified in Sect. 4.1 with physical coefficients set (33). We then change one of the coefficients and apply the same methodology as the aforementioned example, i.e. calculating \(\frac{|c|_{\infty , H=24}}{|c|_{\infty , H=12}}\) for various values of the coefficient and find the first interval that \(\frac{|c|_{\infty , H=24}}{|c|_{\infty , H=12}}\) is significantly away from 1. In Table 4, we summarize the results by examples that change one of the coefficients in (33). From Table 4, we call tell that the following factors have a positive correlation with the critical mass (suppressing blow-up): mobility of the bacteria \(\mu \) and chemical decay constant k. In contrast, the following factors have a negative correlation (promoting blow-up): chemo-sensitivity of the bacteria \(\chi \) and time scale of chemotaxis \(\epsilon \).

Table 4 Dependence of critical mass on KS physical and biological parameters
Full size table

4.4 Aggregation Behaviors from Non-radial Initial Data

In this subsection, we investigate aggregation behaviors in more general distributions. To this end, we present two experiments: (1) the initial distribution is not centered at the origin or even not a standard Fourier collocation point; (2) a more practical scenario where the initial distribution \(\rho \) models several separated clusters of organisms.

Shifted initial distribution. In this example, we re-do the computation to find the critical mass in parameter setup (33) as the same approach as Sect. 4.3, while shifting the initial distribution concentrated at \([1/6,1/6,1/6]^T\). The concentration point is selected away from the standard Fourier collocation point in our maximal resolution, namely \([-4,4]^3\) domain with each direction discretized by \(H=24\) Fourier mode. Obviously, in this setup, the critical mass shall remain the same as the non-shifted one, namely (47.6, 47.8).

In Fig. 10a, we compute the ration \(\frac{|c|_{\infty , H=24}}{|c|_{\infty , H=12}}\) as in Fig. 8b. The projected critical mass falls between 47.6 and 47.8 identical as the experiment result in Sect. 4.3. In Fig. 10b we show the location of particles representing \(\rho \) at \(T=1\). More specifically, we compute a histogram of particles projected to x-y plane and zoom in to \([-1/3,1/3]^2\). Notice that \([-1/3,1/3]\) covers an interval of three Fourier collocation points under the resolution of computation.

Results in Fig. 10 confirm the capability of our method in predicting critical solutions that may focus on any location.

In essence, our algorithm approximates c with Fourier series (18) and update of c in Alg.1, always shifting the c series relative to the position of X during collocation, see (28). Hence our method avoids the potential numerical inaccuracy in [8] as pointed out in [34].

Fig. 10
figure 10

Results in the example of shifted initial distribution: SIPF predicts blow up with concentration not on the regular Fourier grids

Full size image

Multi-clustered initial distribution. This is a more practical scenario where the initial distribution \(\rho \) models several separated clusters of organisms. The mass in each cluster is below the critical mass, while the total mass is super-critical. To be more specific, we assume the initial distribution is a uniform distribution on four balls with a radius of 0.5. These balls are centered at four vertices of a regular tetrahedron, namely,

$$\begin{aligned} \begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}, \quad \begin{pmatrix} -\frac{1}{2}\\ \frac{\sqrt{3}}{2}\\ 0 \end{pmatrix}, \quad \begin{pmatrix} -\frac{1}{2}\\ frac{\sqrt{3}}{2}\\ 0 \end{pmatrix}, \quad \begin{pmatrix} 0\\ 0\\ \sqrt{2} \end{pmatrix}. \end{aligned}$$
(46)

See also Fig. 11a for the scatter plot of particles representing the initial distribution. In this case, we assume the total mass is \(M_0=80\) and each cluster has a mass of 20, which is below the critical mass for a ball with a radius of \(r=0.5\). Next, we apply the algorithm to compute the KS system up to \(T=0.5\), with two different spatial discretizations (\(H=24\) and \(H=12\)), while keeping the rest of the configurations. In Fig. 11b, we calculate the ratio between the maxima of c versus time for the two different spatial discretizations. We can see the singularities formed in the system at around \(T=0.3\).

Fig. 11
figure 11

Identifying the formation of a finite time singularity at \(t\approx 0.3\) in non-radial solutions

Full size image

In Fig. 12, we present the scatter plot of particles between the time \(T=0.1\) and \(T=0.4\). By comparing Fig. 11a with Fig. 12a, we can see diffusive behavior. This behavior is a result of the mass in each cluster being below the critical mass. The diffusive behavior persists until approximately \(T=0.2\), as depicted in Fig. 11b, where the active particles form a single larger cluster. The mass of this new cluster, centered at the origin, is \(M_0=80\). In Fig. 12c, we can observe the aggregation process starting to form a singularity. This can also be seen from the sharp increase in the ratio of the maximum of c in Fig. 11b. Finally, in Fig. 12d, we can directly identify the possible blow-up at the origin through the scatter plot.

Fig. 12
figure 12

Particle scatter plot at \(T=0.1:0.1:0.4\): three cluster merging and a singularity formation

Full size image

4.5 Critical Mass and Blowup in Parabolic-parabolic KS

As the last example, we examine the singular solutions in the fully parabolic KS systems. For the purpose of exposition, we set \(\epsilon =0.1\) in (1), while keeping the remaining physical parameters constant. The initial condition is assumed to be a uniform distribution on a ball with a radius of 0.8 and \(c(x,0)=0\).

From Fig. 9, we can determine that the critical mass is approximately \(M_0=39\). We apply the same computational configuration as in Fig. 9, except we enlarge the domain to \(L=12\) to accommodate possible diffusive behaviors. We test our algorithm for two cases, \(M_0=40\) and \(M_0=160\), respectively.

The behaviors of the system are reported in Fig. 13. In Fig. 13a and b, we present the scatter plot of the particles representing the density \(\rho \) with \(M_0=40\) and \(M_0=160\), respectively.

We find that despite the initial mass \(M_0=40\) being larger than the critical mass in the case of \(\epsilon =10^{-4}\), the system does not blow up. We report that the variance of the particles grows linearly in computational time T, with diffusion coefficients fitted to be 1.696. In the absence of the chemical attractant, namely \(\chi =0\), the diffusion coefficient is expected to be \(4\mu =4\). This implies that the parabolic-parabolic KS systems with mass below critical mass are effectively diffusive, with diffusion suppressed by chemical attraction.

However, for \(M_0=160\), the system exhibits a possible singularity at the origin. In Fig. 13c, we present the ratio of \(|c|_\infty \) under \(H=24\) and \(H=12\) for both initial masses. Similar to the observation in Fig. 13a and b, the blow-up behavior crucially depends on a critical level of the initial mass.

Fig. 13
figure 13

Effects of initial mass \(M_0\) on focusing behavior (finite time blowup)

Full size image

5 Concluding Remarks

In this paper, we developed a stochastic interacting particle and field algorithm, observed its convergence, and demonstrated its efficacy in computing blowup dynamics of fully parabolic KS systems in 3D from general non-radial initial data. The algorithm is recursive and does not have any history dependence, and the field variable is computed using Fourier series. Since the field variable (concentration) is smoother than the density, the series approach works well with only a few Fourier modes. The aggregation or focusing behavior in the density variable is resolved by using 10k particles. The algorithm successfully detects blowup through the field variable based on the critical amount of initial mass. The algorithm is self-adaptive and does not rely on any assumption about the blowup behavior, which is unknown except in the case of the parabolic-elliptic KS system. A potential weakness of the algorithm is the high cost of series expansion in 3D when a large number of Fourier modes is required for high-resolution computation near the blowup time. We will study this issue in future work.