Global solutions to a haptotaxis system with a potentially degenerate diffusion tensor in two and three dimensions

We consider the potentially degenerate haptotaxis system \begin{equation*} \left\{ \begin{aligned} u_t&= \nabla \cdot (\mathbb{D} \nabla u + u \nabla \cdot \mathbb{D}) - \chi \nabla \cdot (u\mathbb{D}\nabla w) + \mu u(1-u^{r- 1}), \\ w_t&= - uw \end{aligned} \right. \end{equation*} in a smooth bounded domain $\Omega \subseteq \mathbb{R}^n$, $n \in \{2,3\}$, with a no-flux boundary condition, positive initial data $u_0$, $w_0$ and parameters $\chi>0$, $\mu>0$, $r \geq 2$ and $\mathbb{D}: \overline{\Omega} \rightarrow \mathbb{R}^{n\times n}$, $\mathbb{D}$ positive semidefinite on $\overline{\Omega}$. Our main result regarding the above system is the construction of weak solutions under fairly mild assumptions on $\mathbb{D}$ as well as the initial data, encompassing scenarios of degenerate diffusion in the first equation. As a step in this construction as well as a result of potential independent interest, we further construct classical solutions for the same system under a global positivity assumption for $\mathbb{D}$, which ensures the full regularizing influence of its associated diffusion operator. In both constructions, we naturally rely on the regularizing properties of a sufficiently strong logistic source term in the first equation.


Introduction
As the movement of cells plays a significant role in many biological systems and processes, analyzing the underlying mechanisms can prove useful in understanding these systems and processes themselves. One such process, which is naturally of extensive interest, is the invasive movement of tumor cells into healthy tissue along gradients of tissue density during the progression of certain types of cancer, which is governed by a mechanism generally called haptotaxis (cf. [8]). Similarly to the efforts made to understand the related process of chemotaxis (cf. [3]), which models movement along gradients of a diffusive chemical as opposed to non-diffusive tissue, mathematical modeling of haptotaxis has proven to be a fruitful area of study. In both cases, by far the most attention at this point has been paid to approaches employing a Fickian diffusive movement model for the organisms in question, which assumes some homogeneity of the underlying medium. But bolstered by experiments regarding cell aggregation near interfaces between grey and white matter in mouse brains (cf. [6]), it has recently been suggested that especially in more heterogeneous environments, such as brain tissue, cell movement might be better described by non-Fickian diffusion (cf. [4]), which is far less mathematically studied in these taxis settings.
In an effort to add to the base of knowledge in this area, we will focus our efforts here on a haptotaxis model of cancer invasion featuring such non-Fickian myopic diffusion, which was introduced in [11]. More specifically, we consider the system u t = ∇ · (D∇u + u∇ · D) − χ∇ · (uD∇w) + µu(1 − u r−1 ), in a smooth bounded domain Ω ⊆ R n , n ∈ {2, 3}, with a no-flux boundary condition and appropriate parameters χ > 0, µ > 0, r ≥ 2 and D : Ω → R n×n , D positive semidefinite on Ω. The first equation models the invading cancer cells moving according to the aforementioned myopic diffusion, which is represented by the term ∇ · (D∇u + u∇ · D), as well as according to haptotaxis, which is represented by the term −χ∇ · (uD∇w). Apart from this, the equation further incorporates a logistic source. The second equation models the remaining healthy tissue cells and only features a consumption term.
The key feature of interest in the above system from both an application as well as a mathematical perspective is of course the parameter matrix D, which represents a space dependent coupled diffusion and taxis tensor. In practice, this tensor can be derived from the underlying tissue structure by employing direct imaging methods (cf. [11]) and represents the influence of said underlying structure on the movement of cells through it. To account for situations of both locally very dense as well as locally very sparse tissue, which both occur in concrete applications and hinder cell movement significantly, we allow D to be potentially degenerate.
Notably in one dimension, solutions to a closely related system with degenerate diffusion have already been shown to reflect the aggregation behavior in interface regions seen in experiments (cf. [6], [39]) while to our knowledge in systems with non-degenerate diffusion long-time behavior results seem to generally be restricted to homogenization (cf. e.g. [37]). This seems to indicate that models of this kind featuring degenerate diffusion could potentially be a better representation of real world behavior. As such, the development of methods to cope with the challenges of this model, namely the reduced regularizing effect of the degenerate diffusion operator and the destabilizing effects of taxis, while still allowing for a sufficiently large class of matrices to enable the modeling of real world scenarios seems to be a worthwhile endeavor, which to our knowledge has thus far only been addressed in one dimension. Thus, the aim of this paper is the investigation of the apparently still open question whether solutions exist in two and three dimensions even if the diffusion operator is degenerate.

Results.
Our main results regarding the haptotaxis model described above are twofold: First, we establish the existence of global classical solutions given a uniform positivity condition for D, which allows us to basically treat it as we would any other elliptic diffusion operator, as well as a condition ensuring sufficient regularizing influence of the logistic source term. Second, we establish that it is still possible to construct fairly standard weak solutions under much more relaxed conditions for D. More specifically, we drop the assumption that D must be globally positive in Ω and replace it with a set of assumptions much more tailored to our methods for constructing said weak solutions, which are strictly weaker than the prior positivity assumption in allowing for matrices that are in some (small) parts of Ω only positive semidefinite.
Given that the definitions necessary to properly formulate the above results take up significant space, we will not go into more detail here but instead refer the reader to the very next section for the pertinent details regarding said results. As an addition to stating precise versions of our results in the next section, we will further discuss the prototypical examples of a matrix with a single point degeneracy as well as a matrix with degeneracies on a manifold of higher dimension and derive some conditions, under which they still allow for the construction of weak solutions. We do this to help build some additional intuition in parallel to the rather abstract regularity properties introduced in said section as well as to illustrate that our results can work for some scenarios with real world relevance such as e.g. a domain divided by an impenetrable membrane.
Approach. Let us now give a brief sketch of the methods employed to achieve our two main results.
For our classical existence result, we begin by using standard contraction mapping methods to gain local solutions with an associated blow-up criterion as the operator in the first equation is strictly elliptic for globally positive matrices D. We then immediately transition to analyzing the function a := ue −χw , which together with w solves the closely related problem (3.3). We do this because, in a sense, this transformation eliminates the problematic cross-diffusive term from the first equation by integrating it into the function a and its diffusion operator. Using a fairly classical Moser-type iteration argument, we then establish an L ∞ (Ω) bound for a, which translates back to u. Using this bound combined with two testing procedures then yields a further W 1,4 (Ω) bound for w, which together with the already established bound for u is sufficient to ensure that finite time blow-up is in fact impossible in two and three dimensions and thus completes the proof of our first result.
Regarding our second result, we begin by approximating the initial data, the matrix D as well as the logistic source term in such a way as to make the already established global classical existence result applicable to the in this way approximated versions of (1.1). For the family of solutions (u ε , w ε ), ε ∈ (0, 1), gained in this fashion, we then establish a bound of the form by way of an energy-type inequality, which already proved useful in the one-dimensional case discussed in [39]. Using this as a baseline, we then derive the bounds necessary for applications of the Aubin-Lions compact embedding lemma to gain our desired weak solutions as limits of the approximate ones.
Prior work. As haptotaxis models (cf. [36] for a general survey) as well as the closely related chemotaxis models (cf. [3] for a general survey) have been extensively studied in many possible variations since the introduction of their progenitor in the seminal 1970 paper by Keller and Segel (cf. [16]), there is of course a lot of prior art available regarding global existence theory for said models. While it is certainly out of scope for this paper to cover prior results in their entirety, we will nonetheless give an overview of some notable ones.
Let us first note that for the one-dimensional case, where D simplifies to a real-valued function, there are already some results available for a variant of our scenario without a logistic source term (including potential spacial degeneracy) dealing with existence theory as well as long time behavior (cf. [39], [40] and [41]). Weak solutions have also been constructed in very similar haptotaxis systems featuring porous-medium type and signal-dependent degeneracies as opposed to spacial ones (cf. [45]).
Regarding haptotaxis system with non-degenerate diffusion operators, e.g. D ≡ 1 in our system, global existence and sometimes boundedness theory has been studied in various closely related settings (cf. [7], [19], [21], [29], [34], [35], [42]). Notably, these systems often feature an additional equation modeling a diffusive (potentially attractive) chemical and the fixed parameter choice r = 2 for the logistic term in addition to the more regular diffusion. In many of these scenarios, it has further been established that solutions converge to their constant steady states (cf. [20], [21], [25], [29], [37], [44]) under varied but sometimes restrictive assumptions. There has also been some analysis of haptotaxis with tissue remodeling, which is represented in the model by some additional source terms in the equation for w (cf. [24], [27], [30]).
Apart from haptotaxis models, there has also been significant analysis of chemotaxis models featuring degenerate diffusion (cf. [10], [18], [43] including degeneracies depending on the cell density itself).
Lastly, let us just briefly mention that the regularizing effects of logistic source terms we rely on in this paper have already been very well-documented in various chemotaxis systems (cf. [17], [38] among many others) as well as haptotaxis systems (cf. [31]).

Main Results and Related Definitions
As already alluded to in the introduction, we will focus our attention in this paper on the system in a smooth bounded domain Ω ⊆ R n , n ∈ {2, 3}, with parameters χ > 0, µ > 0, r ≥ 2, D : Ω → R n×n , D positive semidefinite on Ω, and some initial data u 0 , w 0 : Ω → [0, ∞).
Our results concerning this system are twofold. We will first derive the following existence result concerning global classical solutions in two and three dimensions under the assumptions that D and the initial data are sufficiently regular, D is positive definite on Ω and the logistic source term is sufficiently strong.
This result, while of course also of independent interest, will then serve as a building block for the construction of weak solutions to the same system under much more relaxed restrictions on D and the initial data.
Chiefly, global positivity of the matrix D is not necessarily needed anymore and is instead replaced by a set of much weaker but more specific regularity assumptions.
The first such regularity property concerns the divergence of D (applied column-wise) and how it can be estimated by the (potentially degenerate) scalar product induced by D.

Definition 2.2.
Let Ω ⊆ R n , n ∈ N, be a bounded domain with a smooth boundary. We then say a positive semidefinite for all Φ ∈ C 0 (Ω; R n ).

Remark 2.4.
It is fairly easy to verify that any smooth, positive definite D allows for such an estimate with the optimal exponent β = 1 2 . Let us therefore now briefly illustrate that the above property is also achievable for less regular D, which are e.g. at some points in Ω only positive semidefinite, by giving some examples. While we will not necessarily fully explore these examples and leave out some of the more cumbersome corner cases for ease of presentation, they will accompany us throughout this section as a tool to give some intuition for later introduced definitions as well as to give concrete examples for degenerate cases in which weak solutions can still be constructed.
We will first take a look at the prototypical case of a matrix-valued function D 1 on a ball with a single degenerate point in the origin, or more precisely we will consider D 1 (x) := |x| s I on Ω := B 1 (0) ⊆ R n , n ∈ N, I being the identity matrix and s being some positive real number.
To illustrate that our framework also supports analysis of singularities occurring on higher dimensional manifolds, let us further consider the similar prototypical example D 2 (x 1 , . . . , x n ) := |x 1 | s I on the same set Ω with s now being a real number greater than 1. As here ∇ · D 2 (x 1 , . . . , x n ) = (s|x 1 | s−2 x 1 , 0, . . . , 0) almost everywhere, we gain that D 2 has the property laid out in Definition 2.2 for all β ∈ ( 1 s , 1) ∩ ( 1 2 , 1) by a similar argument as for the previous example.
As to be expected in both cases, smaller values of s result in the divergence estimate only holding for ever larger exponents β. As we will see in our theorem regarding the existence of weak solutions at the end of this section, these larger values of β will necessitate stronger regularizing influence from the logistic source term to compensate.
Before we can now approach the second regularity property of this section as well as properly defining what we in fact mean by weak solutions in this paper, we need to first introduce a set of function spaces. Said spaces are generally fairly straightforward generalizations of standard Sobolev and Lebesgue spaces incorporating D as well as some spaces derived from them, which are more specific to our setting. For a more thorough discussion of e.g. the degenerate Sobolev spaces introduced below, we refer the reader to [26].
We will further take the introduction of said spaces as an opportunity to present some of their most important properties for our purposes immediately after defining them.

Definition 2.5.
Let Ω ⊆ R n , n ∈ N, be a bounded domain with a smooth boundary and p ∈ [1, ∞).
We then define the Sobolev-type space Let now D ∈ C 0 (Ω; R n×n ) be positive semidefinite everywhere. We then define the Lebesgue-type space L p D (Ω) as the set of all measurable R n -valued functions Φ on Ω with finite seminorm Furthermore, we define the Sobolev-type spaces W 1,p D (Ω) as the completion of C ∞ (Ω) in the norm in the same vain as the standard Sobolev spaces. It is straightforward to see that each space W 1,p D (Ω) can be interpreted as a subspace of L p (Ω) × L p D (Ω) in a natural way and thus elements of these spaces can be written as tuples (ϕ, Φ). As such, there exist the natural continuous projections associated with this representation.
Remark 2.6. For a more comprehensive exploration of these spaces and their properties see e.g. [26].
We will now give a brief overview of the properties the above spaces retain from the standard Sobolev and Lebesgue spaces as well as some of the differences. As most of the proofs translate directly from standard Sobolev theory or are laid out in [26], we will only list the properties we are interested in without extensive argument.
First of all by construction, W 1,p div (Ω; R n×n ), L p D (Ω) and W 1,p D (Ω) are Banach spaces, which are reflexive if p ∈ (1, ∞), by essentially the same arguments as for the standard Sobolev and Lebesgue spaces and, for p = 2, they are in fact Hilbert spaces with the natural inner products. It is further easy to see that, if (ϕ, Φ) is a strong or weak limit of a sequence (ϕ n , Φ n ) n∈N ⊆ W 1,p D (Ω), the function ϕ ∈ L p (Ω) coincides with the pointwise almost everywhere limit of the sequence (ϕ n ) n∈N if it exists due to P 1 being continuous regarding both topologies and well-known results about strong and weak convergence in L p (Ω).
As opposed to the classical Sobolev spaces, the spaces W 1,p D (Ω) can not necessarily be understood as subspaces of the spaces L p (Ω) because their equivalents to the weak gradients in the classical Sobolev spaces are not necessarily unique here, meaning essentially that P 1 is not always injective. (For an example of this, see [26, p. 1877]). Given that this can be problematic when deriving analogues to the (compact) embedding properties of Sobolev spaces for our weaker variants, let us now briefly note that, under sufficient regularity assumptions for D, the spaces W 1,p D (Ω) do in fact embed into the spaces L p (Ω) again. In particular if p = 2, which is the parameter choice we are most interested here, this is the case if √ D ∈ W 1,2 div (Ω; R n×n ) according to Lemma 8 from [26].
While it presents a slight abuse of notation, we will in a similar fashion to [26] use ϕ to mean P 1 (ϕ) ∈ L p (Ω) for elements ϕ ∈ W 1,p D (Ω) when unambiguous and generally use the convention ∇ϕ = P 2 (ϕ) even if ∇ϕ is not necessarily the actual weak derivative. We will further often write to simplify the notation in later arguments. If ϕ is additionally an element of C 1 (Ω), we will always assume ∇ϕ to be equal to the classical derivative, of course.
Having established these function spaces, we can now clearly state the second and last regularity property for D we are interested in. It is a simple compact embedding property, which is mainly used in this paper to facilitate application of the well-known Aubin-Lions lemma.

Definition 2.7.
Let Ω ⊆ R n , n ∈ N, be a bounded domain with a smooth boundary. We say a positive semidefinite D ∈ C 0 (Ω; R n×n ) allows for a compact L 1 (Ω) embedding if W 1,2 D (Ω) embeds compactly into L 1 (Ω).

Remark 2.8.
Let us briefly note that any D, which is equal to zero on any open subset U of Ω, cannot fulfill the property laid out in Definition 2.7 as it is well documented that L 2 (U ), which is equal to W 1,2 D (U ) in this case, does not embed compactly into L 1 (U ).
We will now give some additional criteria for the above compact embedding property to not only make our results easier to use in application but also to help us prove that both of the examples discussed in Remark 2.4 in fact fulfill it.

Lemma 2.9.
Let Ω ⊆ R n , n ∈ N, be a bounded domain with a smooth boundary and N ⊆ Ω be a relatively closed set in Ω with measure zero. Let then then D allows for a compact L 1 (Ω) embedding.
Proof. Due to our assumption that √ D ∈ W 1,2 div (Ω; R n×n ), Lemma 8 from [26] immediately yields that the projection P 1 : W 1,2 D (Ω) → L 2 (Ω) ⊆ L 1 (Ω) from Definition 2.5 is injective and thus provides us with a continuous embedding of W 1,2 D (Ω) into L 1 (Ω). It thus only remains to show that this embedding is in fact compact given the various criteria outlined above.
To do this, we first fix a bounded sequence (ϕ k ) k∈N ⊆ W 1,2 D (Ω). We then only need to construct a subsequence of (ϕ k ) k∈N that converges in L 1 (Ω) to some function ϕ to prove our desired outcome. As it is further possible to find another sequence ( , we can further assume that ϕ k ∈ C ∞ (Ω) for all k without loss of generality. If W 1,2 D (Ω) now embeds compactly into L 1 loc (Ω \ N ), we can choose a subsequence (ϕ kj ) j∈N and function ϕ : Ω → R such that ϕ kj → ϕ in all L 1 (Ω N,ε ), ε > 0, as j → ∞. As by our assumptions N ∪ ∂Ω is closed and thus Ω \ N = k∈N Ω N,1/k , we can then employ a standard diagonal sequence argument to gain yet another subsequence, which we will again call (ϕ kj ) j∈N for convenience, with the property that ϕ kj → ϕ almost everywhere in Ω \ N and thus almost everywhere in all of Ω as j → ∞ because N is a null set. Given that the thus constructed subsequence is further bounded in L 2 (Ω) due to it being bounded in W 1,2 D (Ω), we can use Vitali's theorem and the de La Valleé Poussin criterion for uniform integrability (cf. [9, pp. 23-24]) to conclude that ϕ kj → ϕ in L 1 (Ω) as well, yielding the first part of our result.
If D is positive definite on Ω \ N , then for every ε > 0 there exists K(ε) > 0 such that D > K(ε) on Ω N,ε due to the continuity of D and the fact that Ω N,ε ⊆ Ω \ N is compact. Thus, the norms of the spaces W 1,2 (Ω N,ε ) and W 1,2 D (Ω N,ε ) are equivalent. As such, the sequence (ϕ k ) k∈N is bounded in all of the spaces W 1,2 (Ω N,ε ), ε > 0. Due to our further assumption that there exists ε 0 > 0 such that W 1,2 (Ω N,ε ) embeds compactly into L 1 (Ω N,ε ) for all ε ∈ (0, ε 0 ) and the fact that any compact set K ⊆ Ω \ N is a subset of some Ω N,ε as another consequence of Ω \ N being equal to ε>0 Ω N,ε , a standard diagonal sequence argument yields a subsequence along which the functions ϕ k converge to some ϕ in L 1 loc (Ω \ N ). Combining this with the arguments from the previous paragraph then yields the second part of our result.
To now complete the proof, we first note that [1, Theorem 6.3] states that a Lipschitz boundary condition for the sets Ω N,ε ensures the Sobolev embedding necessary for our second result and thus the third result follows directly from the second. √ D 1 , √ D 2 ∈ W 1,2 div (Ω; R n×n ) in dimensions two or higher. Furthermore due to the fairly straightforward geometry of the degeneracy set N in both cases, it is easy to verify that both examples also fulfill the third criterion in Lemma 2.9 and thus both D 1 and D 2 allow for a compact L 1 (Ω) embedding in accordance with Definition 2.7.
While we have now invested some effort into formalizing the restrictions on D necessary for our later construction of weak solutions, we have yet to clarify what we in fact mean by a weak solution to (2.1). Let us now rectify this in the following definition.

Definition 2.11.
Let Ω ⊆ R n , n ∈ N, be a bounded domain with a smooth boundary and let χ > 0, µ > 0 We then call a tuple of functions ). As we have at this point clearly defined the target and some of the preconditions, let us now outright state the second main theorem we endeavor to prove in this paper.
be positive semidefinite everywhere. Let further D allow for a divergence estimate with exponent β (cf. Definition 2.2) and let D allow for a compact L 1 (Ω) embedding (cf. Definition 2.7). Finally, let u 0 ∈ L z[ln(z)]+ (Ω) and w 0 ∈ C 0 (Ω) be some initial data with Then there exist a.e. non-negative functions that are a weak solution to (2.1) in the sense of Definition 2.11.  the same holds true for D 2 . Further due to the arguments presented in Remark 2.10, both D 1 and D 2 allow for a compact L 1 (Ω) embedding. Therefore, the above theorem means that, for sufficiently regular initial data u 0 , w 0 and if either D = D 1 and r and s satisfy (2.5) or D = D 2 and r and s satisfy (2.6), weak solutions to (2.1) in fact exist in two and three dimensions.

Existence of Classical Solutions
As the existence of classical solutions to (2.1), apart from being an interesting result by its own merits, plays an important role in our construction of their weak counterparts, we will in this section first focus on their derivation. As such, our ultimate goal for this section will be the proof of our first main result, namely Theorem 2.1. The methods presented here will in many ways mirror those for similar systems with a standard Laplacian as diffusion operator. We mainly verify that the differing elements in our systems do not impede said methods.
Comparing the very strong regularity assumptions for D in this section to the much weaker ones in the following section devoted to the construction of weak solutions, the question why the gap in assumed regularity between these sections is as large as it is naturally presents itself. Let us therefore briefly address this issue. It is certainly possible to derive most of the a priori estimates, which are used in this section to argue that blow-up of local solutions is impossible, under similarly specific regularity assumptions as seen in Definition 2.2 or Definition 2.7 (albeit with some additions). But generalizing the theory employed by us to first gain said local solutions with less regular D would necessitate Schauder and semigroup theory for potentially very degenerate operators, which is out of scope for this paper. Furthermore, we think that this result is already of interest in and of itself.

Existence of Local Solutions
After this introductory paragraph giving our rational for the assumptions about D in this section, we will now focus on the construction of local solutions to the system (2.1) as a first step in constructing global ones. As for a positive definite matrix D, the diffusion operator in the first equation is strictly elliptic and therefore accessible to most of the same existence and regularity theory as the Laplacian, we will not go into detail concerning the construction of local solutions but rather refer the reader to a local existence result for a similar haptotaxis system with our operator replaced by the Laplacian in [31].
is a classical solution to (2.1) on Ω × (0, T max ) with initial data (u 0 , w 0 ) and satisfies the following blow-up criterion: For ease of further discussion, we now fix such a maximal local solution (u, w) on (0, T max ) with initial data (u 0 , w 0 ) and the parameters as stated in the above introductory paragraphs.
Before diving into the derivation of more substantial bounds for the above solution, we derive a straightforward mass bound for the first solution component as well as an L ∞ (Ω) bound for the second solution component. These bounds will not only prove useful when ruling out blow-up in this section but also serve as a baseline for bounds derived in our later efforts focused on the construction of weak solutions.
Proof. Integrating the first equation in (2.1) and applying partial integration yields for all t ∈ (0, T ) and therefore immediately give us the first half of our result by time integration. Given that further w t ≤ 0 due to the second equation in (2.1), the second half of our result follows directly as well.

A Priori Estimates
The next natural step after establishing local solutions with an associated blow-up criterion is of course arguing that finite-time blow-up is impossible and the maximal local solutions were in fact global all along.
To do this, we will devote this section to a set of a priori estimates, which increase in strength as the section goes on until they rule out blow-up of both u and w.
As is not uncommon in the analysis of these kinds of haptotaxis systems (cf. [31]), we will from now consider the function a := ue −χw defined on Ω × [0, T max ) and its associated initial data a 0 := u 0 e −χw0 defined on Ω in addition to the actual solutions components u and w themselves. A simple computation then shows that (a, w) is a classical solution of the following related system:  The key property of the above system, which makes it so useful for our purposes, is that it in a sense eliminates the taxis term or at least the explicit gradient of w from the first equation (by in a sense integrating it into a and its diffusion operator). This alleviates many of the normal problems associated with the taxis term in testing or semigroup based approaches used to derive a priori estimates. A second useful property of this transformation is that, by definition, bounds that do not involve derivatives are easily translated back from a to u as we will see later. Note however that, as soon as we want to back propagate bounds about the gradient of a to u, the complications introduced by the taxis term come back into play, making this transformation much less useful for endeavors of this kind.
We now begin by translating the baseline estimates given in Lemma 3.2 to our newly defined function a as we will henceforth focus on (a, w) as our central object of analysis for quite some time. We will further for the foreseeable future work under the assumption that T max < ∞ as this is exactly the case we want to rule out by leading this assumption to a contradiction with the blow-up criterion.

Proof.
As Ω a = Ω ue χw ≤ e χ w L ∞ (Ω) Ω u, this is a direct consequence of Lemma 3.2 if T max < ∞. In preparation for a later Moser-type iteration argument for the first solution component a (cf. [2] and [23] for some early as well as [14] and [28] for some more contemporary examples of this technique), which will later be used to rule out its finite-time blow-up, we will now derive a recursive inequality for terms of the form Ω a p . This recursion will in fact allow us to estimate each term of the form Ω a p by terms of the form ( Ω a p 2 ) 2 with constants independent of p, which will prove sufficient to later gain an L ∞ (Ω) bound for a. The method employed to gain said recursion is testing the first equation in (3.3) with e χw a p−1 followed by some estimates based on the Gagliardo-Nirenberg inequality.
To facilitate this derivation of said recursion, we will from now on assume that the regularizing influence of the logistic source term in the first equation of (2.1) is sufficiently strong, or more precisely we assume that either r > 2 or µ is sufficiently large in comparison to χ and the L ∞ (Ω) norm of w 0 . However at this point and therefore for the whole of the Moser-type iteration argument, we will not use our assumed restriction to two or three dimensions just yet.
Proof. We test the first equation in (3.3) with e χw a p−1 and apply partial integration to see that for all t ∈ (0, T max ) and p ≥ 2. Given our assumptions for D in (3.1), we can use Young's inequality to further estimate that as well as more elementary that for all t ∈ (0, T max ) and p ≥ 2, which when applied to (3.4) results in for all t ∈ (0, T max ) and p ≥ 2. If r > 2, we can now further estimate that for all t ∈ (0, T max ) and p ≥ 2 by Young's inequality. If, however, r = 2 and µ ≥ χ w 0 L ∞ (Ω) , it is immediately obvious that with K 1 := 1 for all t ∈ (0, T max ) and p ≥ 2. As such, we can in both cases conclude from (3.5) that with K 2 := (µ + 2M 3 )e χ w0 L ∞ (Ω) for all t ∈ (0, T max ) and p ≥ 2.
We can now use the Gagliardo-Nirenberg inequality to fix a constant K 3 > 0 such that for all t ∈ (0, T max ) and p ≥ 2 with . Applying this to (3.6) then implies for all t ∈ (0, T max ) and p ≥ 2. Time integration then yields for all t ∈ (0, T max ) and p ≥ 2 as T max < ∞, which after estimating the sum on the right-hand side by thrice the maximum of its summands completes the proof.
We will now proceed to give the actual iteration argument yielding an L ∞ (Ω)-type bound for a and therefore u, which is sufficient to rule out finite-time blow-up for the first solution component u.
Proof. Let p i := 2 i , i ∈ N 0 , and J i := sup t∈(0,Tmax) Ω a pi (·, t) 1 p i . Then J 0 is finite because of Corollary 3.3 and the fact that p 0 = 1. We further know that Due to Lemma 3.4, we can conclude that there exists a constant K 2 ≥ 1 such that the numbers J i conform to the following recursion: Iterating this recursion finitely many times ensures that all J i are finite.
If there exists an incrementing sequence of indices i ∈ N, along of which J i ≤ max(K 1 K 2 , K 3 2 ), we immediately gain our desired result by taking the limit of J i along said sequence. As such, we can now assume that there exists i 0 ∈ N with to cover the remaining case. Given these assumptions, the above recursion simplifies to for all i ≥ i 0 with some K 3 > 0 (only depending on K 2 ) as the function z → (zK 2 ) K 2 √ z is bounded on [1, ∞). By now again iterating this recursion finitely many times, we gain that for all i ≥ i 0 due to the series on the right side being of geometric type, we can conclude from (3.7) that the sequence J i is uniformly bounded. Therefore, taking the limit i → ∞ gives us our desired bound for a. As u = ae χw , the corresponding bound for u follows directly from this and Lemma 3.2.
To now establish that finite-time blow-up of the second solution component w is equally as impossible, we will begin by testing the first equation in (3.3) with −∇·(D∇a) and combining the result with the differential equation associated with d dt Ω |∇w| 4 . The key to extracting a sufficiently strong bound for w is to then use the strength of the absorptive terms originating from the fully elliptic operator −∇ · (D∇·) to counteract the influence of potentially destabilizing terms due to the haptotaxis interaction. Note that the ellipticity of the operator is ensured because we assume that D is positive definite everywhere in Ω. Lemma 3.6. If T max < ∞ and further r > 2 or µ ≥ χ w 0 L ∞ (Ω) , then there exists a constant C > 0 such that ∇w(·, t) L 4 (Ω) ≤ C for all t ∈ (0, T max ).
Proof. Given Lemma 3.5, we can fix a constant K 1 ≥ 1 such that a(·, t) L ∞ (Ω) ≤ K 1 and Ω a 2 (·, t) + a 2r (·, t) + a 4 (·, t) ≤ K 1 for all t ∈ (0, T max ). (3.8) Using the Gagliardo-Nirenberg inequality and standard regularity estimates (cf. [ This in turn implies that for all t ∈ (0, T max ) with K 3 := K 3 1 K 2 . After establishing these preliminaries, we now note that the first equation in (3.3) can also be written as We then test this variant of said equation with −∇ · (D∇a) and employ partial integration (using the fact that (∇ · D) · ν = 0 on ∂Ω) as well as Young's inequality to conclude that for all t ∈ (0, T max ) with K 4 := 8 max µ, µe χ(r−1) w0 L ∞ (Ω) , χ w 0 L ∞ (Ω) e χ w0 L ∞ (Ω) 2 . Using the bounds outlined in (3.1) and (3.8), we can now further derive that for all t ∈ (0, T max ). Applying these three estimates combined with the second bound in (3.8) to (3.10) then yields 1 2 As our second step, we now obtain the following estimate for the time derivative of certain gradient terms of the second solution component w as follows: for all t ∈ (0, T max ) with K 7 := w 0 L ∞ (Ω) e χ w0 L ∞ (Ω) . Now combining this with (3.11) (using an appropriate scaling factor) we gain for all t ∈ (0, T max ) with K 8 := K 5 + 1 4K3 . The application of (3.9) to the inequality above then yields with K 9 := 16K 3 K 7 K 8 for all t ∈ (0, T max ), which, by a standard comparison argument and the assumption that T max is finite, directly gives us our desired result.
Remark 3.7. The result of the above lemma only ensures that finite-time blow-up of the second solution component is impossible in two and three dimensions according to our blow-up criterion (3.2). As such, it is at this point and only this point in this section, where our restriction to two or three dimensions becomes necessary. This, of course, in turn means that any extension of the results of this section to a higher dimensional setting would only need to extend the above argument to one providing better bounds for the gradient of w.
Given that Lemma 3.5 and Lemma 3.6 rule out any kind of finite-time blow-up for our local solutions, the proof of the first central result of this paper can now be stated quite succinctly.
Proof of Theorem 2.1. If we assume T max < ∞, Lemma 3.5 and Lemma 3.6 in combination contradict the consequence of the blow-up criterion (3.2) in this case. Therefore, T max = ∞ and thus the local solutions constructed in Lemma 3.1 must be in fact global. This is sufficient to prove Theorem 2.1 as the fixed assumptions of this section were in fact identical to those of said theorem.
Remark 3.8. It is also possible to construct classical solutions in the two dimensional case without relying on logistic influences by using some methods that have previously been used when for example dealing with standard diffusion and some slightly modified versions of our arguments (cf. [3]).
Essentially, the argument boils down to using an estimate of the form with ε being potentially arbitrarily small (cf. [5, p.1199]) in combination with an additional baseline Ω u ln(u) estimate based on an energy-type inequality (cf. Lemma 4.2) to establish an L 2 (Ω) estimate. From there, the arguments are very similar to the Moser-type iteration argument presented above, only with some slight complications added, which are easily surmountable. Lemma 3.6 translates basically verbatim.
We decided not to present this result here as it will not be needed for our later construction of weak solutions and is not appreciably different from what we have done here or has already been done in the classical diffusion case.

Existence of Weak Solutions
We have at this point established all the classical existence theory we want to address in this paper and therefore will now transition to our construction of weak solutions, which is in part based on said classical theory.

Approximate Solutions
As is fairly common, our construction of weak solutions will centrally rely on approximation of said solutions by classical solutions, which solve a suitably regularized version of the original problem. As we already derived global existence of classical solutions for the system (2.1) with very strong assumptions on D, we of course want to construct our weak solutions under much weaker assumptions on D because there would be almost nothing gained otherwise. As such, the central regularization employed by us will be concerned with approximating a potentially quite irregular D by matrices D ε that are sufficiently regular to ensure classical existence of solutions. Apart from this, we will use approximated initial data. We will also slightly modify the logistic source term to ensure r > 2 in our approximated system because we can then further eliminate the assumption concerning the parameters χ and µ needed for the classical theory when r = 2. One central advantage of this approach is that our approximate systems are very close to the system we actually want to construct solutions for and thus our regularizations only minimally interfere with the structures present in the system, which we want to exploit for e.g. a priori information.
As for any β ∈ [ 1 2 , 2 3 ] the condition β 1−β ≤ r is always fulfilled independent of our choice of r ∈ [2, ∞) and as it is easy to see that, if D allows for a divergence estimate in accordance with Definition 2.2, it also allows for a divergence estimate with any larger exponent, we can assume that the parameter β seen in the second of the above properties is in fact an element of [ 2 3 , 1) ⊆ ( 1 2 , 1) without loss of generality. Then according to Remark 2.3, the aforementioned divergence estimate directly implies that with q := 2β 2β−1 . Given these assumptions, we now choose an approximate family (D ε ) ε∈(0,1) ⊆ C 2 (Ω; R n×n ) with D ε positive definite on Ω, (∇ · D ε ) · ν = 0 on ∂Ω for all ε ∈ (0, 1) and We can further choose this family in such a way as to ensure that for all Φ ∈ C 0 (Ω; R n ) and ε ∈ (0, 1). These additional properties for the approximation D ε essentially mean that the regularity properties assumed for D are also valid for said approximation in an ε independent fashion.
Remark 4.1. Let us briefly illustrate how such an approximation of D = (d i,j ) i,j∈{1,...,n} can be achieved. This will be a two-step process. We first approximate D in our desired function space with the appropriate boundary conditions and then, as a second step, we show that, with only slight modification, we can gain the remaining properties from that approximation.
For the initial approximation, we assume without loss of generality that D is smooth. We can do this as it is well-known that a standard convolution argument would give us a smooth approximation of D in our desired space, which we can then approximate again to gain all additional desired properties. In our case, the key property not covered by such a convolution based method is that we want all our approximate matrices to have very specific boundary values. As such, we will now demonstrate how an approximation of a smooth D by matrices with exactly this property can be achieved using the continuity properties of semigroups associated with carefully chosen sectorial operators (cf. [12]).
To this end, we fix functions d ′ i,j such that for all i, j ∈ {1, . . . , n}. As can be easily seen, the functions d ′ i,j are linear combinations of the components of D and therefore smooth as well. We then set d ′ i,j,ε = e εLi,j d ′ i,j , ε ∈ (0, 1), where L i,j is the negative Laplacian on Ω with boundary conditions ∇ϕ·ν+ 1 2 (∂ xi ϕ)ν j + 1 2 (∂ xj ϕ)ν i = 0 and (e tLi,j ) t≥0 is the associated semigroup. Due to the well-known continuity properties of said semigroup (cf. [15], [22], [33]), we know that d ′ i,j,ε → d ′ i,j and therefore d i,j,ε → d i,j in W 1,q (Ω) ∩ C 0 (Ω) as ε ց 0 with d i,j,ε defined in an analogous fashion to (4.4). Thus, D ε := (d i,j,ε ) i,j∈{1,...,n} → D in our desired way. Further, on ∂Ω for all ε ∈ (0, 1) due to the prescribed boundary conditions of the operators L i,j . Thus, we have constructed a suitable approximate family for D with the correct boundary conditions.
Having now presented the full argument used to achieve the boundary condition (4.5), let us briefly note that we introduced the functions d ′ i,j to ensure that the operators L i,j have sufficiently non-tangential boundary conditions and are therefore sectorial (cf. [22], [33]), which is of course necessary for our semigroup based arguments.
As our second step, we will now fix one such family of approximations of D and call it D ′ ε , ε ∈ (0, 1), as we still want to slightly modify it. We can assume that for all ε ∈ (0, 1) without loss of generality. If we then set D ε : for all ε ∈ (0, 1) without affecting any of the desired properties that we already derived as we only modify D ′ ε by adding constants that converge to zero as ε ց 0. This gives us (4.3). To derive the divergence estimate, we first observe that for all ε ∈ (0, 1) and Φ ∈ C 0 (Ω; R n ). We can then further estimate for all ε ∈ (0, 1) and Φ ∈ C 0 (Ω; R n ) using our assumed divergence estimate for D and (4.3). This gives us (4.2) and thus completes the discussion of our construction.
We will now proceed to construct our approximate initial data. To do this, we first fix families (u 0,ε ) ε∈(0,1) , (w ′ 0,ε ) ε∈(0,1) ⊆ C 3 (Ω) of positive functions with (D ε ∇u 0,ε ) · ν = (D ε ∇w ′ 0,ε ) · ν = 0 on ∂Ω and as ε ց 0. These families can again be constructed by using convolutions or by a similar semigroup based method as seen before in the much more challenging case of the family (D ε ) ε∈(0,1) . Positivity of both families can further be achieved by first approximating the function in a non-negative way, which is a property of both convolution and semigroup based methods, and then adding ε to the resulting approximation as a secondary step.
One important consequence of the above approximations is that we can fix a uniform constant M > 0 such that and for all ε ∈ (0, 1).

Uniform A Priori Estimates
We will now derive the bounds necessary to ensure compactness of our families of approximate classical solutions in function spaces conducive to the construction of our desired weak solutions to (2.1) as limits of said approximate solutions along a suitable sequence of ε ∈ (0, 1).
Apart from the baseline established in Lemma 3.2 for the classical existence theory, which can be easily translated to our approximate solutions in an ε-independent fashion, we will now derive some extended bounds based on an energy-type inequality as an additional baseline for later arguments in this section. This type of energy inequality was already used in the one-dimensional case in [39].
By another testing procedure for the first equation in (4.8), which is very similar to the one already used by us in the proof of Lemma 4.2, we will now derive our final preliminary set of bounds for this section. Proof. Fix T > 0.
We first note that as ε = ε j ց 0.
Proof. Given that both the families (u (Ω)) * ) according to Lemma 4.5, we can apply the Aubin-Lions lemma (cf. [32]) to the above families using the triple of embedded spaces W 1,2 D (Ω) ⊆ L 1 (Ω) ⊆ (W n+1,2 (Ω)) * . Note that this is only possible as the first embedding is in fact compact by our assumptions (cf. Definition 2.7). Therefore, there exists a null sequence (ε j ) j∈N ⊆ (0, 1) and functionsũ, w : Ω × [0, T ) → R such that as ε = ε j ց 0. This sequence is constructed by applying the Aubin-Lions lemma countably infinitely many times on time intervals of the form [0, T ], T ∈ N, combined with a straightforward extension and diagonal sequence argument. We can further choose the above sequence in such way as to ensure that u 1 2 ε →ũ and w ε → w pointwise almost everywhere as ε = ε j ց 0 by potentially switching to another subsequence. Due to the family (w ε ) ε∈(0,1) furthermore being uniformly bounded in L ∞ (Ω × (0, ∞)) (cf. (4.7) and Lemma 3.2), the above convergence properties directly imply (4.20) as well as the fact that w is non-negative almost everywhere and w ∈ L ∞ (Ω × (0, ∞)).
We now set u :=ũ 2 and observe that the above almost everywhere pointwise convergence for the already constructed sequences then ensures that As all not yet explicitly established regularity properties for u and w directly follow from the convergence properties and we have at this point proven all said properties, this completes the proof.
For the remainder of this section, we will now fix the functions u, w as well as the sequence (ε j ) j∈N constructed in the preceding lemma. While the convergence properties derived in Lemma 4.6 are in fact already sufficient to allow us to translate the weak solution property from our approximate solutions to our now established solution candidates, we will as a last effort before the proof of Theorem 2.12 derive some more specifically tailored convergence properties to handle some of the more complex terms in the weak solution definition.
We now similarly estimate that