1 Introduction↩︎

In the 1970s, Keller and Segel [1], [2] established a pioneering mathematical framework for chemotaxis. They formulated a system of nonlinear partial differential equations to represent the essential biological mechanism, in which cellular or organismal movement is directed by chemical cues that can be attractive or repulsive. Mathematically, the multi-species (\(d\)-species) Keller–Segel chemotaxis model is to find the cell (or organism) density functions \(\rho_i(\boldsymbol{x},t) ~ (i=1,\ldots,d)\) and the chemoattractant concentration function \(c(\boldsymbol{x},t)\) such that \[\label{model:dsKS} \left\{ \begin{align} &\partial_t \rho_i = \kappa_i\Delta \rho_i-\chi_i \nabla \cdot(\rho_i \nabla c), & \qquad \text{in}~~ \Omega \times(0, T], \\ &\partial_t c = \beta \Delta c- \alpha c+ \sum_{i=1}^d \gamma_i \rho_i, & \qquad \text{in}~~ \Omega \times(0, T]. \end{align} \right.\tag{1}\] Here \(\Omega \subset \mathrm{R}^2\) is assumed to be a two-dimensional convex, bounded and open domain. The parameters \(\kappa_i ~(i=1,\ldots,d)\) and \(\beta\) are positive diffusion coefficients, \(\chi_i>0~ (i=1,\ldots,d)\) is the chemoattractant sensitivity constant, \(\alpha \geq 0\) is the consumption rate of chemoattractant, and \(\gamma_i \geq 0 ~ (i=1,\ldots, d)\) represents the production rate of chemoattractant.

Without loss of generality, we only consider the two-species Keller–Segel chemotaxis model 1 , which involves identifying three real functions \(u = u(\boldsymbol{x}, t)\), \(v = v(\boldsymbol{x}, t)\) and \(c = c(\boldsymbol{x}, t)\) such that \[\label{model:KS} \left\{ \begin{align} & \partial_t u = \kappa_1\Delta u-\chi_1 \nabla \cdot(u \nabla c), & \qquad \text{in}~~ \Omega \times(0, T],\\ &\partial_t v = \kappa_2\Delta v-\chi_2 \nabla \cdot(v \nabla c), & \qquad \text{in}~~ \Omega \times(0, T],\\ & \partial_t c = \beta \Delta c- \alpha c+ \gamma_1 u +\gamma_2 v, & \qquad \text{in}~~ \Omega \times(0, T], \end{align} \right.\tag{2}\] subject to homogeneous Neumann boundary conditions \[\label{model:KS:bc} \begin{align} &\partial_{\boldsymbol{n}} u :=\nabla u \cdot \boldsymbol{n}=0,\quad\partial_{\boldsymbol{n}} v =0, \quad \partial_{\boldsymbol{n}} c=0, & & \text{ on } \partial\Omega \times(0, T], \end{align}\tag{3}\] and initial conditions \[\label{model:KS:ic} \begin{align} &u(\boldsymbol{x},0) =u^0(\boldsymbol{x}), \quad v(\boldsymbol{x},0) =v^0(\boldsymbol{x}), \quad c(\boldsymbol{x},0)=c^0(\boldsymbol{x}), & & \text{ in } \Omega, \end{align}\tag{4}\] where \(\boldsymbol{n}\) represents the unit outer normal vector onto the boundary.

Significantly, the Keller–Segel chemotaxis system 2 –4 obeys the mass conservation law, i.e., \[\label{model:KS:MC} \begin{align} M[u](t) &:= \int_\Omega u(\boldsymbol{x}, t) d\boldsymbol{x} = \int_\Omega u^0(\boldsymbol{x}) d\boldsymbol{x}=M[u](0),\\ M[v](t) &:= \int_\Omega v(\boldsymbol{x}, t) d\boldsymbol{x} = \int_\Omega v^0(\boldsymbol{x}) d\boldsymbol{x}=M[v](0). \end{align}\tag{5}\] Besides, for non-negative regular initial data 4 , i.e., \(u^0(\boldsymbol{x})\ge 0\), \(v^0(\boldsymbol{x})\ge 0\) and \(c^0(\boldsymbol{x})\ge 0\), the Keller–Segel chemotaxis system admits unique solutions with non-negative cell density and chemoattractant concentration, i.e., \[\label{model:KS:positivity} \begin{align} u(\boldsymbol{x},t) \geq 0, \quad v(\boldsymbol{x},t) \geq 0, \quad c(\boldsymbol{x},t) \geq 0, \quad \text{ in } \Omega. \end{align}\tag{6}\] Moreover, the Keller–Segel chemotaxis system can be viewed as a Wasserstein gradient flow driven by the total free energy \[\label{def:Energy} \begin{align} E[u,v,c] := \int_\Omega \left[\frac{\gamma_1 \kappa_1}{\chi_1}f(u) + \frac{\gamma_2 \kappa_2}{\chi_2}f(v) -\gamma_1 u c -\gamma_2 v c +\frac{\beta}{2} |\nabla c|^2 + \frac{\alpha}{2} c^2 \right] d\boldsymbol{x}, \end{align}\tag{7}\] where \(f(\rho) := \rho \log(\rho) - \rho\) with \(\rho \in (0,+\infty)\) for \(\rho=u, v\). It is easy to verify that the two-species Keller–Segel model is energy dissipative, i.e., \[\label{model:KS:Energy} \begin{align} \frac{dE[u,v,c]}{dt} = -\int_\Omega \left[\frac{\chi_1}{\gamma_1} u \Bigl(\nabla \frac{\delta E}{\delta u}\Bigr)^2 + \frac{\chi_2}{\gamma_2} v \Bigl(\nabla \frac{\delta E}{\delta v} \Bigr)^2 + \Bigl(\frac{\delta c}{\delta t}\Bigr)^2 \right] d\boldsymbol{x} \le 0, \end{align}\tag{8}\] where \(\frac{\delta E}{\delta u} = \frac{\gamma_1 \kappa_1}{\chi_1} f'(u) - \gamma_1 c\) and \(\frac{\delta E}{\delta v} = \frac{\gamma_2 \kappa_2}{\chi_2} f'(v) - \gamma_2 c\).

In recent years, considerable efforts have been devoted to the development of structure-preserving numerical methods for the Keller–Segel model that rigorously maintain the physical laws 5 , 6 , and 8 . For example, Chertock et al. [3] developed a fourth-order hybrid finite-volume–finite-difference scheme that can preserve positivity and mass conservation, while demonstrating high-order spatial accuracy and structure-preserving capabilities. Based on the Slotboom formulation [4]–[6], Liu et al. [7] designed second-order central difference schemes that can also preserve both positivity and mass conservation. Although optimal-order error estimates were not established, the authors provided a stability analysis for the proposed scheme. Moreover, inspired by the KKT-based positivity-preserving limiter algorithm, Cheng and Shen [8] developed a positivity/bound-preserving and mass-conservative Lagrange multiplier method for nonlinear parabolic systems. This method bypasses the need for complex nonlinear constrained optimization solvers, thereby substantially improving computational efficiency. Based on the scalar auxiliary variable (SAV) approach and function transformation, Huang and Shen [9] also constructed a positivity/bound-preserving, mass-conservative, and unconditionally modified energy-dissipative high-order time discretization scheme. Recently, Tong and Cai [10] constructed a Crank-Nicolson (CN) type finite difference scheme with second-order accuracy for the Poisson-Nernst-Planck equation, where a novel projection approach is adopted to ensure both positivity and mass conservation. Most importantly, they proved the optimal-order error estimates for the proposed scheme. This approach has also inspired us to develop structure-preserving numerical schemes for the two-species Keller–Segel chemotaxis model.

The block-centered finite difference (BCFD) method [11], also known as the cell-centered finite difference method, has been widely applied in recent years to the solution of various PDE models [11]–[17]. Notably, it can achieve second-order spatial accuracy on non-uniform spatial grids without sacrificing accuracy compared to standard finite difference schemes, and thus has great potential for simulating problems such as the Keller–Segel model with local blow-up solutions. As a positive first step, we recently developed a linearly implicit, fully decoupled CN-BCFD scheme on non-uniform spatial grids [18]. The proposed scheme not only guarantees mass conservation and second-order convergence, but also demonstrates a remarkable ability to capture blow-up phenomenon effectively and accurately. To further improve computational efficiency, it is better to adopt a time-staggered grid discretization that is able to facilitate linearization. In a recent paper [19], Zhang et al. proposed a uniformly time-staggered numerical scheme for the Keller–Segel–Navier–Stokes model, which achieves decoupling and linearization of the cell density and chemoattractant concentration. However, rigorous error estimates for the fully discrete scheme are lacking. Motivated by these observations, we propose a non-uniform time-staggered BCFD scheme to better capture blow-up phenomenon effectively and provide a rigorous theoretical analysis.

In this paper, we primarily focus on preserving positivity 6 and mass conservation 5 for the two-species Keller–Segel system 2 –4 . To this end, we employ the time-staggered BCFD discretization in which the cell densities are first computed via a prediction step, yielding solutions may not be positive; the predicted solutions are then projected onto a function space that enforces both positivity and mass conservation constraints via the standard \(L^2\) projection; and finally, the chemoattractant concentration is solved using the up-to-date cell densities. This approach ensures the desired physical properties, while requiring only the solutions of two simple nonlinear single-variable algebraic equations; see Eq. 21 of Remark 2 for reference. Moreover, the proposed time-staggered BCFD method with \(L^2\) projection is shown to be second-order accuracy in both time and space. In summary, this work presents a linearly implicit second-order time-staggered BCFD prediction-then-projection scheme with three main contributions:

• The proposed scheme is almost linear (expect for the efficient \(L^2\) projection step) and fully decoupled through a variable-step staggered-in-time discretization approach, which significantly enhances computational efficiency.

• The scheme unconditionally preserves positivity and mass conservation of the cell densities at the discrete level. In addition, the non-negativity of the chemoattractant concentration is ensured under a sufficient time-step condition.

• Optimal-order error analysis is rigorously established on non-uniform temporal grids, enabling the use of an efficient adaptive time-stepping strategy to accurately capture the blow-up phenomenon.

The rest of the paper is organized as follows. In Section 2, we propose a fully discrete time-staggered CN-BCFD scheme with an \(L^2\) projection strategy, and prove the positivity-preserving and mass-conservation properties at the discrete level for the two-species Keller–Segel chemotaxis model. Optimal-order error estimates together with the unique solvability of solutions to the proposed scheme are presented in Section 3. In Section 4, we present several numerical experiments to validate the accuracy, physical property-preserving properties, and reliability of the proposed scheme in simulating the blow-up phenomenon. Concluding remarks are given in Section 5. Throughout this paper, we denote by \(K\) with or without subscripts a generic positive constant that is independent of the grid parameters, but may have different values in different occurrences.

2 A non-uniform time-staggered BCFD prediction-then-projection scheme↩︎

This section is devoted to the construction of a time-staggered BCFD prediction-then-projection method for the two-species Keller–Segel chemotaxis system 2 –4 . For simplicity, below we assume \(\Omega :=(a^x, b^x)\times (a^y, b^y)\) and take the physical parameters \(\gamma_i = \kappa_i = \chi_i = \alpha = \beta \equiv 1\).

2.1 Notations and preliminaries↩︎

To fully decouple the concentration variable from the two-species density variables in 2 while enabling linearization, we employ a staggered-in-time variable-step Crank–Nicolson time-marching algorithm coupled with a non-uniform grid BCFD spatial discretization.

First, we introduce two distinct families of non-uniform partitions of \([0,T]\). The primal time levels are defined as \(0 = t_0<t_1<\cdots<t_N = T\), with stepsizes \(\tau_{n-1/2}=t_n-t_{n-1}\) for \(1\leq n \leq N\), and the cell density unknowns are evaluated at these mesh nodes. The staggered time levels are chosen as \(t_{n-1/2} = {(t_{n-1}+t_n)}/{2}\), the midpoint of each subinterval \([t_{n-1}, t_n]\) for \(1\leq n \leq N\), and the the chemoattractant concentration variable is approximated at these staggered nodes. Additionally, we set \(\tau_0 = t_{1/2}-t_0={\tau_{1/2}}/{2}\) and \(\tau_{n} =t_{n+1/2}-t_{n-1/2}={(\tau_{n+1/2} + \tau_{n-1/2})}/{2}\) for \(1\leq n \leq N-1\). Let the maximum stepsize be \(\tau:=\max_{1\leq n\leq N}\tau_{n-1/2}\). Furthermore, we assume that there exist two positive constants \(\sigma_* \le \sigma^*\) such that, for all \(1\leq n \leq N-1\), \[\sigma_* \tau_{n-1/2} \le \tau_{n+1/2} \le \sigma^* \tau_{n-1/2}.\] Given staggered-in-time grid functions \(\{g^n\}_{n\ge 0}\) and \(\{g^{n+1/2}\}_{n\ge 0}\), we define \[\begin{align} & d_\tau g^{n+1/2}= \frac{g^{n+1}-g^{n}}{\tau_{n+1/2}}, ~~ D_\tau g^{n}= \frac{g^{n+1/2}-g^{n-1/2}}{\tau_{n}}. \end{align}\]

Next, let \(N_x\) and \(N_y\) be the numbers of spatial grids along the \(x\)- and \(y\)-coordinates, respectively. Similar to those used in [20], [21], non-uniform staggered spatial grids are introduced. The primal grid points are denoted by \[\begin{align} \Pi_x: ~ a^x=x_{1 / 2}<x_{3 / 2}< \ldots < x_{i-1 / 2}<x_{i+1 / 2} < \ldots < x_{N_x+1 / 2}=b^x,\\ \Pi_y: ~ a^y=y_{1 / 2}<y_{3 / 2}< \ldots < y_{j-1 / 2}<y_{j+1 / 2} < \ldots < y_{N_y+1 / 2}=b^y, \end{align}\] with grid sizes \(h_i=x_{i+1/2}-x_{i-1/2}\) for \(i=1, \ldots, N_x\) and \(k_j=y_{j+1/2}-y_{j-1/2}\) for \(j=1, \ldots, N_y\). Let \(h:=\max \{h_i, k_j\}\). The middle grid points are denoted by \[\begin{align} \Pi_x^*: ~ x_i= (x_{i-1/2}+x_{i+1/2} ) / 2, ~~ i=1, \ldots, N_x,~~ \Pi_y^*: ~ y_j= (y_{j-1/2}+ y_{j+1/2} ) / 2, ~~ j=1, \ldots, N_y, \end{align}\] with grid sizes \(h_{i+1/2}=x_{i+1}-x_i= (h_{i+1} +h_i)/2\) for \(i=1,\ldots, N_x-1\) and \(k_{j+1/2}=y_{j+1}-y_j= (k_{j+1} +k_j)/2\) for \(j=1,\ldots, N_y-1\). Given spatial grid functions \(g=\{g_{i, j}\}\), \(\hat{g}=\{g_{i+1/2, j}\}\) and \(\check{g}=\{g_{i, j+1/2}\}\) defined on \(\Pi_x^*\times \Pi_y^*\), \(\Pi_x\times \Pi_y^*\) and \(\Pi_x^*\times \Pi_y\), respectively, we define \[\begin{align} & [d_x g]_{i+1/2, j}=d_x g_{i+1/2, j}=\frac{g_{i+1, j}-g_{i, j}}{h_{i+1/2}}, ~~ && [d_y g]_{i, j+1/2}=d_y g_{i, j+1/2}=\frac{g_{i, j+1}-g_{i, j}}{k_{j+1/2}}, \\ & \bigl[D_x \hat{g}]_{i, j}=D_x \hat{g}_{i, j}=\frac{\hat{g}_{i+1/2, j}-\hat{g}_{i-1/2, j}}{h_i}, ~~ && \bigl[D_y \check{g}]_{i, j}=D_y \check{g}_{i, j}=\frac{\check{g}_{i, j+1/2}-\check{g}_{i, j-1/2}}{k_j}. \end{align}\] Besides, we introduce the discrete inner products and norms on \(\Pi_x^*\times \Pi_y^*\), \(\Pi_x\times \Pi_y^*\) and \(\Pi_x^*\times \Pi_y\), respectively, as follows: \[\begin{align} & (f, g)_{\rm M}=\sum_{i=1}^{N_x} \sum_{j=1}^{N_y} h_i k_j f_{i, j} g_{i, j}, &&\|f\|_{\rm M}^2 = (f, f)_{\rm M}, \\ & (f, g)_x=\sum_{i=1}^{N_x-1} \sum_{j=1}^{N_y} h_{i+1/2} k_j f_{i+1/2, j} g_{i+1/2, j}, &&\|f\|_x^2 = (f, f)_x, \\ & (f, g)_y=\sum_{i=1}^{N_x} \sum_{j=1}^{N_y-1} h_i k_{j+1/2} f_{i, j+1/2} g_{i, j+1/2}, &&\|f\|_y^2 = (f, f)_y,\\ & (\boldsymbol{f}, \boldsymbol{g})_{\rm TM} = (f^x, g^x )_x+ (f^y, g^y)_y, &&\|\boldsymbol{f}\|_{\rm TM}^2 = (\boldsymbol{f}, \boldsymbol{f})_{\rm TM},\\ & \|\boldsymbol{d}g\|_\infty = \max _{i, j} |[d_x g]_{i+1/2, j}|+ \max _{i, j} |[d_y g]_{i, j+1/2}|. && \end{align}\]

The following two lemmas shall be used in the subsequent analysis.

Finally, given values \(\{p_{i,j} = p(x_i,y_j)\}\), for any points \((x, y) \in [x_i, x_{i+1}]\times [y_j, y_{j+1}]\), \(i=1,\ldots, N_x-1\), \(j=1,\ldots, N_y-1\), we introduce the piecewise bilinear interpolation function \(\ell_h p(x, y)\) by \[\label{Operator:bi} \begin{align} \ell_h p(x, y)& = \frac{(x_{i+1}-x)(y_{j+1}-y)}{(x_{i+1}-x_i)(y_{j+1}-y_j)}p_{i, j} + \frac{(x-x_i)(y_{j+1}-y)}{(x_{i+1}-x_i)(y_{j+1}-y_j)} p_{i+1, j} \\ & \quad + \frac{(x_{i+1}-x)(y-y_j)}{(x_{i+1}-x_i)(y_{j+1}-y_j)} p_{i, j+1} + \frac{(x-x_i)(y-y_j)}{(x_{i+1}-x_i)(y_{j+1}-y_j)} p_{i+1, j+1}. \end{align}\tag{9}\]

2.2 The PP-MC-PBCFD scheme↩︎

At each time step, the spatial discretization of the two-species Keller–Segel chemotaxis system 2 –4 is carried out using the non-uniform grid BCFD method. The primal scalar variables—namely, the cell densities and the chemoattractant concentration—are discretized on the staggered grid \(\Pi_x^*\times \Pi_y^*\), whereas the flux/gradient variables are approximated on \(\Pi_x\times \Pi_y^*\) (for \(x\)-direction) and \(\Pi_x^*\times \Pi_y\) (for \(y\)-direction). For temporal discretization, the cell densities are computed at the staggered time points \(t_{n+1/2}\) using a prediction-then-projection approach, and the chemoattractant concentration is subsequently solved at the primal time levels \(t_{n+1}\).

To be specific, let \(\mathbb{X}_h\) be the set of real-valued grid functions defined on \(\Pi_x^*\times \Pi_y^*\), i.e., \[\mathbb{X}_h:=\left\{v_h \mid v_h= \{v_{i,j}\}, ~~ i=1, \ldots, N_x, ~ j=1, \ldots, N_y\right\}.\] We denote the numerical solutions of the chemotaxis system 2 –4 by \(\{u_h^{n+1/2},v_h^{n+1/2},c_h^{n+1}\}\in \mathbb{X}_h \times \mathbb{X}_h \times \mathbb{X}_h\). Moreover, define a subspace of \(\mathbb{X}_h\) as \[\mathbb{X}_h[q]:=\left\{q_h \in \mathbb{X}_h \mid q_h \ge 0, ~ M_h[q_h^{n+1/2}] = M_h[q_h^0]\right\},\] where \(M_h[q]:= (q, 1 )_{\mathrm{ M }}\) represents the discrete version of mass. In addition, we denote the average density and concentration solutions at the temporal midpoints \((t_{n+1/2}+t_{n-1/2})/{2}\) and \(t_{n+1/2}\), respectively, as \[\overline{u}_h^{n} := \frac{1}{2}({u}_h^{n+1/2} + u_h^{n-1/2}), \quad \overline{v}_h^{n} := \frac{1}{2}({v}_h^{n+1/2} + v_h^{n-1/2}), \quad \overline{c}_h^{n+1/2} := \frac{1}{2}(c_h^{n+1}+c_h^n).\]

Now, let the initial approximations be given by \(\{u_h^{0},v_h^{0},c_h^{0}\}:=\{\ell_h u^{0}, \ell_h v^{0},\ell_h c^{0}\} \in \mathbb{X}_h[u] \times \mathbb{X}_h[v] \times \mathbb{X}_h\). For \(n \ge 0\), we propose the positivity-preserving, mass-conservative, projection-based BCFD scheme, termed the PP-MC-PBCFD scheme, as follows:

Step 1: For \(n=0\), solve the predicted density solutions \(\{\widetilde{u}_h^{1/2},\widetilde{v}_h^{1/2}\}\in \mathbb{X}_h \times \mathbb{X}_h\) at \(t=t_{1/2}\) via the semi-implicit Euler-BCFD method: \[\tag{10} \begin{align} \frac{\widetilde{u}_h^{1/2}-u_h^0}{\tau_0} &= D_x (d_x \widetilde{u}_h^{1/2}) + D_y(d_y \widetilde{u}_h^{1/2} ) - D_x([\ell_h \widetilde{u}_h^{1/2}][d_x c_h^0]) - D_y([\ell_h \widetilde{u}_h^{1/2}][d_y c_h^0]),\tag{11}\\ \frac{\widetilde{v}_h^{1/2}-v_h^0}{\tau_0} &= D_x(d_x \widetilde{v}_h^{1/2}) + D_y(d_y \widetilde{v}_h^{1/2}) - D_x([\ell_h \widetilde{v}_h^{1/2}][d_x c_h^0]) - D_y([\ell_h \widetilde{v}_h^{1/2}][d_y c_h^0]).\tag{12} \end{align}\] For \(n\ge 1\), given \(\{u_h^{n-1/2},v_h^{n-1/2}, c_h^{n-1}, c_h^{n}\}\in \mathbb{X}_h[u] \times \mathbb{X}_h[v] \times \mathbb{X}_h \times \mathbb{X}_h\), solve the predicted density solutions \(\{\widetilde{u}_h^{n+1/2},\widetilde{v}_h^{n+1/2}\}\in \mathbb{X}_h \times \mathbb{X}_h\) at the staggered time point \(t=t_{n+1/2}\) via the semi-implicit CN-BCFD method: \[\tag{13} \begin{align} \frac{\widetilde{u}_h^{n+1/2}-u_h^{n-1/2}}{\tau_{n}} & = D_x(d_x \overline{\widetilde{u}}_h^{n}) + D_y(d_y \overline{\widetilde{u}}_h^{n}) - D_x \bigl([\ell_h \overline{\widetilde{u}}_h^{n}][d_x c_h^{*,n}]\bigr) - D_y \bigl([\ell_h \overline{\widetilde{u}}_h^{n}][d_y c_h^{*,n}]\bigr),\tag{14}\\ \frac{\widetilde{v}_h^{n+1/2}-v_h^{n-1/2}}{\tau_{n}} &= D_x(d_x \overline{\widetilde{v}}_h^{n}) + D_y(d_y \overline{\widetilde{v}}_h^{n}) - D_x \bigl([\ell_h \overline{\widetilde{v}}_h^{n}][d_x c_h^{*,n}]\bigr) - D_y \bigl([\ell_h \overline{\widetilde{v}}_h^{n}][d_y c_h^{*,n}]\bigr),\tag{15} \end{align}\] where \(\overline{\widetilde{u}}_h^{n} := \frac{1}{2}(\widetilde{u}_h^{n+1/2} + u_h^{n-1/2})\) and \(\overline{\widetilde{v}}_h^{n} := \frac{1}{2}(\widetilde{v}_h^{n+1/2} + v_h^{n-1/2})\) represent the average values at the midpoint \(t_*:=\frac{t_{n+1/2}+t_{n-1/2}}{2}\), and the concentration variable in the nonlinear part is approximated by the linear extrapolation formula, i.e., \(c_h^{*,n} := \frac{t_n-t_*}{\tau_{n-1/2}}c_h^{n-1}+\frac{t_*-t_{n-1}}{\tau_{n-1/2}}c_h^n\) for \(n \ge 1\).

Step 2: Project the intermediate density solutions \(\{\widetilde{u}_h^{n+1/2},\widetilde{v}_h^{n+1/2}\}\) from \(\mathbb{X}_h \times \mathbb{X}_h\) to \(\mathbb{X}_h[u] \times \mathbb{X}_h[v]\) via the standard discrete \(L^2\) projection, and obtain the corrected density solutions \(\{u_h^{n+1/2},v_h^{n+1/2}\}\in \mathbb{X}_h[u] \times \mathbb{X}_h[v]\) such that \[\label{model:KS:shemeI2} \begin{align} & \min_{\{u_h^{n+1/2}, v_h^{n+1/2 }\}\in \mathbb{X}_h[u] \times \mathbb{X}_h[v]} \frac{1}{2} \left( \| u_h^{n+1/2} - \widetilde{u}_h^{n+1/2} \|^2 + \| v_h^{n+1/2} - \widetilde{v}_h^{n+1/2} \|^2 \right). \end{align}\tag{16}\]

Step 3: Solve the concentration solution \(c_h^{n+1}\in \mathbb{X}_h\) at the primal time level \(t=t_{n+1}\) via the CN-BCFD scheme: \[\label{model:KS:shemeI3} \begin{align} d_\tau c_h^{n+1/2}= D_x(d_x \overline{c}_h^{n+1/2}) + D_y(d_y \overline{c}_h^{n+1/2}) - \overline{c}_h^{n+1/2} + u_h^{n+1/2} + v_h^{n+1/2}. \end{align}\tag{17}\]

Below, we present a sufficient condition for the characterization of a non-singular M-matrix, which is the main tool for establishing positivity of the chemoattractant concentration variable.

We present the main structure-preserving conclusion for the PP-MC-PBCFD scheme as follows.

3 Unique solvability and error estimates↩︎

In this section, we show that the PP-MC-PBCFD scheme 10 –17 is second-order accurate in both time and space under uniform spatial grids and the following regularity assumptions: \[\label{eq:regular} c(\boldsymbol{x},t),u(\boldsymbol{x},t),v(\boldsymbol{x},t) \in W^{3,\infty}(0,T;W^{2,\infty}) \cap L^{ \infty}(0,T;W^{4,\infty}) \cap W^{2,\infty}(0,T;W^{3,\infty}).\tag{23}\] In particular, we assume that for some positive constant \(K_*\), there holds \[\label{eq:regular:1} \|\nabla c\|_{L^{\infty}(0,T;L^{\infty})} \leq K_*.\tag{24}\]

Let \(c^{n+1} = c(\cdot, t_{n+1})\), and \(u^{n+1/2} = u(\cdot, t_{n+1/2})\) and \(v^{n+1/2} = v(\cdot, t_{n+1/2})\) be the exact solutions of system 2 –4 at time \(t=t_{n+1}\) and \(t=t_{n+1/2}\), respectively. For the purpose of numerical analysis, we introduce two "biased" auxiliary solutions \[{\theta}_{\epsilon}^{n+1/2} := \bigl(1 + \epsilon_{\theta}^{n+1/2}\bigr) \ell_h\theta^{n+1/2},~~ \epsilon_{\theta}^{n+1/2} := \frac{M_h[\theta^0] - M_h[\theta^{n+1/2}] }{M_h[ \theta^{n+1/2}]}, \quad \text{for}~~ \theta=u, v.\] Note that under assumption 23 , for sufficiently small \(h\), we have \({u}_{\epsilon}^{n+1/2}, ~{v}_{\epsilon}^{n+1/2} \geq 0\) and \[\label{eq:epsilon} |\epsilon_{u}^{n+1/2}| + |\epsilon_{v}^{n+1/2}| \leq Kh^2,~~ |\epsilon_{u}^{n+1/2}-\epsilon_{u}^{n-1/2}|+|\epsilon_{v}^{n+1/2}-\epsilon_{v}^{n-1/2}| \leq K\tau_n h^2,\tag{25}\] for \(1\leq n \leq N-1\). In fact, the first part in 25 is a direct consequence of the second-order midpoint rule, triangle inequality and the mass conservation law 5 , i.e., \[\begin{align} |\epsilon_{u}^{n+1/2}| = \Big|\frac{M_h[u^0] - M_h[u^{n+1/2}] }{M_h[u^{n+1/2}]}\Big| &\leq K \big( | M_h[u^0]-M[u^0]| + | M[u^{n+1/2}] - M_h[u^{n+1/2}]| \big)\\ &\leq K \|u\|_{L^{\infty}(0,T;W^{2,1})} h^2. \end{align}\] For the second part in 25 , again by 5 and the standard midpoint rule error estimate, we get \[\begin{align} |\epsilon_{u}^{n+1/2}-\epsilon_{u}^{n-1/2}| &= \Big|M_h[u^0] \frac{M_h[u^{n-1/2}]-M_h[ u^{n+1/2}]}{M_h[u^{n-1/2}] M_h[u^{n+1/2}]}\Big| \leq K\big|M_h[u^{n-1/2}]-M_h[ u^{n+1/2}]\big|\\ &\leq K\big| \big(M[u^{n+1/2}] - M_h[u^{n+1/2}]\big)-\big(M[u^{n-1/2}] - M_h[u^{n-1/2}]\big)\big|\\ &\leq K \|u\|_{W^{1,\infty}(0,T;W^{2,1})} \tau_n h^2. \end{align}\]

Combining Taylor expansion, the Cauchy–Schwarz inequality, and Lemmas 2–3, we are led to estimates for the following local truncation errors: \[\begin{align} R_{u}^{1/2} &:= \frac{{u}_{\epsilon}^{1/2}-u^0}{\tau_0} - D_xd_x {u}_{\epsilon}^{1/2} - D_yd_y {u}_{\epsilon}^{1/2} - D_x \bigl([\ell_h {u}_{\epsilon}^{1/2}][d_xc^0]\bigr) - D_y \bigl([\ell_h {u}_{\epsilon}^{1/2}][d_y c^0]\bigr), \tag{26} \\ R_{v}^{1/2} &:= \frac{{v}_{\epsilon}^{1/2}-v^0}{\tau_0} - D_xd_x {v}_{\epsilon}^{1/2} - D_yd_y {v}_{\epsilon}^{1/2} - D_x \bigl([\ell_h {v}_{\epsilon}^{1/2}][d_x c^0]\bigr) - D_y \bigl([\ell_h {v}_{\epsilon}^{1/2}][d_y c^0]\bigr), \tag{27} \end{align}\] and for \(n \geq 1\), \[\begin{align} R_{u}^{n} &:= \frac{{u}_{\epsilon}^{n+1/2}-{u}_{\epsilon}^{n-1/2}}{\tau_{n}} - D_x d_x \overline{u}_{\epsilon}^{n} - D_y d_y \overline{u}_{\epsilon}^{n} - D_x \bigl([\ell_h \overline{u}_{\epsilon}^{n}][d_x c^{*,n}]\bigr) - D_y \bigl([\ell_h \overline{u}_{\epsilon}^{n}][d_y c^{*,n}]\bigr), \tag{28} \\ R_{v}^{n} &:= \frac{{v}_{\epsilon}^{n+1/2}-{v}_{\epsilon}^{n-1/2}}{\tau_{n}} - D_x d_x \overline{v}_{\epsilon}^{n} - D_y d_y \overline{v}_{\epsilon}^{n} - D_x \bigl([\ell_h \overline{v}_{\epsilon}^{n}][d_x c^{*,n}]\bigr) - D_y \bigl([\ell_h \overline{v}_{\epsilon}^{n}][d_y c^{*,n}]\bigr), \tag{29} \end{align}\] and for \(n \geq 0\), \[\begin{align} R_{c}^{n+1/2} := d_\tau c^{n+1/2}- D_x d_x \overline{c}^{n+1/2} - D_y d_y \overline{c}^{n+1/2} + \overline{c}^{n+1/2} - {u}_{\epsilon}^{n+1/2} - {v}_{\epsilon}^{n+1/2}. \label{eq:Rc95n} \end{align}\tag{30}\]

Next, set \({e}_{w}^{n+1/2} := {w}_{\epsilon}^{n+1/2} - w_h^{n+1/2}\), \(\widetilde{e}_{w}^{n+1/2}:={w}_{\epsilon}^{n+1/2} - \widetilde{w}_h^{n+1/2}\) for \(w=u, v\), and \({e}_{c}^{n+1} := \ell_hc^{n+1} - c_h^{n+1}\). It is straightforward to verify that \[\label{err:trunc:mass} (e_{w}^{n+1/2},1)_{\mathrm{ M }} =(\widetilde{e}_{w}^{n+1/2},1)_{\mathrm{ M }} = 0~~ \text{for} ~w =u,v.\tag{32}\] In addition, we define \({e}_{w}^{0} := {w}^{0} - w_h^{0}\) for \(w=u, v\). Note that \(w^0\) and \(w_h^0\) are equal at the grid points \(\Pi_x^*\times \Pi_y^*\), i.e., \(e_{w,i,j}^0 \equiv 0\) for \(i=1, \ldots, N_x, ~ j=1, \ldots, N_y\).

Now, by combining the results of Lemmas 9 and 10 together, and using the mathematical induction method, we proceed to prove the main convergence result.

4 Numerical results↩︎

In this section, we conduct several numerical experiments using the proposed PP-MC-PBCFD scheme (Eqs. 10 –17 ) to demonstrate its accuracy, verify the preservation of three important physical laws: positivity, mass conservation, and energy dissipation. Meanwhile, the blow‑up dynamics of the multi‑species Keller–Segel chemotaxis system on staggered non‑uniform spatio-temporal grids are also simulated. In the following, we introduce the non-uniform temporal and spatial grid partitions as follows: \[\label{eq:adap} \begin{align} t_{n} &= t_{{\rm fix},n} + \mu \, \Delta t \, (-1 + 2\,rand),\quad n=1,\ldots,N-1,\\ x_{i+1/2} &= x_{{\rm fix},i+1/2} + \nu \, \Delta x\, (-1+2\, rand),\quad i=1,\ldots,N_x-1,\\ y_{j+1/2} &= y_{{\rm fix},j+1/2} + \omega \, \Delta y\, (-1+2\, rand),\quad j=1,\ldots,N_y-1, \end{align}\tag{78}\] where \(t_{{\rm fix},n}= n\Delta t \,(n=0,\ldots,N)\), \(x_{{\rm fix},i+1/2}=a^x+i\Delta x\,(i=0,\ldots,N_x)\), \(y_{{\rm fix},j+1/2}=a^y+j\Delta y\,(j=0,\ldots,N_y)\), with uniform grid sizes \(\Delta t = {T}/{N}\), \(\Delta x= (b^x-a^x)/N_x\), and \(\Delta y= (b^y-a^y)/N_y\). Here, \(\mu\), \(\nu\) and \(\omega\) are small mesh parameters that control the extent of the random mesh perturbations within a specific range. In particular, when \(\mu=\nu=\omega=0\), 78 defines a uniform spatio-temporal partition. The symbol \(rand\) denotes a uniformly distributed random number in \([0,1]\). In what follows, we assume \(N_x=N_y (=N_z) \equiv M\) and \(\mu=\nu=\omega\).

This example is mainly used to test the accuracy of the PP-MC-PBCFD scheme on staggered non‑uniform spatio-temporal grids. To numerically evaluate both the spatial and temporal accuracy, we set the grid sizes in 78 such that \(\Delta x=\Delta t=1/M\), and measure the discrete \(L^2\) errors for the approximations of \(u\), \(v\), \(c\) and \(\nabla c\). As shown in Table ¿tbl:tab:accuracy?, the scheme clearly exhibits second‑order convergence in both time and space on uniform grids (i.e., \(\mu=0\)), which is fully consistent with the theoretical conclusion in Theorem 5. Furthermore, for small mesh perturbations with \(\mu=0.1\) and \(0.3\), although a rigorous convergence proof for non‑uniform spatial partitions is not yet available, the numerical results in Table ¿tbl:tab:accuracy? nonetheless indicate that the scheme maintains second‑order accuracy. This observation is accordance with the theoretical proof established for a linearized non-uniform grid mass-conservative BCFD scheme developed in [18] for the Keller–Segel chemotaxis model.

We use this example to verify the preservation of positivity and mass conservation properties of the PP-MC-PBCFD scheme 10 –17 under both uniform and non‑uniform spatio-temporal grids. Meanwhile, we also test the original energy dissipativity law of the proposed scheme, using the discrete version of energy 7 defined as \[\label{def:disEnergy} \begin{align} &E_h[u_h^n,v_h^n,c_h^n]\\ &\quad :=\big(u_h^n\log u_h^n + v_h^n\log v_h^n - u_h^n - v_h^n - u_h^n c_h^n - v_h^n c_h^n + \frac{1}{2}(c_h^n)^2,\boldsymbol{1}\big)_{\rm M} + \frac{1}{2}\big( (\boldsymbol{d} c_h^n)^2,\boldsymbol{1}\big)_{\rm TM}. \end{align}\tag{79}\] In the following simulation, the computational domain is discretized using \(M=80\) grid points in both the \(x\)- and \(y\)-directions, and the time stepsize in 78 is set to \(\Delta t=2.0\times10^{-3}\). The simulation results on both uniform grids (i.e., \(\mu=0\)) and non-uniform grids (i.e., \(\mu=0.1\)) are summarized in Figs. 1–2, which show the time evolution of the maximum and minimum values of \(u\), \(v\) and \(c\), along with the total mass of \(u\) and \(v\), and the iteration numbers of semismooth Newton solver for computing the \(L^2\) projection 16 . We have the following observations: (i) the maximum values of both \(u\) and \(v\) initially decrease while the minimum values initially increase; thereafter, both quantities gradually approach their steady states; (ii) the densities \(u\) and \(v\) remain positivity, and their masses are always conserved to at least 10 significant figures; (iii) the positivity of the concentration is also preserved throughout the simulation; (iv) the semi-smooth Newton method for solving the mass-conservative multiplier \(\xi\) and \(\theta\) in 21 converges in just one iteration per time step, demonstrating the efficiency of the proposed \(L^2\) projection process; and (v) as shown in Fig. 3, the original discrete energy 79 of the PP-MC-PBCFD scheme exhibits favorable dissipative behavior on both uniform (i.e., \(\mu=0\)) and non-uniform (i.e., \(\mu=0.1\)) grids, although a rigorous proof of such dissipative property is not established herein. To the best of our knowledge, comparable second-order linear schemes for the multi-species Keller–Segel model, together with rigorous positivity-preserving and energy-dissipation analysis, remain scarce in the literature.

To better simulate the blow-up phenomenon, we adopt the following time-adaptive strategy: \[\label{tau95adap} \tau_n := \min\bigg\{\max \bigg\{\frac{\tau_{max}}{\sqrt{1+\zeta\max\{(\|D_{\tau}u_h^{n-1}\|_{\infty},\|D_{\tau}v_h^{n-1}\|_{\infty})\}}},\tau_{min}\bigg\},\varsigma\tau_{n-1}\bigg\},\tag{80}\] where \(\tau_{max}\) and \(\tau_{min}\) denote the maximum and minimum allowable time stepsizes, respectively, and \(\zeta\) and \(\varsigma\) are two positive tunable parameters used to adjust the stepsize of the next time level. In addition, for this simulation, we adopt the following specially designed non-uniform spatial grids, which concentrate significantly more grid nodes near the blow-up point \((0.5,0.5)\) [18]: \[\label{grid:mid} \begin{cases} x_{N_{x}/2+i+1/2}=\frac{1}{2} + \frac{i^2}{2\left(N_x / 2+1\right)^2}, & i=0,1, \ldots, N_x / 2+1, \\ x_{N_{x}/2-i+1/2}=\frac{1}{2} -\frac{i^2}{2\left(N_x / 2+1\right)^2}, & i=1,2, \ldots, N_x / 2+1, \end{cases}\tag{81}\] with the \(y_{j+1/2}~(j=0,\ldots, N_y)\) grid points defined in an analogous manner.

In this test, we set \(M=80\) grid points in each spatial direction and adopt the grid partitions 81 along both the \(x\)- and \(y\)-directions. The adaptive time-stepping strategy 80 is employed with parameters \(\tau_{max}=1.0\times10^{-3}\), \(\tau_{min}=1.0\times10^{-5}\), \(\zeta = 1.0\times10^{-5}\) and \(\varsigma=5\). The simulation results including the time evolution of the extrema of \(u\), \(v\) and \(c\), the total mass of \(u\) and \(v\), and the iteration number of the semismooth Newton solver used in the \(L^2\) projection 16 , as well as the time evolution of the discrete energy are displayed in Figs. 4 and 5 (c), respectively. These results lead to conclusions similar to those drawn in Example 2. In addition, the evolution of the adaptive time stepsize is also illustrated in Fig. 5 (b). As can be observed, the time stepsize initially increases and then gradually decreases until it reaches the minimum value. In fact, it is precisely at this moment that the blow-up phenomenon occurs (see Fig. 5 (a)), and the minimum time stepsize \(\tau_{min}=1.0\times10^{-5}\) is thus required to adequately capture such blow-up phenomenon. This behavior is fully consistent with the expected physical scenario. Furthermore, the time-adaptive strategy 80 simulates the blow-up phenomenon in just \(86\) time steps. In contrast, a uniformly small stepsize \(\tau_{min} = 1.0\times10^{-5}\) requires \(1,860\) time steps. As a result, the adaptive approach is approximately 22 times more efficient. Finally, we also present in Fig. 6 the time evolution of the cell densities \(u\), \(v\), and the chemoattractant concentration \(c\). The results clearly show that both \(u\) and \(v\) undergo blow-up at \(t = 1.86\times10^{-2}\), which is in full agreement with the underlying physical mechanism.

In this test, we set \(M=40\) grid points in each spatial direction and choose the time stepsize \(\Delta t = 2.0\times10^{-3}\) as specified in 78 . Fig. 7 presents the simulation results on uniform grids (i.e., \(\mu=0\)), including the time evolution of the extrema of \(u\), \(v\) and \(c\), the total mass of \(u\) and \(v\), and the iteration number of the semismooth Newton solver used in the \(L^2\) projection 16 . Fig. 8 displays the corresponding results on non-uniform grids with \(\mu=0.1\). In addition, the time evolution of the original discrete energy on both uniform (i.e., \(\mu=0\)) and non-uniform (i.e., \(\mu=0.1\)) grids is shown in Fig. 9, which clearly illustrates its strict dissipative behavior. Furthermore, we also present slices at \(x = 0.5\), \(y = 0.5\) and \(z = 0.5\) of the cell densities \(u,v\) and the chemoattractant concentration \(c\) at four different time instants \(t = 0,\, 8.2 \times 10^{-3},\, 4.0 \times 10^{-2},\, 0.2\) in Fig. 10. Overall, the conclusions drawn here are largely consistent with those from the 2D Example 2.

5 Conclusion↩︎

This paper has introduced a fully decoupled, linearly implicit, positivity-preserving, and time-staggered BCFD prediction-then-projection scheme for the multi-species Keller–Segel chemotaxis system. We demonstrate several key features of the proposed scheme:

• The proposed scheme ensures unconditional positivity and mass conservation for the cell densities at the discrete level. Moreover, under a suitable time-step condition, the non-negativity of the chemoattractant concentration is also rigorously established.

• The unique solvability and optimal-order error estimates of the PP-MC-PBCFD scheme are rigorously established using the mathematical induction method and the discrete energy analysis approach, where the estimate established in Lemma 8 for the discrete \(L^2\) projection 16 plays a crucial role in the error analysis. As proved in Theorem 5, the cell densities \(u\) and \(v\) achieve second-order accuracy in both time and space in the discrete \(L^2\) norm, while the chemoattractant concentration \(c\) attains the same order in the discrete \(H^1\) norm.

• The use of variable time stepsize and time-staggered discretization fully decouples the solutions of the multi-species cell density variables and the chemoattractant concentration variable, while also enabling linearization, thereby significantly enhancing computational efficiency.

• An adaptive time-stepping strategy 80 , driven by the numerical solution evolution behavior, together with the time-staggered BCFD method on non-uniform spatial grids, effectively and accurately captures the blow-up phenomenon.

Furthermore, extensive numerical experiments have validated the accuracy, positive-preserving and mass conservation properties of the proposed scheme for the multi-species Keller–Segel chemotaxis system, while also demonstrating its reliability in simulating the blow-up phenomenon. As mentioned earlier, the model inherently satisfies the physical energy dissipation law [6], [7], [9], [28]–[30], which has been confirmed by our numerical results. However, a theoretical guarantee of this energy dissipation property is unavailable and requires further investigation.

CRediT authorship contribution statement↩︎

Ao Zhang: Methodology, Formal analysis, Software, Writing-Original draft. Bingyin Zhang: Methodology, Formal analysis, Writing-Original draft. Hongfei Fu: Conceptualization, Supervision, Writing-Reviewing and Editing, Methodology, Funding acquisition.

Declaration of competing interest↩︎

The authors declare that they have no competing interests.

Data availability↩︎

Data are available upon reasonable request.

Acknowledgements↩︎

This work was supported in part by the Natural Science Foundation of Shandong Province (No. ZR2024MA023) and by the National Natural Science Foundation of China (No. 12131014).

References↩︎