Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 98, pp. 1–21. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2025.98 BLOW-UP PREVENTION AND RATE OF CONVERGENCE OF SOLUTIONS FOR N-DIMENSIONAL PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION OF CHEMOATTRACTANT JIASHAN ZHENG, YUYING WANG Abstract. This article studies the Neumann-boundary initial-value problem for a parabolic- parabolic chemotaxis-consumption system in a smooth bounded domain. For regular nonnega- tive initial data, we prove that the classical solution to the corresponding no-flux problem re- mains globally and uniformly bounded under structural assumptions. This is achieved through a novel trigonometric-type weight function rather than an exponential one; therefore we not only significantly improve previous results, but also providing a versatile context to resolve per- tinent systems. More importantly, we confirm the convergence of the solution to an equilibrium constant. 1. Introduction Across phylogenetic scales, ranging from unicellular bacteria to macroscopic mammals, or- ganismal survival hinges on the capacity to navigate the complex environment by synthesizing and interpreting multifaceted internal and external chemical signals. This chemotactic motility underpins diverse adaptive behaviors, including nutrient foraging, predator evasion, and sexual re- production. Generally, analogous signal-response migratory patterns are observed in multicellular contexts, where coordinated cell movements play a pivotal role in embryonic organogenesis and adult tissue maintenance. Drawing from the mathematical insights, chemotaxis, an evolutionary conserved biological mechanism, denotes the oriented movement of motile organisms in response to spatial and temporal gradients of chemoattractant, and has spurred an overwhelming amount of interdisciplinary research endeavors. This highly refined sensory-motor integration paradigm endows living systems with the ability to detect and react to the imperceptible fluctuations of chemoattractant concentrations, thereby exerting a critical influence on governing indispensable biological processes, which span microbial foraging behaviors to mammalian immune surveillance. One of the simplest mathematical descriptions of chemotactic cell motility is epitomized in the trailblazing minimal coupled reaction-diffusion equations initiated by Keller and Segel [13] ut = ∆u− χ∇ · (u∇v), x ∈ Ω, t > 0, vt = ∆v − v + u, x ∈ Ω, t > 0 (1.1) in a bounded domain Ω ⊂ RN (N ≥ 1) with smooth boundary, wherein the positive parameter χ, referred to as the chemotaxis sensitivity coefficient, quantifying the cellular response intensity to chemical gradients. For the unknowns u = u(x, t) and v = v(x, t) that correspond to the popula- tion density of migrating cells (or bacteria) and the concentration of chemoattractant substance, respectively, the system (1.1) itself draws attention to the bidirectional interplay between cells and chemical signal, while disregarding the surrounding environment that may influence the chemotac- tic dynamics. Over the past few decades, a growing interest of investigations has been centered on 2020 Mathematics Subject Classification. 35K20, 35K55, 92C17. Key words and phrases. Chemotaxis; blow-up prevention; chemoattractant consumption; global existence; exponential decay. ©2025. This work is licensed under a CC BY 4.0 license. Submitted August 4, 2025. Published October 16, 2025. 1 2 J. ZHENG, Y. WANG EJDE-2025/98 the corresponding problem (1.1) and its various modifications with the intention to describe the processes dominated by chemotaxis effects, resulting in solution existence, asymptotic behavior and blow-up occurrence under multiple restricting hypotheses [2, 3, 4, 6, 7, 18, 19, 20, 21, 38, 42, 44, 47]. Whenever chemotactically migrating cells fail to produce signaling molecules, but instead nav- igate by consuming chemical gradients, considerably different signal evolution paradigms emerge compared to those in (1.1), during which such scenario is biologically exemplified by oxygen tactic motility patterns displayed by swimming aerobic bacteria like Bacillus subtilis. In the context of disregarding hydrodynamic interactions with the surrounding environment, a heuristic paradigmatic [14] takes into explicit account the signal consumption by individual organisms, and constitutes a fluid-free counterpart to the coupled chemotaxis-(Navier-)Stokes system: ut = ∆u− χ∇ · (u∇v), x ∈ Ω, t > 0, vt = ∆v − uv, x ∈ Ω, t > 0. (1.2) As to (1.2), the synergistic influence of diffusion and absorptive depletion serves as an effec- tive mechanism to suppress potential blow-up phenomena in two-dimensional bounded domains with Neumann boundary conditions, which is of significant importance as it ensures the stabil- ity. Compared to the classical Keller-Segel system, the availability of a priori bounds for the initial chemical concentration v in L∞(Ω) indicates a heightened likelihood toward the tendency of global existence and boundedness, leaving unclear ambiguity how this may be used to exert control over the quantity ∇v, which directly dictates chemotactic cell movement, however. And particularly, the author in [22] revealed that the global classical solution (u, v) pertaining to the considered initial-boundary value problem is global and uniformly bounded provided that the condition 0 < ∥v0∥L∞(Ω) ≤ 1 6(N+1)χ holds, and exhibits asymptotic stabilization toward the con- stant equilibrium state ( 1 |Ω| ∫ Ω u0(x)dx, 0) as t → ∞ in the large time limit [10, 23, 40], such that the smoothly bounded classical solution for this problem belongs to C2,1(Ω̄× [T,∞)) in two- dimensional scenario under consideration, where the center to this analysis is the energy inequality d dt {∫ Ω u lnu+ ∫ Ω |∇ √ v|2 } + ∫ Ω |∇u|2 u + ∫ Ω v|D2 ln v|2 ≤ 0 for all t > 0 for smooth positive solutions in bounded convex domains. Subsequently, such finding was fur- ther generalized by Baghaei and Khelghati [1], who refined the smallness condition on the initial chemical concentration to 0 < ∥v0∥L∞(Ω) < π√ 2(N + 1)χ , while still ensuring global boundedness and convergence to the semi-trivial steady state in the regular sense of ∥v0∥L∞(Ω) < δ with some δ ∈ pi/ √ 2(N + 1)χ [16]. What is more, Heihoff [8] has constructed a global-in-time classical solution with the property 0 < ∥v0∥L∞(Ω) < 2 3Nχ , delivering a broader spectrum under which global boundedness can be assured in recent years. Other than that, extending these results to higher-dimensional version, it has been shown that the same system after all permits the existence of some globally defined solutions within weaker concepts of solvability, asymptotically stabilizing toward the associated constant steady states (u, v) → ( 1 |Ω| ∫ Ω u0(x)dx, 0 ) as t → ∞ for each appropriate initial data (u0, v0) [17, 28]. Moreover, the propensity for blow-up prevention is not restricted exclusively to (1.2). Instead, it manifests itself as a more ubiquitous trait of the chemotaxis-consumption model, which is robustly substan- tiated by a plethora of supplementary research outcomes on revealing analogous impacts in related structures, wherein the core attractant degradation mechanism from (1.2) retains intricately inter- twined with other diverse biological complexities, including the multi-species interaction dynamics [12] and fluid dynamic coupling in aqueous environments [33]. In addition to the above, it might be also pertinent to further underline that when the diffusion operator ∆u in (1.2) is modified to a more generalized form of moderately enhanced diffusion, namely ∇ · (D(u)∇u) with D(u) being assumed to satisfy D(u) → +∞ as u → ∞, the associated no-flux initial boundary problem admits globally bounded solutions even in such three dimensions [11, 37]. We encourage readers to the detailed discussions available in [5, 27, 29, 37, 39, 43, 45, 46] to receive more discoveries. EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 3 2. Main results Building upon a comprehensive synthesis of contemporary research in this direction, the evolving trajectory of studies conveys the mounting complexity and analytical exigencies inherent in higher- dimensional systems, particularly those intersecting chemotaxis phenomena and fluid dynamics, from which this progression necessitates deeper scrutiny of the fundamental mechanisms governing emergent behaviors within these frameworks. In alignment with this precedent and motivated by the above literatures, the overarching objective of the current endeavor consists in exploring both global existence and asymptotic stability of the typical chemotaxis process (1.2) itself ut = ∆u− χ∇ · (u∇v), x ∈ Ω, t > 0, vt = ∆v − uv, x ∈ Ω, t > 0 (2.1) under homogeneous Neumann boundary conditions ∂u ∂ν = ∂v ∂ν = 0, x ∈ ∂Ω, t > 0, (2.2) enforcing impermeability of the domain boundary to both bacterial cells and chemoattractant gradients, and the initial conditions u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ Ω (2.3) in a bounded domain Ω ⊂ RN with smooth boundary ∂Ω and nonnegative initial data (u0, v0), where the described parameter χ > 0 quantifies the chemotactic sensitivity coefficient, and where our main attention is devoted to the most relevant circumstance of spatial dimensions N ≥ 2. Resembling, among this type of system, u = u(x, t) is a suitable rescaled variable corresponding to the density of the cells (or organisms) population, while the unknown v = v(x, t) represents the concentration of the chemical substance. More importantly, ∂ ∂ν denotes the differentiation with respect to the outward normal vector field on the boundary ∂Ω. To the extent as we recognize, the cross-diffusive term in the first equation of (2.1) explains the adaptive motility mechanism of individual cells, wherein their movement partially implies a bias in cellular migration favor- ing regions toward increasing oxygen concentration, a behavior consistent with observed aerobic chemotactic strategies in microbial populations. Furthermore, the second term on the right-hand side of the first equation in (2.1) models the chemotactic flux component. Drawing an analogy to the well-established transport laws, such as Fourier’s principle of heat conduction, this term posits that bacterial movement in response to chemical gradients becomes proportional to the gradient magnitude itself under moderate concentration variations. Extending the rationale underlying these classical physical laws, it follows that such gradient-proportional responses are inevitable in weak gradient regimes, provided that threshold activation criteria are satisfied. Turning to the second equation in (2.1), this accounts for the supposition that oxygen is consumed at a constant rate upon bacterial contact, without compensatory production mechanisms, and that the concen- tration of substrate remains always sufficiently high so that the rate of consumption kinetic is governed by bacterial metabolic capacity, as opposed to the availability of substrate. Prior to the presentation of our main findings in these respects, the nonnegative initial data u0 and v0 are presumed to be such that u0 ∈ C0(Ω̄) with u0 ≥ 0 in Ω and u0 ̸≡ 0, x ∈ Ω̄, v0 ∈ W 1,∞(Ω) with v0 ≥ 0 in Ω and v0 ̸≡ 0, x ∈ Ω̄. (2.4) To be more precise, our first result asserting global existence and uniform boundedness reads as follows. Theorem 2.1. Let Ω ⊂ RN (N ≥ 2) be a convex bounded domain with smooth boundary ∂Ω and the parameter χ > 0. Besides, for arbitrary given initial data u0 and v0, suppose that the conditions described in (2.4) are satisfied. Then with the regular hypothesis 0 < ∥v0∥L∞(Ω) < π χ √ 2 N , (2.5) 4 J. ZHENG, Y. WANG EJDE-2025/98 it is confirmed that system (2.1)-(2.3) admits a unique nonnegative globally classical solution (u, v) with u ∈ C0(Ω̄× [0,∞)) ∩ C2,1(Ω̄× (0,∞)), v ∈ C0(Ω̄× [0,∞)) ∩ C2,1(Ω̄× (0,∞)) ∩ L∞ loc([0,∞);W 1,∞(Ω)). More importantly, this solution remains uniform-in-time bounded on Ω × (0,∞) in the situation where a constant C > 0 can be found such that ∥u(·, t)∥L∞(Ω) + ∥v(·, t)∥W 1,∞(Ω) ≤ C for all t > 0. Remark 2.2. As revealed in prior studies [8, 22], the corresponding problem (2.1) possesses globally classical solutions exclusively under the smallness of the chemotaxis sensitivity coefficient χ, whereas Theorem 2.1 demonstrates that the global existence of classical solutions for the coupled system (2.1)-(2.3) can be contingent upon the positive parameter χ within a sufficiently broader criterion. From this perspective, it is obvious that Theorem 2.1 generalizes the results of global existence when contrasted against the precedent research conjectured by [8, 22]. Remark 2.3. T o the best of our understanding, one of the most extremely substantial im- pediments encountered herein hinges on the construction of a suitable trigonometric-form weight functional φ(v) defined on 0 ≤ v ≤ ∥v0∥L∞(Ω), promoting the derivation of a priori informa- tion on Lk-bounds for u whenever k > N/2 (N ≥ 2), which absolutely distinguishes our work apart from previous investigations e.g. [8, 22, 26], where most of these studies were predominantly underpinned via exponential weight functions (as further delineated in Lemma 4.3 below). Remark 2.4. We would like to be mention that the ideas developed presented in the current con- tribution showcase unprecedented adaptability to systematically deal with a wide-ranging spec- trum of pertinent chemotactic models (cf. e.g. [32]-[36] and [43]), accordingly guaranteeing that the analogous outcomes can eventually be derived. Afterwards, our second conclusion concentrates on the exponential convergence properties of solutions to the system (2.1) under consideration. Theorem 2.5. Let Ω be a smoothly bounded domain in RN with N ≥ 2. Also, the initial data (u0, v0) are assumed to conform to the conditions (2.4). Then there exists a positive constant C such that the globally classical solution of the system (2.1)-(2.3) satisfies ∥u(·, t)− ū0∥L∞(Ω) ≤ Ce−ϑ(t−t1) for all t > t1 with ϑ < min{ū0, λ} and ū0 := 1 |Ω| ∫ Ω u0(x)dx, as well as λ representing the first nonzero eigen- value of −∆ in Ω under Neumann boundary conditions, and ∥v(·, t)∥L∞(Ω) ≤ Ce−ϱ(t−t1) for all t > t1, where ϱ < ū0 and t1 > 0 is some fixed time. Remark 2.6. While earlier seminal studies constructed solution existence for chemotaxis-consumption systems [8, 22], we complement the theoretical research context by rigorously developing the char- acterization of large-time asymptotic behavior of solutions in the present contribution, which enables our research to be more meaningful. Remark 2.7. In compliance with a series of analogous approaches from [33, 41] accompanied by slight adaptations, and as a consequence of our enhanced regularity estimates of solutions derived in our subsequent analysis, however, this implies that we are ready to streamline certain verification procedures. For the sake of notational clarity and conciseness, we set ∫ Ω fdx and ∫ t 0 ∫ Ω fdxds as ∫ Ω f and∫ t 0 ∫ Ω f , respectively, without any ambiguity. Moreover, it is worth pointing out that the values of the positive constants symbolized by C, Ci (i = 1, 2, . . . ) over the whole work may potentially EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 5 vary from one line to another, or even in the same line. Simultaneously, in numerous scenarios where it proves convenient, we shall frequently use equivalent abbreviations like, for a function f , ∥f(t)∥Lp(Ω) = ∥f(·, t)∥Lp(Ω) = (∫ Ω |f(x, t)|pdx )1/p . The remaining of this paper is organized as follows. In the present Section 2, mainly in- spired from the pioneering work [13] regarding chemical production mechanisms in the classical Keller-Segel system (1.1), and from [14] on the influence of chemoattractant dynamics in direct chemotaxis-consumption model (1.2), we pay attention to continuing to explore the system (1.2) under the initial-boundary value conditions (2.2)-(2.3). This leads to our primary contributions on both global existence of boundedness and large-time convergent behavior especially for higher- dimensional configurations N ≥ 2. Section 3 constructs the local-in-time existence of solutions and introduces essential preliminary conclusions that form the foundational groundwork to validate the analytical illustration required to support our main results. In Section 4, with the help of an absolutely disparate preference of an appropriate trigonometric-type weight function instead of the conventionally exponential ones, which is regarded as a pivotal technical advancement, we initially establish the Lk-estimates for the cellular density u (k > N/2), where when synergized with the standard Neumann heat semigroup theory under homogeneous boundary conditions, this facili- tates critical boundedness properties for the chemical concentration v in space domain W 1,q(Ω) with N < q < Nk (N−k)+ , becoming a cornerstone of Lp-bounds of the first component u of solu- tions for arbitrary large p > max{N, k} by means of the free conservation ∥u∥L1(Ω) = ∥u0∥L1(Ω). Beyond that, having at hand the above materials and benefiting from the maximum principle alongside the standard parabolic regularity theory or the Lp-Lq estimates for the Neumann heat semigroup (et∆)t≥0 as well as the extensibility criterion of the local existence of solutions, we consequentially authenticate the statements from Theorem 2.1. Finally, in terms of the theoret- ical reasoning as revealed in [33] and [41], Section 5 derives the exponential convergence decay behavior to equilibrium states of solutions, and further, certify the assertion of Theorem 2.5. 3. Preliminaries and local well-posedness As a preparation toward demonstrating the qualitative identities of the classical solution, we originally concentrate on the local solvability and uniqueness of solutions to the system (2.1)-(2.3) drawn upon the well-established fixed-point arguments of the interrelated framework deduced in [25, Lemma 2.3] and [31, Lemma 2.1], meaning that we shall leave out the concrete proof to prevent duplication. Lemma 3.1 (Local existence and uniqueness). Suppose that Ω ⊂ RN (N ≥ 2) is a smoothly bounded domain. Then for arbitrary given initial datum (u0, v0), there exist the maximal existence time Tmax ∈ (0,∞] and a uniquely determined quadruple (u, v) of nonnegative functions, namely, u ∈ C0(Ω̄× [0, Tmax)) ∩ C2,1(Ω̄× (0, Tmax)), v ∈ C0(Ω̄× [0, Tmax)) ∩ C2,1(Ω̄× (0, Tmax)) ∩ L∞ loc([0, Tmax);W 1,∞(Ω)), which solves (2.1)-(2.3) in the classical sense in Ω × (0, Tmax), and such that u, v > 0 in Ω̄ × (0, Tmax). In particular, if the maximal existence time Tmax < ∞, then lim sup t↗Tmax { ∥u(·, t)∥L∞(Ω) + ∥v(·, t)∥W 1,∞(Ω) } = ∞. Before going further, it is imperative to cite the findings related to the asymptotic behavior of the heat semigroup under Neumann boundary conditions, which are crucial for some estimates referred later. Lemma 3.2 ([30, Lemma 1.3]). Let (eτ∆)τ≥0 be the Neumann heat semigroup in Ω, and let µ denote the first nonzero eigenvalue of −∆ in Ω under Neumann boundary conditions. Then there exist constants ki = ki(Ω)(i = 1, 2, 3, 4) depending on Ω only that possess the following properties. 6 J. ZHENG, Y. WANG EJDE-2025/98 (i) Whenever 1 ≤ q ≤ p ≤ ∞, we have ∥eτ∆φ∥Lp(Ω) ≤ k1 ( 1 + τ− N 2 ( 1 q− 1 p ) ) e−µτ∥φ∥Lq(Ω) (3.1) for all τ > 0 for each φ ∈ Lq(Ω) is valid with ∫ Ω φ = 0. (ii) In the setting 1 ≤ q ≤ p ≤ ∞, we have ∥∇eτ∆φ∥Lp(Ω) ≤ k2 ( 1 + τ− 1 2− N 2 ( 1 q− 1 p ) ) e−µτ∥φ∥Lq(Ω) (3.2) for all τ > 0 and any φ ∈ Lq(Ω); (iii) If 1 < q ≤ p < ∞, then ∥eτ∆∇ · φ∥Lp(Ω) ≤ k3 ( 1 + τ− 1 2− N 2 ( 1 q− 1 p ) ) e−µτ∥φ∥Lq(Ω) (3.3) for all τ > 0 and every φ ∈ (C∞ 0 (Ω))N Additionally, the operator eτ∆∇· admits a uniquely determined improvement to an operator from Lq(Ω) into Lp(Ω), with the norm controlled according to (3.3). (iv) When 2 ≤ p < ∞, we have ∥∇eτ∆φ∥Lp(Ω) ≤ k4e −µτ∥∇φ∥Lp(Ω) for all τ > 0 and φ ∈ W 1,p(Ω). (3.4) Apart from that, we list some inequalities frequently applied in the forthcoming statements. Actually, during the proof of the main results, we shall make use of the commonly recognized Young inequality with ε. Lemma 3.3. Assume p and q are given positive numbers obeying 1 < p, q < +∞ and 1 p + 1 q = 1. Then for any ε > 0 and positive constants a and b, it follows that ab ≤ εap + 1 q (εp)− q p bq. Let us then invoke the Gagliardo-Nirenberg interpolation inequality that is employed consis- tently throughout this paper. Lemma 3.4. Let Ω ⊂ RN (N ≥ 1) be a smoothly bounded domain. In addition, suppose p ≥ 1 and 0 < q ≤ p. Then one can determine a constant CGN = C(p, q,N,Ω) > 0 such that for arbitrary r > 0 and each φ ∈ W 1,2(Ω) ∩ Lq(Ω), ∥φ∥Lp(Ω) ≤ CGN ( ∥∇φ∥θL2(Ω)∥φ∥ 1−θ Lq(Ω) + ∥φ∥Lr(Ω) ) , where θ ∈ (0, 1) is given by 1 p = (1 2 − 1 N ) θ + 1− θ q ⇐⇒ θ = N q − N p 1− N 2 + N q . 4. A priori estimates. Proof of Theorem 2.1 This section is dedicated to elaborating on a priori estimates of the local solution as the starting point toward its extension to be a global-in-time one, whence our reasoning consequentially is based on an absolutely different technique, which at its core challenge relies in the Lk-norm regularity of the population density u for k > N 2 and N ≥ 2, and which shall be accomplished through meticulously tracking the evolution undergone by coupled functional of the form∫ Ω ukφ(v), (4.1) coupled with suitable chosen weight function φ(v) contingent upon the chemoattractant concen- tration v, designed to remain uniformly bounded from both above and below by positive constants. In stark contrast to the preceding relevant studies, where similar tactics have been explored in prior works (cf. [8, 22, 26, 29]), the concurrent study introduces a methodological innovation necessitated by the chemoattractant consumption mechanism. Especially, we develop a novel trigonometric-type weight function that means a refinement over exponential one to deal with the structural complexities arising from consumption terms. To achieve this, we henceforth separately EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 7 conduct an analysis in the following rudimentary but pivotal lemmas, which collectively make up the bedrock of our approach. Lemma 4.1. Let (u, v) be the solution to (2.1)-(2.3) and (2.4) be valid. Under the mild as- sumptions of Lemma 3.1, the first component u of the solution satisfies the mass conservation property ∥u(·, t)∥L1(Ω) = ∥u0∥L1(Ω) for all t ∈ (0, Tmax). (4.2) Moreover, one has ∥v(·, t)∥L∞(Ω) ≤ ∥v0∥L∞(Ω) for all t ∈ (0, Tmax). (4.3) Proof. The assertion (4.2) yields an instantaneous conclusion of the integration for the first equa- tion in (2.1), whereas (4.3) results from the standard parabolic maximum principle. □ In the sequel, we implicitly consider a quadruple (u, v) of nonnegative functions obtained in Lemma 3.1 as the solution of equations (2.1) that is described by the maximal existence time Tmax ∈ (0,∞]. In the aftermath of the analytical foundation specified in Lemma 4.1, we are now well-positioned to compile a hierarchy of crucial properties concerning the selected weight functional φ(v) defined on 0 ≤ v ≤ ∥v0∥L∞(Ω) thereof, consequentially paving the way for higher regularity estimates of solutions. Lemma 4.2. Suppose that Ω is a bounded domain with smooth boundary in RN and N ≥ 2. Furthermore, let Theorem 2.1 be valid. Then one introduces a weight function φ(v) := { cos (√kχ 2 v )}−(k−1) for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.4) such that φ(v) > 0 for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.5) φ′(v) > 0 for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.6) φ′′(v) > 0 for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.7) with k > N/2 and sufficiently approaching to N/2. Proof. At first, whenever 0 ≤ v ≤ ∥v0∥L∞(Ω), we plan to demonstrate that 0 ≤ √ k 2 χv ≤ π 2 . Through taking advantage of (2.5), this allows for a choice of a certain constant L > 1 such that ∥v0∥L∞(Ω) := π Lχ √ 2 N < π χ √ 2 N , whence for each δ ∈ (1, L], letting k := N 2 δ 2, and a series of straightforward calculations immediately reveals that 0 ≤ √ k 2 χv ≤ √ k 2 χ∥v0∥L∞(Ω) = δ 2L π ≤ π 2 . So that as observed in (4.4), we abbreviate for completeness and convenience that z(v) := lnφ(v) = −(k − 1) ln { cos (√kχ 2 v )} for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.8) which in turn enables us to quickly discover φ′(v) = φ(v)z′(v) for all 0 < v ≤ ∥v0∥L∞(Ω) (4.9) as well as φ′′(v) = φ(v)[ ( z′(v) )2 + z′′(v)] for all 0 < v ≤ ∥v0∥L∞(Ω). (4.10) Followed by the weight function and the supposition provided in (4.4) and (2.5) respectively, this readily entails (4.5). As to the claim of (4.6)-(4.7), the identity of (4.8) applies so as to warrant z′(v) = √ k(k − 1)χ 2 tan (√kχ 2 v ) > 0 for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.11) because tan( √ kχ 2 v) > 0, by (2.5), and z′′(v) = k(k − 1)χ2 4 sec2 (√kχ 2 v ) > 0 for all 0 ≤ v ≤ ∥v0∥L∞(Ω), (4.12) 8 J. ZHENG, Y. WANG EJDE-2025/98 whereupon taking into consideration this together with (4.5), (4.9)-(4.11) trivially contributes to not only (4.6) but also (4.7), and moreover, correspondingly completes the proof of Lemma 4.2. □ Bases on Lemma 4.2, it suffices to affirm the Lk-norm bounds of the solution u(·, t) achieved via the functional (4.1) in the scenario k > N 2 and N ≥ 2 for all t ∈ (0, Tmax). Lemma 4.3. Let Theorem 2.1 hold. Upon theabove-mentioned hypotheses of Lemma 4.2, there exists a suitable positive constant C fulfilling ∥u(·, t)∥Lk(Ω) ≤ C for all t ∈ (0, Tmax), (4.13) provided that k > N/2 is adequately close to N/2. Proof. Going back to the first and second equations of (2.1), performing the simple differentiation and integrating it by parts across the domain Ω, one has 1 k d dt ∫ Ω ukφ(v) = ∫ Ω uk−1φ(v)ut + 1 k ∫ Ω ukφ′(v)vt = ∫ Ω uk−1φ(v) [∆u− χ∇ · (u∇v)] + 1 k ∫ Ω ukφ′(v) (∆v − uv) = − ∫ Ω [ (k − 1)uk−2φ(v)∇u+ uk−1φ′(v)∇v ] ∇u+ χ ∫ Ω [ (k − 1)uk−2φ(v)∇u+ uk−1φ′(v)∇v ] u∇v − 1 k ∫ Ω ( kuk−1φ′(v)∇u+ ukφ′′(v)∇v ) ∇v − 1 k ∫ Ω uk+1vφ′(v) = −(k − 1) ∫ Ω uk−2φ(v)|∇u|2 − ∫ Ω ( 1 k φ′′(v)− χφ′(v) ) uk|∇v|2 + ∫ Ω [(k − 1)χφ(v)− 2φ′(v)]uk−1∇u · ∇v − 1 k ∫ Ω uk+1vφ′(v) for all t ∈ (0, Tmax). Noticing the nonnegativity of the solutions to the associated problem (2.1) and (4.6) then reveals 1 k ∫ Ω uk+1vφ′(v) > 0 for all t ∈ (0, Tmax), that is, −1 k ∫ Ω uk+1vφ′(v) < 0 for all t ∈ (0, Tmax), therefore yielding 1 k d dt ∫ Ω ukφ(v) ≤ −(k − 1) ∫ Ω uk−2φ(v)|∇u|2 − ∫ Ω ( 1 k φ′′(v)− χφ′(v) ) uk|∇v|2 + ∫ Ω [(k − 1)χφ(v)− 2φ′(v)]uk−1∇u · ∇v := J1 + J2 + J3 for all t ∈ (0, Tmax). (4.14) On the basis of (4.11)-(4.12), we interpolate some elementary calculations to certainly achieve z′′(v) = k(k − 1)χ2 4 sec2 (√kχ 2 v ) = k(k − 1)χ2 4 { 1 + tan2 (√kχ 2 v )} = 1 k − 1 (z′(v)) 2 + 1 4 k(k − 1)χ2 for all 0 ≤ v ≤ M0, (4.15) whereas the identities (4.9)-(4.10) make sure that through an appropriate combination of (4.14), the second term becomes J2 := − ∫ Ω (1 k φ′′(v)− χφ′(v) ) uk|∇v|2 EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 9 = − ∫ Ω {1 k [ (z′(v)) 2 + z′′(v) ] − χz′(v) } ukφ(v)|∇v|2 for all t ∈ (0, Tmax), and by means of (4.15), the above equation gives rise to J2 = − ∫ Ω {1 k [ (z′(v)) 2 + 1 k − 1 (z′(v)) 2 + 1 4 k(k − 1)χ2 ] − χz′(v) } ukφ(v)|∇v|2 = − ∫ Ω { 1 k − 1 (z′(v)) 2 + 1 4 (k − 1)χ2 − χz′(v) } ukφ(v)|∇v|2 = − 1 4(k − 1) ∫ Ω [2z′(v)− (k − 1)χ]2ukφ(v)|∇v|2 for all t ∈ (0, Tmax). (4.16) From this, depending on another application of the basic inequality a2 + b2 ≥ 2|a| · |b| for all a, b ∈ R, we acquire from the first term amalgamating with the second integral on the right-hand side of (4.14) as well as (4.16) that J1 + J2 := −(k − 1) ∫ Ω uk−2φ(v)|∇u|2 − ∫ Ω (1 k φ′′(v)− χφ′(v) ) uk|∇v|2 = −(k − 1) ∫ Ω uk−2φ(v)|∇u|2 − 1 4(k − 1) ∫ Ω [2z′(v)− (k − 1)χ] 2 ukφ(v)|∇v|2 ≤ − ∫ Ω |2z′(v)− (k − 1)χ|uk−1φ(v)|∇u||∇v| for all t ∈ (0, Tmax), (4.17) whence again by (4.9), the third term of (4.14) assuredly contributes to J3 := ∫ Ω [(k − 1)χφ(v)− 2φ′(v)]uk−1∇u · ∇v = ∫ Ω [(k − 1)χ− 2z′(v)]uk−1φ(v)∇u · ∇v for all t ∈ (0, Tmax), (4.18) which in conjunction with (4.17)-(4.18) inserted into (4.14) occurs 1 k d dt ∫ Ω ukφ(v) ≤ 0 for all t ∈ (0, Tmax), and meanwhile, accompanied by an integration from 0 to t infers that for a suitable positive constant C1 := ∫ Ω uk 0φ(v0), ∫ Ω ukφ(v) ≤ C1 for all t ∈ (0, Tmax), henceforth yielding that there exists C2 > 0 complying with∫ Ω uk ≤ ∫ Ω ukφ(v) ≤ C2 for all t ∈ (0, Tmax) (4.19) according to the boundedness of φ(v), 1 ≤ φ(v) ≤ { cos (√kχ 2 ∥v0∥L∞(Ω) )}−(k−1) := C3 > 1 for all 0 ≤ v ≤ ∥v0∥L∞(Ω) since { cos (√kχ 2 ∥v0∥L∞(Ω) )}k−1 < 1, culminating in the confirmation of (4.13). □ Building upon Lemma 4.3 as a baseline, we pay attention to the forthcoming auxiliary lemma, which underscores to be instrumental in the further estimates of the solution (u, v) to the coupled system (2.1)-(2.3). 10 J. ZHENG, Y. WANG EJDE-2025/98 Lemma 4.4. Let Lemma 4.3 hold. One has a constant C > 0 such that the second component of solutions to system (2.1)-(2.3) satisfies ∥v(·, t)∥W 1,q(Ω) ≤ C for all t ∈ (0, Tmax), (4.20) where N < q < Nk (N − k)+ , (4.21) and k is taken from Lemma 4.3. Proof. Given that the condition k > N 2 warrants Nk (N−k)+ > N , the variation-of-constants formula imposed to the v-equation in (2.1) becomes applicable to allow v(·, t) = et∆v0 − ∫ t 0 e(t−s)∆u(·, s)v(·, s)ds for all t ∈ (0, Tmax). As a result of the known regularization properties of (et∆)t≥0 in Lemma 3.2, this allows for the choice of two positive constants C1 and C2 complying with ∥v(·, t)∥W 1,q(Ω) ≤ ∥∇et∆v0∥Lq(Ω) + ∫ t 0 ∥∇e(t−s)∆u(·, s)v(·, s)∥Lq(Ω)ds ≤ C1 + C2 ∫ t 0 ( 1 + (t− s)− 1 2− N 2 ( 1 k− 1 q ) ) e−λ1(t−s)∥u(·, s)v(·, s)∥Lk(Ω) (4.22) for all t ∈ (0, Tmax), where λ1 is the first positive eigenvalue of −∆ under homogeneous Neumann boundary conditions. From (4.21)and the identity −1 2 − N 2 (1 k − 1 q ) > −1 (4.23) we obtain ∫ t 0 ( 1 + (t− s)− 1 2− N 2 ( 1 k− 1 q ) ) e−λ1(t−s)ds ≤ ∫ ∞ 0 ( 1 + σ− 1 2− N 2 ( 1 k− 1 q ) ) e−λ1σdσ < +∞ for all t ∈ (0, Tmax), (4.24) and that for some C3 > 0, the latter on the right-hand side of (4.22) bringing together (4.23)-(4.24) also trivially provides C2 ∫ t 0 ( 1 + (t− s)− 1 2− N 2 ( 1 k− 1 q ) ) e−λ1(t−s)∥u(·, s)v(·, s)∥Lk(Ω) ≤ C3∥u(·, s)∥Lk(Ω)∥v(·, s)∥L∞(Ω) for all t ∈ (0, Tmax), (4.25) which in light of the standard parabolic maximum principle (4.3) and Lemma 4.3 yields the existence of a positive constant C4 such that C2 ∫ t 0 ( 1 + (t− s)− 1 2− N 2 ( 1 k− 1 q ) ) e−λ1(t−s)∥u(·, s)v(·, s)∥Lk(Ω) ≤ C4 for all t ∈ (0, Tmax), (4.26) whence, taking into account a substitution of (4.26) into (4.22), we arrive at the desired assertion (4.20). □ In accordance with tho bounds stipulated by Lemma 4.3 and Lemma 4.4, we next check a vital estimate of the population density u in the spaces Lp(Ω), for p > max{N, k} being taken arbitrarily large, which is an essential component for demonstrating higher regularity of solutions and in turn ultimately underpins the global existence arguments. Lemma 4.5. Assume Theorem 2.1 and Lemmas 4.3-4.4 are valid. Also assume that (u, v) is a nonnegative solution of problem (2.1), and the initial data u0 as well as v0 fulfill (2.4). Then for any choice of p > max{N, k} and k > N/2 originated from Lemma 4.3, there exists a constant C > 0 such that ∥u(·, t)∥Lp(Ω) ≤ C for all t ∈ (0, Tmax). (4.27) EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 11 Proof. Via applying the homogeneous Neumann boundary conditions, we multiply both sides of the first equation in (2.1) by up−1 for p > max{N, k}, and then infer upon the differentiation and integration by parts over the general domain Ω, that 1 p d dt ∫ Ω up = ∫ Ω up−1ut = ∫ Ω up−1 [∆u− χ∇ · (u∇v)] = −(p− 1) ∫ Ω up−2|∇u|2 + (p− 1)χ ∫ Ω up−1∇u · ∇v for all t ∈ (0, Tmax), (4.28) where we interpolate using the Young inequality to the second integral of (4.28) to obtain 1 p d dt ∫ Ω up ≤ −p− 1 2 ∫ Ω up−2|∇u|2 + 1 2 (p− 1)χ2 ∫ Ω up|∇v|2 for all t ∈ (0, Tmax). (4.29) After that, we devote our attention to absorbing each term on the right-hand side of (4.29). For the first one, under a series of simple calculations, it shows that −p− 1 2 ∫ Ω up−2|∇u|2 = −2(p− 1) p2 ∥∇up/2∥2L2(Ω) for all t ∈ (0, Tmax). (4.30) Regarding to the rightmost integral, in light of the Hölder inequality and Lemma 4.4, there is a suitable positive constant C1 such that 1 2 (p− 1)χ2 ∫ Ω up|∇v|2 ≤ 1 2 (p− 1)χ2 (∫ Ω u pq q−2 ) q−2 q (∫ Ω |∇v|q )2/q ≤ C1∥up/2∥2 L 2q q−2 (Ω) for all t ∈ (0, Tmax), from which it follows from the well-known Gagliardo-Nirenberg inequality that one introduces C2 > 0 satisfying 1 2 (p− 1)χ2 ∫ Ω up|∇v|2 ≤ C2 { ∥∇up/2∥ Np 2 −N(q−2) 2q 1−N 2 + Np 2 L2(Ω) ∥up/2∥ 1− Np 2 −N(q−2) 2q 1−N 2 + Np 2 L 2 p (Ω) + ∥up/2∥ L 2 p (Ω) }2 for all t ∈ (0, Tmax), (4.31) whence applying the elementary inequality (a + b)2 ≤ 2(a2 + b2) for all a, b ∈ R, some suitable constant C3 > 0 can be picked such that (4.31) transforms into 1 2 (p− 1)χ2 ∫ Ω up|∇v|2 ≤ C3 { ∥∇up/2∥ 2· Np 2 −N(q−2) 2q 1−N 2 + Np 2 L2(Ω) ∥up/2∥ 2−2· Np 2 −N(q−2) 2q 1−N 2 + Np 2 L 2 p (Ω) + ∥up/2∥2 L 2 p (Ω) } ≤ C3 { ∥∇up/2∥ 2· Np 2 −N(q−2) 2q 1−N 2 + Np 2 L2(Ω) ∥u0∥ p−p· Np 2 −N(q−2) 2q 1−N 2 + Np 2 L1(Ω) + ∥u0∥pL1(Ω) } for all t ∈ (0, Tmax), (4.32) and again, by means of the Young inequality and the mass conservation (4.2), it is not difficult to ascertain the existence of a positive constant C4 owing to ∥u0∥L1(Ω) > 0 complying with 1 2 (p− 1)χ2 ∫ Ω up|∇v|2 ≤ p− 1 p2 ∥∇up/2∥2L2(Ω) + C4 for all t ∈ (0, Tmax) (4.33) as a result of the restriction, (4.21) implies 2 · Np 2 − N(q−2) 2q 1− N 2 + Np 2 < 2. 12 J. ZHENG, Y. WANG EJDE-2025/98 Substituting (4.30) and (4.33) back into (4.29) correspondingly contributes to 1 p d dt ∫ Ω up ≤ −p− 1 p2 ∫ Ω |∇up/2|2 + C4 for all t ∈ (0, Tmax), (4.34) where adding ∫ Ω up for all p > max{N, k} and k > N 2 (N ≥ 2) on both sides of (4.34) indicates 1 p d dt ∫ Ω up + ∫ Ω up ≤ −p− 1 p2 ∫ Ω |∇up/2|2 + ∫ Ω up + C4 for all t ∈ (0, Tmax). (4.35) Apart from that, by the analogous procedure as in the arguments of (4.31)-(4.32), and once more recalling the Gagliardo-Nirenberg inequality, this means that there exists certain positive constants Ci(i = 5, 6, 7) such that∫ Ω up = ∥up/2∥2L2(Ω) ≤ C5 { ∥∇up/2∥ 2· Np 2 −N 2 1−N 2 + Np 2 L2(Ω) ∥up/2∥ 2−2· Np 2 −N 2 1−N 2 + Np 2 L 2 p (Ω) + ∥up/2∥2 L 2 p (Ω) } for all t ∈ (0, Tmax), (4.36) which on account of 2 · Np 2 − N 2 1− N 2 + Np 2 < 2 and in view of the Young inequality yields ∥∇up/2∥ 2· Np 2 −N 2 1−N 2 + Np 2 L2(Ω) ≤ C6∥∇up/2∥2L2(Ω) + C7 for all t ∈ (0, Tmax), (4.37) and which whenever connected with (4.37) inserted into (4.36) yields that for some C8 > 0,∫ Ω up ≤ p− 1 p2 ∥∇up/2∥2L2(Ω) + C8 for all t ∈ (0, Tmax). The aforementioned property bringing together the substitution of (4.35) trivially develops into 1 p d dt ∫ Ω up + ∫ Ω up ≤ C9 for all t ∈ (0, Tmax). Thanks to the standard ODE comparison, the affirmation of Lemma 4.5 is obtained. □ Taking advantage of the preliminary information and revisiting the second equation of (2.1), we improve our knowledge on constructing the W 1,∞-boundedness for the second component v of solutions to the corresponding question (2.1)-(2.3). Lemma 4.6. In the regular hypothesis of Theorem 2.1, it holds that we find a suitable constant C > 0 such that ∥v(·, t)∥W 1,∞(Ω) ≤ C for all t ∈ (0, Tmax). (4.38) Proof. Because ∥u(·, t)∥Lp(Ω) is bounded for any large p > max{N, k} with k > N 2 , we can obtain, based on the asymptotic behavior of the heat semigroup under Neumann boundary conditions (similar to Lemma 4.4) or the standard regularity theory of parabolic equations (refer to [15],) that (4.38) is obvious. □ With the above technical preparation at hand, we intend to use results in Lemmas 4.5-4.6 to introduce an upper limit for the solution, which is related to u within the space L∞(Ω). Lemma 4.7. Let (u, v) be the solution of system (2.1)-(2.3). Also, the initial data u0 and v0 be such that the conditions (2.4) hold. Then there is a constant C > 0 fulfilling ∥u(·, t)∥L∞(Ω) ≤ C for all t ∈ (0, Tmax). (4.39) Proof. From the boundedness of u(·, t) in the space Lp(Ω) for arbitrary large p > max{N, k} estimated in Lemma 4.5, a straightforward application of the well-known Moser-Alikakos iteration procedure (see [24, Lemma A.1]) generates (4.39), therefore verifying Lemma 4.7. □ EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 13 Now the proof of our first theorem proceeds seamlessly with the compilation of all these prepara- tory steps. 4.1. Proof of Theorem 2.1. As observed in the extensibility criterion of the local existence in Lemma 3.1, it shows that the maximal existence time Tmax = ∞. Correspondingly, Theorem 2.1 follows from connecting with the conclusions through Lemmas 4.6 and 4.7. In addition, depending on the verified Theorem 2.1, we improve our knowledge to summarize the subsequent corollary as an important analytical instrument to demonstrate the asymptotic stability hereafter. Corollary 4.8. Under the assumptions of Lemma 4.5 and Tmax = ∞, one can determine a positive constant C such that ∫ T 0 ∫ Ω |∇u|2 ≤ C for all T > 0. (4.40) Proof. Given (4.34), then the desired conclusion results upon choosing p = 2 and an integration from 0 to T with respect to t. □ 5. Exponential decay. Proof of Theorem 2.5 Having resolved the uniform boundedness concerns so far, we then devote ourselves to describing the large-time asymptotic behavior, with particular emphasis on demonstrating the exponential decay properties of solutions in this section. First of all, encouraged by the framework outlined in [33], we plan to systematically certify that the first solution component u stabilizes to the spatially uniform equilibrium state, namely, u(·, t) → ū0 given that ū0 := 1 |Ω| ∫ Ω u0(x)dx, which starts with a sequence of essential lemmas. Lemma 5.1. Suppose that Ω ⊂ RN (N ≥ 2) is a smoothly bounded domain. Then the first component u of solutions to system (2.1)-(2.3) fulfills (u(·, t))t>3 is relatively compact in C0(Ω̄). Proof. Rewriting the second equation in (2.1) gives rise to ut = ∆u− u+ u− χ∇ · (u∇v), whence according to an associated variation-of-constants formula we have u(·, t) = e−tAu(·, s0) + ∫ t s0 e−(t−s)Au(·, s)ds− χ ∫ t s0 e−(t−s)A∇ · (u(·, s)∇v(·, s))ds (5.1) for all t > s0, with s0 ∈ {2, 3, . . . } and Ameasuring the sectorial extension of−∆+1 in Lp(Ω) under homogeneous Neumann boundary conditions for arbitrary large p > max{N, k}. Furthermore, one then employs Aα under the restriction of 0 < α < 1/2 on (5.1) to indicate that ∥Aαu(·, t)∥Lp(Ω) ≤ ∥Aαe−tAu(·, s0)∥Lp(Ω) + ∫ t s0 ∥Aαe−(t−s)Au(·, s)∥Lp(Ω)ds + χ ∫ t s0 ∥Aαe−(t−s)A∇ · (u(·, s)∇v(·, s))∥Lp(Ω)ds for all t > s0. (5.2) Moving forward, we are in a position to study each term on the right-hand side of (5.2). First of all, by the regularity estimates ([9]) alongside the Lp-boundedness constructed in Lemma 4.5 and the mass conservation property (4.2), there are suitable constants λ2, Ci(i = 1, 2, 3, 4) > 0 such that ∥Aαe−tAu(·, s0)∥Lp(Ω) ≤ C1t −αe−λ2t∥u(·, s0)∥Lp(Ω) ≤ C2 for all t > s0 (5.3) and ∫ t s0 ∥Aαe−(t−s)Au(·, s)∥Lp(Ω)ds ≤ C3 ∫ t s0 (t− s)−αe−λ2(t−s)∥u(·, s)∥Lp(Ω)ds ≤ C4 (5.4) 14 J. ZHENG, Y. WANG EJDE-2025/98 for all t > s0, from which once more recalling Lemma 4.6 and [9, Lemma 2.1] ascertain the existence of certain positive constants C5 and C6 such tat χ ∫ t s0 ∥Aαe−(t−s)A∇ · (u(·, s)∇v(·, s))∥Lp(Ω)ds ≤ C5 ∫ t s0 (t− s)− 1 2−αe−λ2(t−s)∥u(·, s)∇v(·, s)∥Lp(Ω)ds ≤ C5 ∫ t s0 (t− s)− 1 2−αe−λ2(t−s)∥u(·, s)∥Lp(Ω)∥∇v(·, s)∥L∞(Ω)ds ≤ C6 for all t > s0 (5.5) since 0 < α < 1/2, unswervingly contributing to − 1 2 − α > −1 and then,∫ t s0 (t− s)− 1 2−αe−λ2(t−s)ds ≤ ∫ ∞ 0 σ− 1 2−αe−λ2σdσ < +∞. Putting the above three estimates (5.3)-(5.5) together, we obtain a constant C7 > 0 that satisfies ∥Aαu(·, t)∥Lp(Ω) ≤ C7 for all t > s0. (5.6) As a result of p > max{N, k} for k > N 2 and noticing that α ∈ (0, 1 2 ) by assumption, this particularly enables us to choose α ∈ (0, 1 2 ) fulfilling α > N 2p , and guarantee that one may introduce 0 < δ < 2α− N p such that D(Aα) ↪→ Cδ(Ω̄), and such that ∥u(·, t)∥Cδ(Ω̄) ≤ C8∥Aαu(·, t)∥Lp(Ω) ≤ C9 for all t ≥ 3 for positive constants C8 and C9. From this and (5.6), using the Arzelá-Ascoli theorem easily implies the statement of Lemma 5.1. □ Based on the Lemma 4.5 to Lemma 4.6, let us then embark on deducing the pivotal estimation below, which is a cornerstone for time evolution dynamics. Lemma 5.2. Under assumptions of Theorem 2.1, there exists a positive constant C such that∫ T 0 ∫ Ω |u∇v|2 ≤ C for all T > 0. (5.7) Proof. By considering Lemmas 4.5-4.6, combined with the Hölder inequality, we can obtain C1 > 0 such that ∫ T 0 ∫ Ω |u∇v|2 ≤ (∫ T 0 ∫ Ω u4 )1/2(∫ T 0 ∫ Ω |∇v|4 )1/2 ≤ C1 for all T > 0, arriving to (5.7). □ In the sequel, for the purpose of constructing asymptotic convergence in the full time horizon t → ∞, we focus on an elementary but relatively weak decay estimate of ut for the time derivatives. Lemma 5.3. Let Ω ⊂ RN (N ≥ 2) be a smoothly bounded domain and Lemma 5.2 hold. Then there exists C > 0 such that∫ T 1 ∥∂tu(·, t)∥2(W 1,2(Ω))∗dt ≤ C for all T > 1. Proof. Multiplying both sides of the first equation in (1.2) by φ ∈ W 1,2(Ω), integrating by parts over the domain and performing some direct calculations, it follows from the Neumann boundary conditions and the Hölder inequality that∣∣ ∫ Ω ∂tu(·, t)φ ∣∣ = ∣∣ ∫ Ω [∆u− χ∇ · (u∇v)]φ ∣∣ = ∣∣− ∫ Ω ∇u · ∇φ+ χ ∫ Ω u∇v · ∇φ ∣∣ ≤ ( ∥∇u∥L2(Ω) + χ∥u∇v∥L2(Ω) ) ∥∇φ∥L2(Ω) for all t > 1. EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 15 An application of the Young inequality, Corollary 4.8, the identity φ ∈ W 1,2(Ω) and Lemma 5.2 allows for a selection of positive constants Ci(i ∈ {1, 2}) such that∫ T 1 ∥∂tu(·, t)∥2(W 1,2(Ω))∗dt ≤ C1 {∫ T 1 ∫ Ω |∇u|2 + ∫ T 1 ∫ Ω |u∇v|2 } ≤ C2 for all T > 1, and the proof of Lemma 5.3 is complete. □ According to the preceding lemmas, the quantitative information encapsulated in (4.40) can be transformed into L∞-norm stabilization for the first component u of solutions. Specifically, proceeding a similar way as in [33, Lemma 8.2], we have that the first component component u of the solution converges uniformly to the equilibrium state ū0; thereby completing the preparation of the exponential convergence decay. Lemma 5.4. Let (u, v) be a solution of system (2.1). Then the first component u of solutions satisfies ∥u(·, t)− ū0∥L∞(Ω) → 0 as t → ∞, (5.8) where ū0 := 1 |Ω| ∫ Ω u0(x)dx. Proof. By Lemma 5.1, we just need to make sure that the initial data u0 constitutes the sole element of the corresponding ω-limit set on u. In other words, it is sufficient to show that whenever (tk)k∈N ⊂ (3,∞) and u∞ ∈ C0(Ω̄) are such that tk → ∞ and u(·, tk) → u∞ in C0(Ω̄) as k → ∞, (5.9) we necessarily confirm u∞ ≡ ū0. At this juncture, given any (tk)k∈N and u∞, we define uk(x, s) := u(x, tk + s), x ∈ Ω, s ∈ (0, 1), k ∈ N. Then let the positive constant CP > 0 represent the Poincaré constant that satisfies∥∥φ− 1 |Ω| ∫ Ω φ ∥∥2 L∞(Ω) ≤ CP ∥∇φ∥2L2(Ω) for all φ ∈ W 1,2(Ω), which upon relying on routine manipulations connected with Corollary 4.8 and (4.2) evidently leads to ∥uk − ū0∥2L2(Ω×(0,1)) = ∫ tk+1 tk ∫ Ω |u(x, t)− ū0|2 dx dt ≤ CP ∫ tk+1 tk ∫ Ω |∇u(x, t)|2 dx dt → 0 as k → ∞. (5.10) Apart form that, we write ũ∞(x, s) := u∞(x) for all (x, s) ∈ Ω× (0, 1), (5.11) and henceforth approximate ∥uk − ũ∞∥2L2((0,1);(W 1,2(Ω))∗) = ∫ tk+1 tk ∥u(·, t)− u∞∥2(W 1,2(Ω))∗dt = ∫ tk+1 tk ∥u(·, t)− u(·, tk) + u(·, tk)− u∞∥2(W 1,2(Ω))∗dt ≤ 2 ∫ tk+1 tk ∥u(·, t)− u(·, tk)∥2(W 1,2(Ω))∗dt+ 2 ∫ tk+1 tk ∥u(·, tk)− u∞∥2(W 1,2(Ω))∗dt (5.12) for all k ∈ N, whereby we deduce through (5.9) and the embedding theory C0(Ω̄) ↪→ (W 1,2(Ω))∗ that the rightmost integral provides 2 ∫ tk+1 tk ∥u(·, tk)− u∞∥2(W 1,2(Ω))∗dt = 2∥u(·, tk)− u∞∥2(W 1,2(Ω))∗ → 0 for all k → ∞, (5.13) 16 J. ZHENG, Y. WANG EJDE-2025/98 and for the first one, in consideration of the Hölder inequality and Lemma 5.3 we certainly achieve 2 ∫ tk+1 tk ∥u(·, t)− u(·, tk)∥2(W 1,2(Ω))∗dt = 2 ∫ tk+1 tk ∥∥∫ t tk ut(·, s)ds ∥∥2 (W 1,2(Ω))∗ dt ≤ 2 ∫ tk+1 tk (∫ t tk ∥ut(·, s)∥2(W 1,2(Ω))∗ds ) dt ≤ 2 ∫ tk+1 tk (∫ t tk ∥ut(·, s)∥2(W 1,2(Ω))∗ds ) (t− tk)dt ≤ 2 ∫ ∞ tk ∥ut(·, s)∥2(W 1,2(Ω))∗ds → 0 as k → ∞. (5.14) Collecting (5.13)-(5.14) and substituting into (5.12) consequently indicates that ∥uk − ũ∞∥2L2((0,1);(W 1,2(Ω))∗) → 0 as k → ∞. Whence returning to (5.10)-(5.11), it follows that ũ∞ = ū0 in Ω× (0, 1). Therefore promoting that u∞ ≡ ū0 on Ω, and (5.8) holds. □ Drawing upon Lemma 5.4 and according to standard parabolic comparison arguments, we derive the rate of exponential decay for the second component v of solutions to the system (2.1). Lemma 5.5. Assume Ω is a bounded domain with smooth boundary in RN (N ≥ 2). Then for all choices of ε ∈ (0, ū0) and some fixed s1 > 0, one has ∥v(·, t)∥L∞(Ω) ≤ ∥v0∥L∞(Ω)e −(ū0−ε)(t−s1) for all t > s1. Proof. Following Lemma 5.4, it is not difficult to discover that u(·, t) → ū0 as t → ∞ uniformly in Ω̄, where one introduces ε ∈ (0, ū0) satisfying u(·, t) ≥ ū0 − ε for all t > s1. (5.15) Once more revisting the second equation of (2.1), we obtain from the positivity of the solution component v along with (5.15) that vt ≤ ∆v − (ū0 − ε)v for all t > s1. Let y(t) denote the solution of the problem y′(t) + (ū0 − ε)y(t) = 0 for all t > s1 (5.16) with the initial condition y(s1) = ∥v(·, s1)∥L∞(Ω). (5.17) Note that [33, Lemma 2.1] allows us to deduce, for all t > 0, that t 7→ ∥v(·, t)∥L∞(Ω) is non- increasing. Consequently, an application of the comparison principle and (5.16) that upon being integrated from s1 to t yields v(·, t) ≤ y(t) ≤ y(s1)e −(ū0−ε)(t−s1) for all t > s1, whence returning to (5.17) combined with the parabolic maximum principle, we arrive at v(·, t) ≤ ∥v0∥L∞(Ω)e −(ū0−ε)(t−s1) for all t > s1, resulting in the validity of Lemma 5.5. □ On the basis of Lemma 5.5, we turn our attention to a exponential decay dynamics for the chemical concentration v in space domain W 1,p(Ω), under the requirements p > max{N, k} and k > N 2 (N ≥ 2). This plays an important role in confirming the exponentially asymptotic conver- gence of u. EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 17 Lemma 5.6. Let α0 < min{ū0 − ε, λ3} with λ3 being the first nonzero eigenvalue of −∆ in Ω under Neumann boundary conditions. For each ε ∈ (0, ū0), in the structural assumptions of p > max{N, k 2} with k as introduced in Lemma 4.2, there exists a positive constant C such that for some s2 > 0, ∥∇v(·, t)∥Lp(Ω) ≤ Ce−α0(t−t0) for all t > s2. Proof. By an associated variation-of-constants formula to the second equation in (2.1), we repre- sent v as v(·, t) = e(t−s2)∆v(·, s2)− ∫ t s2 e(t−s)∆u(·, s)v(·, s)ds for all t > s2, which in turn indicates that ∥∇v(·, t)∥Lp(Ω) ≤ ∥∇e(t−s2)∆v(·, s2)∥Lp(Ω) + ∫ t s2 ∥∇e(t−s)∆u(·, s)v(·, s)∥Lp(Ω)ds (5.18) for all t > s2. Furthermore, when in conjunction with the known results for the Neumann heat semigroup (3.4), Lemma 4.6 becomes applicable to provide constants λ3, C1, C2 > 0 such that the first term on the right-hand side of (5.18) can be estimated as ∥∇e(t−t0)∆v(·, s2)∥Lp(Ω) ≤ C1e −λ3(t−s2)∥∇v(·, s2)∥Lp(Ω) ≤ C2 for all t > s2. (5.19) Continuing, utilization of the Hölder inequality and the smoothing Lp-Lq estimates of (et∆)t≥0 in (3.2) with q := p, we identify positive constants C3 and C4 satisfying∫ t s2 ∥∇e(t−s)∆u(·, s)v(·, s)∥Lp(Ω)ds ≤ C3 ∫ t s2 ( 1 + (t− s)−1/2 ) e−λ3(t−s)∥u(·, s)v(·, s)∥Lp(Ω)ds ≤ C4|Ω|1/p ∫ t s2 ( 1 + (t− s)−1/2 ) e−λ3(t−s)∥u(·, s)∥L∞(Ω)∥v(·, s)∥L∞(Ω)ds for all t > s2. Taking advantage of Lemma 4.7 and Lemma 5.5, the above then generates the existence of a positive constant C5 such that∫ t s2 ∥∇e(t−s)∆u(·, s)v(·, s)∥Lp(Ω)ds ≤ C5∥v0∥L∞(Ω) ∫ t s2 ( 1 + (t− s)−1/2 ) e−λ3(t−s)e−α0(s−s2)ds = C5∥v0∥L∞(Ω)e −α0(t−s2) ∫ t−s2 0 ( 1 + σ−1/2 ) e−(λ3−α0)σdσ for all t > s2. (5.20) We then determine an adequately large number M0 > 0 complying with M0 ≥ 2 { C2 + C5∥v0∥L∞(Ω) ∫ t−s2 0 ( 1 + σ−1/2 ) e−(λ3−α0)σdσ } , (5.21) and let T̃0 := sup { T0 ≥ s2| ∥∇v(·, t)∥Lp(Ω) ≤ M0e −α0(t−s2) for all t ∈ [s2, T0] } , (5.22) which is well-defined and positive. A substitution of (5.19) and (5.20) into (5.18), and using (5.21) leads to ∥∇v(·, t)∥Lp(Ω) ≤ M0 2 e−α0(t−s2) for all t ∈ [s2, T̃0), actually inferring that T̃0 cannot be finite, namely, T̃0 = ∞, where in addition, we have invoked the fact that ∫ t−t0 0 ( 1 + σ−1/2 ) e−(λ3−α0)σdσ < +∞ on account of −1/2 > −1, and meanwhile, implying (5.6). □ 18 J. ZHENG, Y. WANG EJDE-2025/98 Taking into account Lemmas 5.5-5.6, we arrive at the exponential convergence rate for the first component u of solutions to the system (2.1). Lemma 5.7. For each ε ∈ (0, ū0), let α0 < min{ū0 − ε, λ3}. Then a constant C > 0 can be selected such that the solution of system (2.1) satisfies ∥u(·, t)− ū0∥L∞(Ω) ≤ Ce−α0(t−s3) for all t > s3 with some fixed s3 > 0. Proof. The variation-of-constants representation applied to the u-equation in (2.1) yields u(·, t) = e(t−s3)∆u(·, t1)− χ ∫ t s3 e(t−s)∆∇ · (u(·, s)∇v(·, s))ds for all t > s3. from elementary computations, it follows that ∥u(·, t)− ū0∥L∞(Ω) ≤ ∥e(t−s3)∆(u(·, s3)− ū0)∥L∞(Ω) + χ ∫ t s3 ∥e(t−s)∆∇ · (u(·, s)∇v(·, s))∥L∞(Ω)ds (5.23) for all t > s3. Given (3.1), there exists C1 > 0 in such a way that the first identity transforms into ∥e(t−s3)∆(u(·, s3)− ū0)∥L∞(Ω) ≤ C1e −λ3(t−s3)∥(u(·, s3)− ū0)∥L∞(Ω) for all t > s3, (5.24) while (3.3) allows for the existence of positive constants C2, C3 abiding by χ ∫ t s3 ∥e(t−s)∆∇ · (u(·, s)∇v(·, s))∥L∞(Ω)ds ≤ C2 ∫ t s3 ( 1 + (t− s)− 1 2− N 2 · 1p ) e−λ3(t−s)∥u(·, s)∇v(·, s)∥Lp(Ω)ds ≤ C3 ∫ t s3 ( 1 + (t− s)− 1 2− N 2 · 1p ) e−λ3(t−s)∥u(·, s)∥L∞(Ω)∥∇v(·, s)∥Lp(Ω)ds for all t > s3, (5.25) whence as a consequence of Lemma 5.6 and again, Lemma 4.7, one figure out C4 > 0 such that (5.25) can be further deduced that χ ∫ t s3 ∥e(t−s)∆∇ · (u(·, s)∇v(·, s))∥L∞(Ω)ds ≤ C4 ∫ t s3 ( 1 + (t− s)− 1 2− N 2 · 1p ) e−λ3(t−s)e−α0(s−s3)ds = C4e −α0(t−s3) ∫ t−s3 0 ( 1 + σ− 1 2− N 2 · 1p ) e−(λ3−α0)σdσ for all t > s3. (5.26) In light of our restriction p > max{N, k} for k > N 2 with N ≥ 2, this correspondingly means that − 1 2 − N 2 · 1 p > −1, subsequently guaranteeing that∫ t−s3 0 ( 1 + σ− 1 2− N 2 · 1p ) e−(λ3−α0)σdσ < +∞. (5.27) Thereafter, applying exactly the same arguments as in (5.21)-(5.22) enables us to collect M1 > 0 sufficiently large such that M1 ≥ 2 { C1e −λ3(t−s3)∥(u(·, s3)− ū0)∥L∞(Ω) + C4e −α0(t−s3) ∫ t−s3 0 ( 1 + σ− 1 2− N 2 · 1p ) e−(λ3−α0)σdσ } . (5.28) Furthermore, we define T̃1 := sup { T1 ≥ s3| ∥u(·, t)− ū0∥L∞(Ω) ≤ M1e −α0(t−s3) for all t ∈ [s3, T1] } . (5.29) EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 19 The above properties (5.27)-(5.29), Lemma 5.4, all the estimates provided in (5.24), and (5.26) inserted into (5.23) give rise to ∥u(·, t)− ū0∥L∞(Ω) ≤ M1 2 e−α0(t−s3) for all t ∈ [s3, T̃1). This is instrumental in T̃1 = ∞ and the verification of Lemma 5.7. □ Following rearrangement of intermediate processes above, we are now in a position to prove Theorem 2.5. 5.1. Proof of Theorem 2.5. The assertion can be readily acquired by collecting the above explanations from Lemma 5.5 and Lemma 5.7. Acknowledgements. Jiashan Zheng was supported by Shandong Provincial Natural Science Foundation (Nos. ZR2022JQ06, ZR2025MS14), and by the National Natural Science Foundation of China (Nos. 11601215, 12561035). Yuying Wang was supported by the Graduate Scientific Research and Innovation Foundation of Yantai University (No. GGIFYTU2514). The authors are grateful to the anonymous reviewers for their constructive comments, which significantly improve the quality of this article. Authors’ contributions: J. Zheng provided guidance in the theoretical aspects, and Y. Wang took charge of drafting this article. Both authors participated in the final review of the manuscript and contributed equally to the completion of this work. References [1] K. Baghaei, A. Khelghati; Boundedness of classical solutions for a chemotaxis model with consumption of chemoattractant, C. R. Math., 355 (2017), 633–639. [2] N. Bellomo, A. Bellouquid, Y. Tao, M. Winkler; Toward a mathematical theory of Keller-Segel models of pattern formation in biological tissues, Math. Models Methods Appl. Sci., 25 (2015), 1663–1763. [3] X. Cao, M. Fuest; Finite-time blow-up in fully parabolic quasilinear Keller-Segel systems with supercritical exponents, Calc. Var. Part. Differ. Equ., 64 (2025), 89. [4] F. Dai, B. Liu; Global solvability and asymptotic stabilization in a three-dimensional Keller-Segel-Navier- Stokes system with indirect signal production, Math. Models Methods Appl. Sci., 31 (2021), 2091–2163. [5] Y. Dong, S. Zhang; Global weak solvability in a self-consistent chemotaxis-Navier-Stokes system involving Dirichlet boundary conditions for the signal, Commun. Contemp. Math., (2024), 2450022. [6] K. Gao, S. Jhang, S. Shi, J. Zheng; Blow-up prevention by subquadratic and quadratic degradations in a three- dimensional Keller-Segel-Stokes system with indirect signal production, J. Differ. Equ., 423 (2025), 324–376. [7] Z. Hassan, W. Shen, Y. Zhang; Global existence of classical solutions of chemotaxis systems with logistic source and consumption or linear signal production on RN , J. Differ. Equ., 413 (2024), 497–556. [8] F. Heihoff; Can a chemotaxis-consumption system recover from a measure-type aggregation state in arbitrary dimension?, Proc. Amer. Math. Soc., 152 (2024), 5229–5247. [9] D. Horstmann, M. Winkler; Boundedness vs. blow-up in a chemotaxis system, J. Differ. Equ., 215 (2005), 52–107. [10] J. Jiang, H. Wu, S. Zheng; Global existence and asymptotic behavior of solutions to a chemotaxis-fluid system on general bounded domains, Asymptot Anal., 92 (2015), 249–258. [11] C. Jin; Global boundedness and eventual regularity of chemotaxis-fluid model driven by porous medium diffu- sion, Commun. Math. Sci., 22 (2024), 1167–1193. [12] H. Jin, Z. Wang; Global stability of prey-taxis systems, J. Differ. Equ., 262 (2017), 1257–1290. [13] E. Keller, L. Segel; Initiation of slime mold aggregation viewed as an instability, J. Theoret. Biol., 26 (1970), 399–415. [14] E. Keller, L. Segel; Traveling bands of chemotactic bacteria: a theoretical analysis, J. Theoret. Biol., 30 (1971), 235–248. [15] O. Ladyzhenskaia, V. Solonnikov, N. Ural’tseva; Linear and quasi-linear equations of parabolic type, Amer. Math. Soc., 23 (1968). [16] J. Lankeit, M. Winkler; Depleting the signal: Analysis of chemotaxis-consumption models-A survey, Stud. Appl. Math., 151(2023), 1197–1229. [17] J. Lankeit, M. Winkler, Chemotaxis-consumption interaction: solvability and asymptotics in general high- dimensional domains, (2025), arXiv:2502.17338. [18] L. Liu, A note on the global existence and boundedness of an N-dimensional parabolic-elliptic predator-prey system with indirect pursuit-evasion interaction, Open Math., 23 (2025), 20240122. [19] L. Liu, Boundedness and global existence in a higher-dimensional parabolic-elliptic-ODE chemotaxis-haptotaxis model with remodeling of non-diffusible attractant, J. Math. Anal. Appl., 549 (2025), 129473. 20 J. ZHENG, Y. WANG EJDE-2025/98 [20] L. Liu, Global well-posedness to a multidimensional parabolic-elliptic-elliptic attraction-repulsion chemotaxis system, Electron. J. Differ. Equ., 2025 (2025), No. 26, 1–20. [21] G. Ren, B. Liu; Global boundedness and asymptotic behavior in a quasilinear attraction-repulsion chemotaxis model with nonlinear signal production and logistic-type source, Math. Models Methods Appl. Sci., 30 (2020), 2619–2689. [22] Y. Tao; Boundedness in a chemotaxis model with oxygen consumption by bacteria, J. Math. Anal. Appl., 38 1(2011), 521–529. [23] Y. Tao, M. Winkler; Eventual smoothness and stabilization of large-data solutions in a three-dimensional chemotaxis system with consumption of chemoattractant, J. Differ. Equ., 252 (2012), 2520–2543. [24] Y. Tao, M. Winkler; Boundedness in a quasilinear parabolic-parabolic Keller-Segel system with subcritical sensitivity, J. Differ. Equ., 252 (2012), 692–715. [25] Y. Tao, M. Winkler; Global solutions to a Keller-Segel-consumption system involving singularly signal- dependent motilities in domains of arbitrary dimension, J. Differ. Equ., 343(2023), 390–418. [26] Y. Tao, M. Winkler; Stabilization in a chemotaxis system modelling T-cell dynamics with simultaneous pro- duction and consumption of signals, European J. Appl. Math., (2024), 1–14. [27] Y. Tian, Z. Xiang; Global boundedness to a 3D chemotaxis-Stokes system with porous medium cell diffusion and general sensitivity, Adv. Nonlinear Anal., 12(2022), 23–53. [28] H. Wang, Y. Li; Renormalized solutionsto a chemotaxissystem with consumption of chemoattractant, Electron. J. Differ. Equ., 2019 (2019), no. 38. 1–19. [29] Y. Wang, L. Pu, J. Zheng; Global existence and boundedness of classical solutions in chemotaxis-(Navier- )Stokes system with singular sensitivity and self-consistent term, Appl. Math. Lett., 166 (2025), 109518. [30] M. Winkler; Aggregation vs. global diffusive behavior in the higher-dimensional Keller-Segel model, J. Differ. Equ., 248 (2010), 2889–2905. [31] M. Winkler, Absence of collapse in a parabolic chemotaxis system with signal-dependent sensitivity, Math. Nachr., 283 (2010), 1664–1673. [32] M. Winkler; Global large-data solutions in a chemotaxis-(Navier-)Stokes system modeling cellular swimming in fluid drops, Comm. Part. Differ. Equ., 37 (2012), 319–351. [33] M. Winkler; Stabilization in a two-dimensional chemotaxis-Navier-Stokes system, Arch. Ration. Mech. Anal., 211(2014), 455–487. [34] M. Winkler; Boundedness and large time behavior in a three-dimensional chemotaxis-Stokes system with nonlinear diffusion and general sensitivity, Calc. Var. Part. Differ. Equ., 54 (2015), 3789–3828. [35] M. Winkler; Global weak solutions in a three-dimensional chemotaxis-Navier-Stokes system, Ann. Inst. H. Poincaré Anal. Non Linéaire, 33 (2016), 1329–1352. [36] M. Winkler; Global existence and stabilization in a degenerate chemotaxis-Stokes system with mildly strong diffusion enhancement, J. Differ. Equ., 264 (2018), 6109–6151. [37] M. Winkler; Chemotaxis-Stokes interaction with very weak diffusion enhancement: blow-up exclusion via detection of absorption-induced entropy structures involving multiplicative couplings, Adv. Nonlinear Stud., 22(2022), 88–117. [38] M. Winkler; Effects of degeneracies in taxis-driven evolution, Math. Models Methods Appl. Sci., 2025. [39] J. Xie, J. Zheng; A new result on existence of global bounded classical solution to a attraction-repulsion chemotaxis system with logistic source, J. Differ. Equ., 298 (2021), 159–181. [40] Q. Zhang, Y. Li; Stabilization and convergence rate in a chemotaxis system with consumption of chemoattrac- tant, J. Math. Phys., 56 (2015), 081506. [41] Q. Zhang, Y. Li; Convergence rates of solutions for a two-dimensional chemotaxis-Navier-Stokes system, Discrete Contin. Dyn. Syst. Ser. B, 20 (2015), 2751–2759. [42] J. Zheng; Eventual smoothness and stabilization in a three-dimensional Keller-Segel-Navier-Stokes system with rotational flux, Calc. Var. Part. Differ. Equ., 61 (2022), 52. [43] J. Zheng, Y. Ke; Further study on the global existence and boundedness of the weak solution in a three- dimensional chemotaxis-Stokes system with nonlinear diffusion and general sensitivity, Commun. Nonlinear Sci. Numer. Simul., 115 (2022), 106732. [44] J. Zheng, Y. Ke; Global existence and eventual smoothness for a 2-dimensional Keller-Segel system with nonlinear logistic source, J. Differ. Equ., 437 (2025), 113293. [45] J. Zheng, D. Qi; Global existence and boundedness in an N-dimensional chemotaxis-Navier-Stokes system with nonlinear diffusion and rotation, J. Differ. Equ., 335 (2022), 347–397. [46] J. Zheng, D. Qi, Y. Ke, Global existence, regularity and boundedness in a higher-dimensional chemotaxis- Navier-Stokes system with nonlinear diffusion and general sensitivity, Calc. Var. Part. Differ. Equ., 61 (2022), 150. [47] P. Zheng, H. Zhang, Global stability in a fully parabolic three-species chemotaxis-competition system, J. Math. Phys., 66, 2025. Jiashan Zheng School of Mathematics and Information Sciences, Yantai University, Yantai 264005, China Email address: zhengjiashan2008@163.com EJDE-2025/98 PARABOLIC-PARABOLIC SYSTEMS WITH CONSUMPTION 21 Yuying Wang (corresponding author) School of Mathematics and Information Sciences, Yantai University, Yantai 264005, China Email address: wangyuy2024@163.com 1. Introduction 2. Main results 3. Preliminaries and local well-posedness 4. A priori estimates. Proof of Theorem 2.1 4.1. Proof of Theorem 2.1 5. Exponential decay. Proof of Theorem 2.5 5.1. Proof of Theorem 2.5 Acknowledgements References