Stabilization and pattern formation in chemotaxis models with acceleration and logistic source

1.   Introduction

This study considers the system

in a bounded domain $\Omega\subset \mathbb{R}^{n}$ ($n\geq1$) with smooth boundary, where the parameters $\chi > 0$, $\alpha > 1$, $\gamma\geq0$ and $\tau_0\in\left\{{0, 1}\right\}$. The functions $u(x, t)$ and $v(x, t)$ denote the organism density and the concentration of the signal emitted by the organism itself respectively, the vector function $\mathbf{w}(x, t) = ({w}_1, {w}_2, \cdots, {w}_n)(x, t)$ represents the velocity of $u$. The system (1.1) is a generalized version of the classical Keller-Segel model (cf. ^[1])

which was used to describe the aggregation movement of cells. For half a century, this model and its variants have attracted a lot of attention and yielded rich results. (cf. the surveys ^[2,3,4,5] for instance).

In the classical chemotaxis model or other models containing taxis mechanisms such as the prey-taxis model (cf. ^[6]), it is generally assumed that the advective velocity of the organism is proportional to the concentration gradient of the stimulus/resource. However, this kind of assumption would neglect some details of the movement and is sometimes inappropriate to describe some observations, such as schooling fish (cf. ^[7,8]), swarms of flying insects (cf. ^[9]) and flea-beetles (cf. ^[10]), where the advective velocity variation (i.e., acceleration) rather than the velocity itself of the organism is supposed to be proportional to the gradient of stimuli/resources. This acceleration assumption was used in mathematical modeling of some predator-prey models, such as ^[11,12], where the model admits Spatio-temporal patterns consistent with experimental observations, while classical predator-prey models can not (cf. ^[13,14]). By introducing the acceleration assumption into chemotaxis models, ^[15] considered (1.1) with $\tau_0 = \gamma = 0$ under the boundary condition

where $\bf n$ is the outward normal vector on $\partial\Omega$ and $\partial_{\bf n} \mathbf{w}: = (\partial_{\bf n}{w}_1, \partial_{\bf n}{w}_2, \cdots, \partial_{\bf n}{w}_n)$. It is proved in ^[15] that for reasonable initial data, the system (1.1) with the above boundary condition admits global bounded solutions for $n\leq3$. This result is significantly different from classical parabolic-parabolic/parabolic-elliptic chemotaxis models, which have the well-known critical mass phenomenon in the two-dimensional case (cf. ^{[16,17,18,19,20]}) and admits blow-up solutions in three-dimensional case (cf. ^[18,21]). For critical mass phenomena in generalized versions of chemotaxis models, we refer the reader to ^{[22,23,24,25,26,27,28,29,30,31,32,33]} and the references therein.

Our main purpose is to study whether the fully parabolic ($\tau_0 = 1$) system (1.1) can preclude blow-up solutions (and has no critical mass phenomena for $n\leq3$) under the boundary condition

Specifically, we shall consider two cases: with/without the logistic source term. In the case of $\gamma = 0$, the main task focuses on whether there exist global bounded solutions for all reasonable initial data in two/three dimensions. When $\gamma > 0$, motivated by the works of the classical chemotaxis growth system (see Remark 1.1), our goal is to find the condition on $\alpha$ such that the model precludes blow-up solutions in high dimensions ($n\geq4)$. Moreover, for global bounded solutions, we shall investigate the large time behavior.

We assume that initial data $(u_0, v_0, \mathbf{w}_0)(x): = (u, v, \mathbf{w})(x, 0)$ conforms to

The first result is about the global boundedness of solutions, as follows.

Theorem 1.1. Suppose that $\tau_0 = 1$ and (1.3) holds with

then the system (1.1)-(1.2) possesses a unique classical solution $(u, v, \mathbf{w})$ satisfying

and $u, v > 0$ in ${\Omega} \times(0, \infty)$. Moreover, the solution is uniformly bounded in time, i.e.,

where $\beta > 1$ and $C > 0$ are two constants independent of $\chi$ and $t$.

Remark 1.1. Recall the classical chemotaxis growth system

with the homogeneous Neumann boundary condition and reasonable initial data. The parameters $h\in \mathbb{R}$, $\chi, \alpha > 0$ and $\tau_0\in\left\{{0, 1}\right\}$, and $\Omega\subset \mathbb{R}^n$ is a smooth bounded domain. In the case of $\alpha = 2$, (1.5) admits global bounded solutions for any $\mu > 0$ if $n = 2$ and for $\mu > \mu_0(\Omega)$ with some positive number $\mu_0(\Omega)$ if $n\geq3$ (cf. ^{[34,35,36,37,38]}). In the case of $\tau_0 = 0$, it is proved in ^[39] that there exist initial data such that solutions blow up in finite time if

From this point of view, Theorem 1.1 shows that the introduction of the acceleration assumption into chemotaxis models can improve the condition for $\alpha$ in some sense, such as (1.1) admits global bounded solutions for $\alpha > 1$ (or $\gamma = 0$) in three dimensions. Indeed, the condition for $\alpha$ in Theorem 1.1 is more correlated with that in the indirect chemotaxis model, see Remark 1.2 for a detailed reason.

Remark 1.2. Assume $\tau_0 = 1$. Now we consider a special case: $\mathbf{w}$ is a conservative vector field with $\mathbf{w} = \chi\nabla\phi$ for some scalar potential function $\phi\in {C^{1}(\overline\Omega)}$. Then the system (1.1) is related to the chemotaxis model with indirect signal production

If $\gamma = 0$, ^[40] shows that $n = 4$ is the critical dimension for (1.6) with the homogeneous Neumann boundary condition or a mixed boundary condition: no blow-up occurs if $\int_\Omega u_0 < \frac{(8\pi^2)}\chi$ (the radially symmetric setting is needed if the boundary condition is of the Neumann type). When $n\geq2$, $\gamma > 0$ and $\alpha > \frac12+\frac n4$, ^[41] proved that (1.6) admits global bounded solutions which convergence to $(1, 1, 1)$ in $L^{\infty}(\Omega)$ under the homogeneous Neumann boundary condition. For the system (1.1)-(1.2), Theorem 1.1 has the same restriction for $\alpha$ when $n\geq4$, but when $n = 2, 3$, the restriction can be improved to $\alpha > 1$ (or $\gamma = 0$). When ignore the random diffusion of $v$ in (1.6) (i.e., the third equation turns into $v_t = -v+u$), assuming $n\geq2$ and $\alpha > \frac n2$, ^[42] obtained global bounded solutions to (1.6) under the boundary condition $\frac{\partial u}{\partial\bf n} = \frac{\partial \phi}{\partial\bf n} = 0$ on $\partial\Omega$.

The global boundedness of solutions is given in Theorem 1.1. Next, we are dedicated to investigating the large time behavior of the obtained solution. The system (1.1)-(1.2) has a non-trivial constant equilibrium $\left(m, m, \bf0\right)$ with

We have the following asymptotic dynamics of the system (1.1)-(1.2).

Theorem 1.2. Suppose that the conditions in Theorem 1.1 hold, then there exists a number $\tilde{\chi} > 0$ such that for all $\chi\in(0, \tilde{\chi})$, the solution $(u, v, \mathbf{w})$ of the system (1.1)-(1.2) obtained in Theorem 1.1 satisfies

where $C$, $\sigma$ are two positive constants independent of $t$. Moreover, if $\gamma > 0$ and $\alpha\geq2$, then we can take $\tilde{\chi} = \sqrt\gamma$.

The above theorem gives the global stability of the equilibrium $(m, m, 0)$ under the condition that $\chi$ is small. However, outside the stable parameter regime, the behavior of solutions is unknown. It is well-known that the one-dimensional classical chemotaxis growth model can generate rich patterns, such as stationary, periodic, and chaotic (cf. ^[43]). To investigate whether the system (1.1)-(1.2) can generate similar patterns, we shall conduct linear analysis, weakly nonlinear analysis, and perform numerical simulations to see possible patterns.

Remark 1.3. The weakly nonlinear analysis in Sect. 4.2 shows that the asymmetric system (1.1)-(1.2) admits pitchfork bifurcations (which occur generically in symmetric systems).

Notations. For brevity, we shall use $C$ or $C_{i}$ $(i = 1, 2, 3, \cdots)$ denote a generic positive constant which may vary from line to line in the context, and use $C_P$ to denote the well-known Poincaré constant

The rest of this paper is organized as follows. Sect. 2 is devoted to establishing global bounded solutions to (1.1) and the proof of Theorem 1.1. By constructing appropriate Lyapunov functionals, the global stability stated in Theorem 1.2 is proved in Sect. 3. Outside the stable parameters regime, in Sect. 4, we first conduct linear analysis to specify the condition on parameters that allow pattern formation, then use a standard perturbation expansion approach (cf. [44,Sect. 6]) to conduct weakly nonlinear analysis to show pitchfork bifurcations, and finally perform numerical simulations for fully nonlinear parameter regimes to show that stationary patterns, spatially inhomogeneous time-periodic patterns, merging and emerging chaotic patterns are all possible. Moreover, we discussed some open questions for further research.

2.   Global boundedness

This section contributes to the proof of Theorem 1.1. To begin with, we give the existence of local solutions.

Lemma 2.1. Suppose that the conditions in Theorem 1.1 hold. Then there exists $T_{max} \in(0, \infty]$ such that the system (1.1) has a unique classical solution $(u, v, \mathbf{w})$ satisfying

and $u, v > 0$ in $\Omega \times(0, T_{max})$. Moreover,

Proof. This proof is based on standard arguments involving the contraction mapping principle. We omit the detailed proof here and refer the reader to ^[15] (the main difference here is that $\mathbf{w}$ satisfies a parabolic equation, however, the homogeneous Dirichlet boundary for $\mathbf{w}$ will provide sufficient regularities in the proof, see Lemma 2.4 below).

The following facts are basic but crucial.

Lemma 2.2. Suppose that the conditions in Theorem 1.1 hold. If $\gamma = 0$,

If $\gamma > 0$,

and

where $s: = \min\left\{{1, \frac{T_{max}}2}\right\}$.

Proof. Integrating the first equation of (1.1) with respect to $x\in\Omega$, one has

If $\gamma = 0$, (2.1) follows immediately. If $\gamma > 0$, by Hölder's inequality $\|u\|_{L^{1}(\Omega)}\leq |\Omega|^{1-\frac1\alpha}\|u\|_{L^{\alpha}(\Omega)}$ thanks to $\alpha > 1$, we obtain

which proves (2.2). Integrating (2.4) in time and using (2.2), one has (2.3).

To derive a priori estimates for $v$, we shall need the following property (cf. [46,Lemma 4.1] and [49,Lemma 1]), which is derived based on smooth properties of the Neumann heat semigroup.

Lemma 2.3. Let $\Omega$ be a bounded domain in $\mathbb{R}^{n}$, $n\geq1$, with smooth boundary. Let $T \in(0, \infty]$ and suppose that $z \in C^{0}(\overline{\Omega} \times[0, T)) \cap C^{2, 1}(\overline{\Omega} \times(0, T))$ is a solution of

where $g \in L^{\infty}\left((0, T); L^{q}(\Omega)\right)$ with some $q\geq1$. Then there exists a constant $C > 0$ such that

In addition to the well-known smooth properties of the Neumann heat semigroup [63,Lemma 1.3], we will use the $L^{p}$-$L^{q}$-estimates of the Dirichlet heat semigroup (cf. [55,Lemma 2.4 (i)], where the proof is based on [58,Proposition 48.4*,48.5 and 48.7*] and similar arguments as in the proof of [63,Lemma 1.3]).

Lemma 2.4. ([55,Lemma 2.4 (i)]). Let $e^{t {\Delta}}$ be the Dirichlet heat semigroup in $\Omega\subset \mathbb{R}^n$ $(n\geq1)$, $\lambda_1 > 0$ denote the first nonzero eigenvalue of $-\Delta$ in $\Omega$ under the Dirichlet boundary condition. If $1 \leq q\leq p \leq \infty$, then for any $z \in L^{q}(\Omega)$, it holds that

and

We are in a position to derive the following boundedness criterion.

Lemma 2.5. Suppose that $n\geq1$, $\gamma\geq0$ and $\alpha > 1$. If there exists $p > n$ such that

for some positive constant $c_0$, then there exist two constants $\beta > 1$ and $C > 0$ independent of $\chi$ and $t$ such that

for all $t\in(0, T_{max})$.

Proof. This proof can follow similar arguments developed in [59,Lemma 4.2] (see also [2,Lemma 3.2], [47,Lemma 3.1] and [48,Lemma 3.2] for instance). However, in order to give a clear dependency between the boundedness of $\sup_{t \in(0, T_{max})}\|u(\cdot, t)\|_{L^{\infty}}$ and $\chi$, we shall sketch the proof here. Now for each $T \in\left(0, T_{\max }\right)$, define

which is finite thanks to the local existence results in Lemma 2.1. Since $p > n$, for any $q\in(n, p)$, it holds that

Then by Hölder's inequality, the interpolation inequality, Lemma 2.2, (2.5) and (2.7), we have

For $t_0: = (t-1)_+$ and each $t\in(0, T)$, the variation-of-constants formula for $u(\cdot, t)$ gives

Since $\alpha > 1$, we have $\gamma\left({u-u^\alpha}\right)\leq \gamma$, then using (2.8), the maximum principle, $L^p$-$L^q$ estimates of the Neumann heat semigroup (cf. [63,Lemma 1.3]) and $\frac12+\frac n{2q} < 1$, one has

where $(e^{t \Delta})_{t \geq 0}$ is the Neumann heat semigroup in $\Omega$ and $\lambda_2$ denotes the first nonzero eigenvalue of $-\Delta$ in $\Omega$ under the homogeneous Neumann boundary condition. In view of (2.9), we arrived at

which along with $\delta\in(0, 1)$ gives

Since the above $C_i$ ($i = 1, 2, 3$) is independent of $T$, we have

With (2.10), we can apply Lemma 2.3 to prove

which together with Lemma 2.4 shows that

The combination of (2.10)–(2.12) proves (2.6) by letting $\beta = \frac1{1-\delta}$.

The following conclusion follows from Lemma 2.2 and Lemma 2.5 immediately.

Corollary 2.6. If $\mathbf{(H1)}$ holds, then there exist two constants $\beta > 1$ and $C > 0$ independent of $\chi$ and $t$ such that (2.6) holds.

Proof. It holds that $\frac n{n-1}\geq\frac 32$ thanks to $n\leq3$. This along with Lemma 2.2 and Lemma 2.3 indicates that

The variation-of-constants formula for $\mathbf{w}$ yields

where $(e^{t \Delta})_{t \geq 0}$ is the Dirichlet heat semigroup in $\Omega$. Using Lemma 2.4, (2.13) and the fact that $\frac n2\left({\frac23-\frac14}\right) < 1$ for $n\leq3$, for all $t\in(0, T_{max})$, we have

In view of Lemma 2.5 and $n\leq3 < 4$, the proof is obtained immediately from (2.15).

Lemma 2.7. Suppose that $n\geq2$, $\gamma > 0$, $\alpha > \frac12+\frac n4$, then there exist two constants $\beta > 1$ and $C > 0$ independent of $\chi$ and $t$ such that (2.6) holds.

Proof. For $\varepsilon > 0$, multiplying (2.4) by $e^{ \varepsilon t}$ yields

For $t_0\in(0, T_{max})$, integrating the above equation in time and using (2.2), we have

which implies that

Based on (2.16) and the smoothing properties of the Neumann heat semigroup, we can use similar arguments as in [41,Lemma 3.4] to find some $\rho_1 > \frac n2$ with $\alpha > \frac{(n+2)\rho_1}{n+2\rho_1}$ such that

Letting $\rho_2: = \frac {n\rho_1}{(n-\rho_1)_+} > n$, then using (2.14), (2.17) and Lemma 2.4, we have

This alongside Lemma 2.5 completes the proof.

Proof of Theorem 1.1. If $\mathbf{(H1)}$ holds, then by Corollary 2.6 we arrive at (2.6). If $\mathbf{(H2)}$ holds, then an application of Lemma 2.7 also gives (2.6). In view of Lemma 2.1 and (2.6), we obtain $T_{max} = \infty$. Then (1.4) follows from (2.6) immediately, which completes the proof.

3.   Global stability

In this section, we shall construct appropriate Lyapunov functionals to prove the large time behavior of solutions stated in Theorem 1.2. In the following, we suppose that the conditions in Theorem 1.1 hold and $(u, v, \mathbf{w})$ is the solution of the system (1.1)-(1.2) obtained in Theorem 1.1. The proof of Theorem 1.2 is mainly divided into two cases: with/without the logistic source term. We first give the following higher order regularities of the solution.

Lemma 3.1. For any $\theta\in(0, 1)$, there exists a positive constant $C(\theta)$ depending on $\theta$ such that

Proof. The conclusion can be proved by standard bootstrap arguments based on (1.4), the $L^{p}$ estimate and the Schauder estimate (cf. ^[51]). We omit the details here for brevity and refer the reader to [60,Theorem 2.1] for instance. The following auxiliary lemma will be used to prove global stability.

Lemma 3.2. ([61,Lemma 1.1]) Let $a \geq 0$ and $b > 0$ be two constants, $F(t) \geq 0, \int_{a}^{\infty} H(t) \mathrm{d} t < \infty$. Assume that $E \in C^{1}([a, \infty)), E$ is bounded from below and satisfies

If either $F \in C^{1}([a, \infty))$ and $F^{\prime}(t) \leq k$ in $[a, \infty)$ for some $k > 0,$ or $F \in C^{\theta}([a, \infty))$ and $\|F\|_{C^{\theta}([a, \infty))} \leq k$ for some constants $0 < \theta < 1$ and $k > 0$, then

With the aid of the Poincaré inequality (1.9), we can construct the following Lyapunov functional in the case of $\gamma = 0$.

Lemma 3.3. Suppose that $\mathbf{(H1)}$ holds with $\gamma = 0$, then there exists a positive constant $\tilde\chi_1 > 0$ such that for all $\chi\in(0, \tilde\chi_1)$, the function

satisfies

where $\eta_1: = \min\left\{{\frac1{3C_P}, \frac12}\right\}$, $m$ and $C_P$ are given by (1.7) and (1.9) respectively.

Proof. Integrating the equations in (1.1) by parts, using (1.9) and Young's inequality, for all $t > 0$, we obtain

and

as well as

The combination of (3.2)–(3.4) implies that

According to (1.4), we know there exists $\tilde\chi_1 > 0$ such that

which indicates that

This proves (3.1).

In the presence of logistics sources, i.e., $\gamma > 0$, we can construct the following Lyapunov functional.

Lemma 3.4. Suppose that $\gamma > 0$ and either $\mathbf{(H1)}$ or $\mathbf{(H2)}$ holds, then there exists a positive constant $\tilde\chi_2 > 0$ such that for all $\chi\in(0, \tilde\chi_2)$, the function

satisfies

where

Proof. Multiplying the first equation in (1.1) by $1-\frac1u$, integrating the result by parts and using Young's inequality, for all $t > 0$, we obtain

where $\varphi(s, \alpha): = \frac{s^{\alpha-1}-1}{s-1}$ for $s\geq0$ and $\alpha > 1$ with $\varphi(1, \alpha): = \lim_{s\rightarrow 1}\varphi(s, \alpha) = \alpha-1$. Similarly, as in deriving (3.3), we have

For all $t > 0$, the combination of (3.4), (3.6) and (3.7) implies that

{Case1: $1 < \alpha < 2$}. It holds that $\varphi(u, \alpha) = (\alpha-1)s_0^{\alpha-2}$ with some $s_0$ between $u$ and $1$. By (1.4) we know that there exists a constant $C_1 > 1$ independent of $\chi\in(0, 1)$ such that $u(x, t) < C_1$ for all $(x, t)\in\Omega\times(0, \infty)$. Therefore, if we take $\tilde \chi_2 = \min\left\{{1, (\gamma(\alpha-1)C_1^{\alpha-2})^\frac12}\right\}$, then

which along with (3.8) proves (3.5).

{Case2: $\alpha\geq 2$}. We first note that $\varphi(0, \alpha) = 1$ and $\varphi(1, \alpha) = \alpha-1\geq1$. Since $\frac d{d\alpha}\varphi(s, \alpha) = \frac{s^{\alpha-1}\ln s}{s-1} > 0$ for $s\neq0, 1$, we have

Let $\tilde \chi_2 = \sqrt\gamma$, then

which alongside (3.8) also arrives at (3.5).

We are now in a position to prove the convergence rate of solutions stated in Theorem 1.2.

Proof of Theorem 1.2. If $\mathbf{(H1)}$ holds with $\gamma = 0$, then it follows from Lemma 3.3 that (3.1) holds for $0 < \chi < \tilde\chi: = \tilde\chi_1$, hence we have

We have from (3.11) that

which along with Lemma 3.1 and the Gagliardo-Nirenberg inequality

proves (1.8) with $\sigma: = \frac{2\eta_1}{n+4}$.

We next assume that $\gamma > 0$ and either $\mathbf{(H1)}$ or $\mathbf{(H2)}$ holds, then Lemma 3.4 indicates that (3.5) holds for $0 < \chi < \tilde\chi: = \tilde\chi_2$. Noting that $\mathcal E_2(t)\geq0$ for all $t > 0$ due to $s-1-\ln s\geq0$ for $s\geq0$. Using Lemma 3.2, Lemma 3.1 and Lemma 3.4, we obtain $\mathcal F(t)\rightarrow0$ as $t\rightarrow \infty$, which implies that

The L'H{$\rm\hat{o}$}pital's rule shows that $\lim\limits_{u \rightarrow 1} \frac{u -1- \ln u}{\left(u -1\right)^{2}} = \lim\limits_{u \rightarrow 1} \frac{1}{2 u } = \frac{1}{2 },$ hence we can find $t_1 > 1$ such that

This together with Lemma 3.4 shows that there exists a constant $\eta_2 > 0$ such that

We are now at the same position as in (3.11), which also proves (1.8) with $\sigma: = \frac{2\eta_2}{n+4}$. Moreover, the proof of Lemma (3.4) (see (3.9) and (3.10)) implies that if $\gamma > 0$ and $\alpha\geq2$, then $\tilde{\chi}$ can be taken as $\tilde{\chi} = \sqrt{\gamma}$.

4.   Linear stability analysis and spatio-temporal patterns

A variety of patterns (including stationary, periodic, and chaotic) are observed in the one-dimensional classical Keller-Segel model with the logistic term (cf. ^[43]). Theorem 1.2 shows that the equilibrium $(m, m, \textbf0)$ of (1.1) is globally asymptotic stable for appropriately small $\chi > 0$. This raises a natural question, whether the considered system (1.1) with the acceleration assumption can generate aggregated patterns for large $\chi$? To this end, we shall first conduct linear analysis outside the stable parameters regime to specify the condition on parameters such that the equilibrium $(m, m, \textbf0)$ is linearly unstable, and then perform numerical simulations to see possible patterns. Throughout this section, we assume that $n = 1$, $\Omega = (0, L)$ is an interval of length $L > 0$ and $\tau_0 = 1$ in (1.1).

4.1.   Linear stability analysis

It is well known that the Neumann and Dirichlet eigenvalue problems

have eigenvalues $0 = \rho^{(N)}_0 < \rho^{(N)}_1\leq \rho^{(N)}_2\leq \cdots\leq\rho^{(N)}_k\leq\cdots$ and $0 < \rho^{(D)}_1\leq \rho^{(D)}_2\leq \cdots\leq\rho^{(D)}_k\leq\cdots$ with corresponding eigenvectors $\left\{{\varphi_k}\right\}_{k = 0}^\infty$ and $\left\{{\psi_k}\right\}_{k = 0}^\infty$ respectively, where $\left\{{\varphi_k}\right\}_{k = 0}^\infty$ and $\left\{{\psi_k}\right\}_{k = 1}^\infty$ respectively form an orthonormal basis of ${L^{2}(\Omega)}$, and $\varphi_0 = \sqrt{\frac1L}$. Moreover,

where $\omega_k: = \frac{k\pi}L\text{ for all }k\in\mathbb N$. Then a solution of (1.1)-(1.2) can be written as

where $m$ is given by (1.7), and we have supplemented the definition

so that the forms of three component solutions are uniform (noting that $\sum_{k = 0}^\infty d_k(t)\psi_k(x) = \sum_{k = 1}^\infty d_k(t)\psi_k(x)$ for any function $d_0(t)$ of $t$).

For linear stability, we assume $|a_k|$, $|b_k|$, $|d_{k+1}|\ll1$ for each mode $k\in\mathbb N$ to represent $(u, v, w)$ has a small perturbation to the equilibrium $(m, m, 0)$ and substitute (4.1) into (1.1) to obtain the following ODE

where

Clearly, by (1.7) we have

The matrix $\mathcal J_k$ has three eigenvalues, denoted by $\lambda_k$ ($i = 1, 2, 3$), then

where $D_i = D_i(k^2)$, $i = 0, 1, 2$, are given by

Define

If the real parts $Re(\lambda_k) < 0$ for all $k\in\mathbb N$, then the equilibrium $(m, m, 0)$ is linearly stable. It follows from (4.4) that $I_1\geq0$, which alongside (4.5) and (4.6) implies that

Therefore, by the Routh-Hurwitz criterion (cf. [56,Appendix B]) we know that $(m, m, 0)$ is linearly stable if and only if $D_0(k^2) > 0$ for all $k\in\mathbb N$. If there exist some $k_0\in \mathbb N$ such that $D_0(k_0^2) < 0$, then $(m, m, 0)$ is linearly unstable. Specifically, we have the following linear stability.

Lemma 4.1. For the system (1.1)-(1.2) with $\chi > 0$, $\gamma\geq0$, $\alpha > 1$ and $m$ given by (1.7), the constant solution $(m, m, 0)$ is linearly unstable if and only if

where $I_1$ is given by (4.3).

Proof. By (4.7) and the Routh-Hurwitz criterion, we see that $(m, m, 0)$ is linearly unstable if and only if $D_0(k^2) < 0$ for some $k\in\mathbb N^+$ (noting that $D_0(k^2)\mid_{k = 0} = I_1\geq0$). The proof is completed by noting that $D_0(k^2) < 0$ if and only if $\chi > \chi^*_k$ for $k\in\mathbb N^+$.

4.2.   Weakly nonlinear analysis

In this subsection, we assume that

with $0 < \varepsilon\ll1$ (the simpler case $\gamma = 0$ can be discussed similarly), then the system (1.1) becomes

Our goal is to perform nonlinear analysis in a weak sense under the condition $0 < \varepsilon\ll 1$. In this respect, we shall derive the amplitude equations for (4.10) by following a standard perturbation expansion approach (cf. [44,Sect. 6]), see also [50,Sect. 3.2] for an application in a symmetric two-species competition chemotaxis model).

We have from (4.3) and Lemma 4.1 that $I_1 = 1$ and hence (4.8) gives $\chi^* = \inf_{k\in\mathbb N^+}\frac{\left(\omega_k^2+1\right)^3}{\omega_k^2}$. Since $\varphi_1(s): = \frac{(s^2+1)^3}{s^2}$ is a convex function in $s > 0$, and $\min\limits_{s > 0}\varphi_1(s) = \varphi_1(\frac{\sqrt2}2) = \frac{27}{4}$. Therefore, for any $L > 0$, $\chi^* = \inf\limits_{k\in\mathbb N^+}\chi^*_k$ attains its infimum at either one or two modes. Without loss of generality, we assume there is only one mode $k^*\in\mathbb N^+$ such that $\chi^* = \chi^*_{k^*} = \inf\limits_{k\in\mathbb N^+}\chi^*_k$. Furthermore, for brevity (see Remark 4.2), we assume

then $\frac{k\pi}L > \frac{\sqrt2}2$ for all $k\in\mathbb N^+$, which indicates that $k^* = 1$ and

Moreover, the following result is a direct consequence of Lemma 4.1.

Corollary 4.2. $\mathcal J^*_k$ is invertible with three eigenvalues have negative real parts for all $k\in\mathbb N$ except at $k = 1$, where $\mathcal J^*_1$ has an eigenvalue 0 with the corresponding eigenvector

Proof. Using (4.3) with $k = 1$, $\chi = \chi^*$ and (4.9), one can obtain three eigenvalues of $\mathcal J^*_1$: 0 and $-\frac{\sqrt3 \left(1+\omega_1^2\right)\left({\sqrt3\pm\rm i}\right)}{2}$, and the corresponding eigenvector of the eigenvalue 0 is $\bm\xi$ given by (4.13). The rest of the proof follows from the proof of Lemma 4.1 immediately.

For $0 < \varepsilon\ll1$ and $\chi = \chi^*+ \varepsilon$, we define the slow time variable $\tau: = \varepsilon t$ and $(U, V, W)(x, \tau): = (u, v, w)(x, t)$, then by (4.10) we know that $(U, V, W)(x, \tau)$ satisfies

We shall use the following asymptotic expansions of $(U, V, W)$ for $0 < \varepsilon\ll1$:

The reason the expansion to be powers of $\varepsilon^\frac 12$ is that a pitchfork bifurcation can be expected although there is no symmetry of the system, see Remark 4.1 for a detailed reason.

Remark 4.1. Pitchfork bifurcations occur generally for symmetric models, such as symmetric chemotaxis models as in ^[50,52]. Here, for the asymmetric system (4.14), a pitchfork bifurcation may also occur at $\varepsilon = 0$. To see this, we denote

with $\mathcal L_1(U, W): = \partial_{xx} U- \partial_{x} (UW)+U-U^2$, and assume that for small $0 < \varepsilon\ll1$, there exists a steady state which can be appropriately expressed as $(\overline U_f, \overline V_g, \overline W_h): = (1+f(x), 1+g(x), h(x))$ for some small $f, g, h$, i.e., $\mathcal L \left({(\overline U_f, \overline V_g, \overline W_h)^T}\right) = (O(\varepsilon), O(\varepsilon), O(\varepsilon)).$ We claim that if

then

which suggests a pitchfork bifurcation. By linearity we only need to verify $\mathcal L_1(\overline U_f, \overline W_h) = \mathcal L_1(\overline U_{-f}, \overline W_{-h})$. In fact, it follows from (4.16) that

which confirms our claim. There is $(f(x), g(x), h(x))$ satisfies (4.16), for example,

with positive constants $p_j$, $j = 0, 1, 2, 3, 4$. The above analysis suggests the possible existence of a pitchfork bifurcation, and we shall see that it indeed exists (see (4.37)).

As in deriving (4.1), we have the extension

where

Noting that there is no need to specify $d_0^{(j)}(\tau)$ for all $j\in\mathbb N^+$ due to (4.2). Substituting (4.18) into (4.14), then we can collect the terms of order $\varepsilon^\frac j2$ for $j\in\mathbb N^+$ (noting that there are no terms of order $\varepsilon^0$). To begin with, we consider the case of $j = 1$.

$\bullet$ Order $\varepsilon^\frac12$ terms: collecting the terms of order $\varepsilon^\frac 12$, one has

Substituting (4.18) with $j = 1$ into (4.20), one has

where $\mathcal M_k$ is defined by (4.19) and $\mathcal J^*_k: = \mathcal J_k\mid_{\chi = \chi^*}$ in (4.3) with $\chi^*$ given by (4.12). It follows from Corollary 4.2 and the orthogonality of sine and cosine series that

For $k = 1$, $(a_1^{(1)}(\tau), b_1^{(1)}(\tau), d_1^{(1)}(\tau))^T$ is a multiple of $\bm\xi$ given in (4.13), this along with (4.2), (4.18) and (4.22) indicates that

where $P(\tau)$ is the unknown amplitude function of $\tau$. Our goal is to find the differential equation that $P(\tau)$ satisfies. To this end, we shall derive the equations that higher order $\varepsilon$ terms satisfy.

$\bullet$ Order $\varepsilon$ terms: the order $\varepsilon$ terms satisfy the following equations

As in deriving (4.21), substituting (4.18) with $j = 2$ into (4.24) and using (4.23), we have

where $(\xi_1, \xi_2, \xi_3)$ is given by (4.13) and

Again by Corollary 4.2 and the orthogonality of sine and cosine series, one has

where $Q(\tau)$ is some unknown function. For $k = 0$ and $k = 2$, (4.25) turns into the following two equations

respectively. Using (4.3) with $\chi = \chi^*$ given in (4.12), one can easily obtain $(\mathcal J^*_0)^{-1}$ and $(\mathcal J^*_2)^{-1}$, which alongside (4.2), (4.18), (4.26) and (4.27) gives

where $\bm\xi$ and $\mathcal M_k$ ($k\in\mathbb N$) are given by (4.13) and (4.19) respectively, and

$\bullet$ Order $\varepsilon^\frac32$ terms: we are now in a position to collect the terms of order $\varepsilon^\frac32$ which provides $\frac d{d\tau}P(\tau)$ and allows us to establish an ODE for $P(\tau$). As before, we have

which along with similar procedures as in deriving (4.21) and (4.25) implies

Substituting the expressions (4.23) and (4.28) for the lower oder terms into the above equation, and using $\omega_{2} = 2\omega_1$ and the orthogonality of sine and cosine series, we can collect the $\cos(\omega_1x)$ terms for the first two components of (4.30) and the $\sin(\omega_1x)$ terms for the third component of (4.30) to arrive at

where

and $\bm\xi$, $\bm\eta$ are given by (4.13) and (4.29) respectively. Let

then $\bm\beta$ is in the left kernel of $\mathcal J^*_1$. Multiplying (4.31) from left by $\bm\beta$ and using $\bm\beta\cdot\bm\xi = 3$, one has

with

where $c_3 < 0$ use the fact that $L^6-12 \pi ^4 L^2-20 \pi ^6 < 0$ due to (4.11).

By introducing $A(t) = \varepsilon^\frac12P(\tau)$, (4.32) leads to

where $c_3, c_4$ are given by (4.33). Since $c_3 < 0$ and $c_4 > 0$, we know that (4.34) has a supercritical pitchfork bifurcation at $\varepsilon = 0$ and (4.34) has three steady states:

Now plugging (4.23) into (4.15) and using $A(t) = \varepsilon^\frac12P(\tau)$, we can obtain the appropriated solution for $0 < \varepsilon\ll1$:

where $A(t)$ solves (4.34). For $0 < \varepsilon\ll1$, we have from $\lim\limits_{t\rightarrow +\infty}\tau = +\infty$ (due to $\tau = \varepsilon t$), (4.35) and (4.36) that besides the unstable equilibrium $(1, 1, 0)$, the system (4.10) also has two other stable steady states with very small amplitude in the form of

This shows that a pitchfork bifurcation indeed occurs for the asymmetric (4.10) at $\varepsilon = 0$. Moreover, it is easy to verify that (4.17) is satisfied thanks to

Remark 4.2. In this subsection, we have only discussed the case $L < \sqrt2\pi$ in (4.11) which is devoted to letting $k^* = 1$ and giving a brevity discussion. For a general $L > 0$, $\chi^*$ in (4.12) is instead by $\chi^* = \frac{\left(\omega_{k^*}^2+1\right)^3}{\omega_{k^*}^2 }$ for some $k^*\in\mathbb N^+$ (which can be specified according to the given $L > 0$). Then the rest of the discussion is similar to the above.

4.3.   Numerical simulations for the fully nonlinear regime

In this subsection, we shall perform numerical simulations to demonstrate possible patterns generated by the system (1.1)-(1.2) for large $\chi > \chi^*$. We shall consider two cases: with/without the logistic source term. To begin with, we consider the case of $\gamma = 0$.

$\bullet$ Case 1: $\gamma = 0$. With $\Omega = (0, L)$, Lemma 4.1 shows that patterns can be expected for $\chi > \chi_* = \frac{\left(1+\omega_1^2\right)^2 }m$ with $\omega_1 = \frac\pi L$. We let $L = 50$ and the initial data $(u_0, v_0, w_0)$ be a small random perturbation of the equilibrium $\left({m, m, 0}\right)$ (without loss of generality, the species mass $m$ is assumed to be 1):

where $R\in(-0.01, 0.01)$ is a random number generated by MATLAB, then we have from (4.8) that

We observe two kinds of patterns generating from the equilibrium $\left({1, 1, 0}\right)$ for $\chi > \chi_*$, as shown in Figure 1, where (a) shows that stationary patterns appear for $\chi = 2$ (not too large than $\chi^*\approx1.008$), while single merging aggregation patterns appear for $\chi = 20$. Moreover, the later kind of patterns typically appear for large $\chi$, such as $\chi = 50$, we omit the pattern here since it is similar to Figure 1 (b) (and the difference is that the larger $\chi$ is, the more aggregation of $(u, v, w)$). Similar stationary patterns as shown in Figure 1 (a) are observed for (1.1) with $\tau_0 = \gamma = 0$ in one-dimensional case (cf. ^[15]). Actually, these two kinds of patterns shown in Figure 1 have been also observed in classical chemotaxis models with/without the logistic source term such as ^[43,57,62], where the works indicate that besides the above two kinds, more patterns will appear in the existence of the logistic source term. From this point of view, we shall next consider the case $\gamma = 1$ and numerically investigate whether other types of patterns (such as periodic and chaotic) can arise for the system (1.1)-(1.2).

$\bullet$ Case 2: $\gamma > 0$. Without loss of generality, we take $\gamma = m = 1$ and $\alpha = 2$ as in (4.9). We first take $L = 50$, then we have from (4.12) that

Figure 2 shows four kinds of patterns generated from $(1, 1, 0)$ of the system (1.1)-(1.2) for $\chi = 8, 10, 12, 30$ respectively. Figure 2 (a) shows stationary patterns for $\chi = 8$, while the case of $\chi = 10$ gives single merging aggregation patterns in Figure 2 (b). These two patterns are similar as that shown in Figure 1 for the case of $\gamma = 0$. As $\chi$ goes on increasing, other kinds of patterns are observed as desired. When $\chi = 12$, Figure 2 (c1) shows the spatially inhomogeneous periodic patterns which are "almost" strip patterns (as the strip patterns shown in Figure 2 (a)) with the eventual time periodicity (at a fixed position $x = 25$) shown explicitly in Figure 2 (c2). Figure 2 (d) considers the case of $\chi = 30$ and shows chaotic patterns corresponding to merging (two local maxima join together if they are close to each other) and emerging (a new maximum form in the space between two local maxima separated from each other, cf. ^[45]). The persistent pattern formation process (for large domain) of merging and emerging as shown in Figure 2 (d) is called "coarsening process" in ^[34].

We then let $L = 20$, then (4.12) implies

In the supercritical case $\chi = 30$, Figure 3 (a) shows a spatially inhomogeneous time-periodic pattern which is more obvious than that shown in Figure 2 (c1). Although these patterns are spatially inhomogeneous, however, at a fixed time, say $t = 600$, we see that the organism $u$ prefers to cluster in places where the signal is highly concentrated, see Figure 3 (b). Furthermore, when $x$ is fixed, e.g., $x = 10$, it is clear that Figure 3 (c) shows that the solution is eventually periodic in time.

Remark 4.3. The above numerical simulations show that for the chemotaxis system (1.1)-(1.2) which is considered under acceleration assumption, large $\chi$ would not destroy the aggregation patterns (two kinds of patterns shown in Figure 1) generated from the nontrivial constant steady state, and can also not produce spatially inhomogeneous periodic patterns or chaotic patterns without the logistic source term. While the introduction of the logistic source term can largely enrich the diversity of patterns as shown in Figure 2 and Figure 3. These observations are consistent with the results on pattern formation for classical chemotaxis models with/without logistic source terms (or volume-filling effect), see ^{[43,53,54,57,62]} for instance.

4.4.   Discussion

The above numerical simulations show that a variety of patterns are observed for the system (1.1)-(1.2). The rich dynamics of solutions leave many interesting questions to be explored, below we discuss some of them.

1. Fully parabolic systems with nonclassical boundary conditions. This paper consider the fully parabolic system (1.1) (for $n\geq1$) under the boundary condition (1.2), however, the global well-posedness (for $n\leq3$) is unknown if the boundary condition instead by

This is also one of the unsolved questions proposed in ^[15], where the key is the lack of parabolic regularity with the above boundary condition.

2. Stationary problems. Theorem 1.2 asserts that the constant equilibrium $(m, m, \textbf{0})$ is exponentially asymptotic stable for small $\chi > 0$, and hence no patterns can be expected in this stable parameter regime. In one-dimensional case $\Omega = (0, L)$ with given $L > 0$, for fixed $\gamma\geq0$, $\alpha > 1$, our linear analysis has specified the critical number $\chi^*$ and patterns can possibly arise for the supercritical case $\chi > \chi^*$. The numerical simulations suggest that the system (1.1)-(1.2) has non-constant steady states whether or not $\gamma = 0$. The existence and stability of non-constant steady states is an interesting and challenging question. In this paper, we only take a very small step in this direction: still in the one-dimensional case, for the weakly nonlinear case $\chi = \chi^*+ \varepsilon$ with $0 < \varepsilon\ll1$, the solutions of the stable non-constant steady states are expressed approximately and the occurrence of pitchfork bifurcations at $\varepsilon = 0$ are confirmed for the asymmetric system (1.1)-(1.2). However, for the fully nonlinear case or for $n\geq2$, a rigorous analysis of the existence and stability of non-constant steady states is blank.

3. Eventually periodic solutions. For $\chi^*\ll\chi$ and $\gamma > 0$, the above numerical simulations indicate that the dynamics of solutions will be more complex, the spatially inhomogeneous time-periodic patterns and chaotic patterns corresponding to merging and emerging are both observed. For the former case, how to give a comprehensive mathematical explanation of the spatially inhomogeneous time-periodic patterns (at least in the one-dimensional case) is a challenging problem.

Acknowledgments

The authors are grateful to the anonymous reviewers for useful comments which improved the exposition of this paper. The research of C. Mu is partially supported by the National Natural Science Foundation of China (Nos. 11971082, 12271064), the Chongqing Talent Support program (No. cstc2022ycjh-bgzxm0169), the Natural Science Foundation of Chongqing (No. cstc2021jcyj-msxmX1051), the Fundamental Research Funds for the Central Universities (Nos. 2019CDJCYJ001, 2020CDJQY-Z001), and Chongqing Key Laboratory of Analytic Mathematics and Applications. The research of W. Tao is partially supported by the National Natural Science Foundation of China (No. 12201082), the fellowship of China Postdoctoral Science Foundation (No. 2020M683232), the Fundamental Research Funds for the Central University (No. 2021CDJQY-047) and PolyU Postdoc Matching Fund Scheme Project ID P0030816/B-Q75G and P0034904 (Primary Work Programme W15F).

Conflict of interest

The authors declare there is no conflict of interest.