INITIAL BOUNDARY VALUE PROBLEMS FOR A SYSTEM OF PARABOLIC CONSERVATION LAWS ARISING FROM CHEMOTAXIS IN MULTI-DIMENSIONS

. We study the qualitative behavior of a system of parabolic conservation laws, derived from a Keller-Segel type chemotaxis model with singular sensitivity, on the unit square or cube subject to various types of boundary conditions. It is shown that for given initial data in H 3 (Ω) , under the assumption that only the entropic energy associated with the initial data is small, there exist global-in-time classical solutions to the initial-boundary value problems of the model subject to the Neumann-Stress-free and Dirichlet-Stress-free type boundary conditions; these solutions converge to equilibrium states, determined from initial and/or boundary data, exponentially rapidly as time goes to inﬁnity. In addition, it is shown that the solutions of the fully dissipative model converge to those of the corresponding partially dissipative model as 2010 Mathematics Subject Classiﬁcation. 35A01, 35B40, 35K61, 35L65, 35Q92, 92C17. the chemical diﬀusion rate tends to zero under the Neumann-Stress-free type boundary conditions. Numerical analysis is performed for a discretization of the model with the Dirichlet-Stress-free type boundary conditions, and a mono- tonic exponential decay to the equilibrium solution (analogous to the continuous case) is proven. Numerical simulations are supplemented to illustrate the exponential decay, test the assumptions of the exponential decay theorem, and to predict boundary layer formation under the Dirichlet-Stress-free type boundary conditions.

the chemical diffusion rate tends to zero under the Neumann-Stress-free type boundary conditions. Numerical analysis is performed for a discretization of the model with the Dirichlet-Stress-free type boundary conditions, and a monotonic exponential decay to the equilibrium solution (analogous to the continuous case) is proven. Numerical simulations are supplemented to illustrate the exponential decay, test the assumptions of the exponential decay theorem, and to predict boundary layer formation under the Dirichlet-Stress-free type boundary conditions. 1. Introduction. We consider the following system of parabolic conservation laws: where Ω = (0, 1) n , n = 2, 3. The purpose of this paper is to study the global wellposedness, long-time behavior and diffusion limit (as ε → 0) of classical solutions to initial-boundary value problems for (1) on the unit square or cube subject to appropriate boundary conditions.
1.1. Background. The system of parabolic conservation laws, (1), is derived from the Keller-Segel type chemotaxis model with singular sensitivity: which appeared as a sub-model in [26] to understand the interaction between vascular endothelial cells (VEC) and vascular endothelial growth factor (VEGF) associated with tumor angiogenesis. Here, the unknown functions u(x, t) and v(x, t) denote the density of VEC and concentration of VEGF, respectively. The parameters D > 0 and ε ≥ 0 are the diffusion coefficients of VEC and VEGF, respectively, and the function µu with µ > 0 measures the degradation rate of VEGF. It is worth mentioning that comparing with D, the magnitude of ε could be small or negligible, due to the observation that the diffusion of VEGF is far less important than its interaction with VEC (cf. [26]). The logarithmic (singular) sensitivity function, ln v describing the signal detection mechanism of VEC with the constant χ > 0 measuring the strength of chemotactic sensitivity, is incorporated to follow the Weber-Fechner's law which states that subjective sensation is proportional to the logarithm of the stimulus intensity. We remark that the logarithmic sensitivity has also appeared in the original Keller-Segel model [22] designed to understand the chemotactic wave propagation of bacterial population observed in the experiment by Adler [1] in addition to its prominent applications in biology (cf. [7,21]), and in some other chemotaxis models describing aggregation (cf. [37,46,12] and references therein). Phenomenologically, since χ > 0, the first equation of (2) indicates that while naturally diffusing, the population of VEC is driven against diffusion by the concentration gradient of VEGF at spatial locations where the chemical signal increases, which suggests that the population of VEC could aggregate over time at certain spatial locations. On the other hand, because the degradation rate of VEGF is proportional to the density of VEC, the signal strength driving the aggregation of VEC decreases rapidly when the population of VEC accumulates. Hence, the density-dependent degradation rate acts effectively to help prevent the population of VEC from aggregating over time. This provides a possible explanation for the role of protease inhibitors in preventing angiogenesis, which, when being applied to angiogenic areas, accelerate the speed of degradation of VEGF and stop the growth of new blood vessels from pre-existing ones, as proposed and numerically observed in [26]. Such a feature, together with the singular sensitivity function, makes the model (2) biologically relevant and mathematically intriguing.
Since the time when the model (2) was initiated, attentions have been attracted from the community of nonlinear PDEs to rigorously study its qualitative behavior. Early evidences, such as the explicit and numerical solutions constructed in [25,26], have led researchers in the field to appreciate that (2) is capable of describing some of the underlying mechanisms elucidating the process of preventing angiogenesis, such as the tendency against aggregative pattern formation leading to the uniform distribution of VEC. However, from the rigorous mathematical perspective, in contrast to a large number of results on the chemotaxis models with linear sensitivity function (see review papers [3,16,17]), existing results for (2) are much less, due to the technical obstacles caused by the logarithmic singularity at v = 0. The well-adapted idea for overcoming such singularity is to take the Hopf-Cole type transformation ( [25]): Q = ∇ ln v, which converts (2) into the system of parabolic conservation laws: x ∈ R n , t > 0. Further applying the scaling to (3) and letting P ≡ u, we obtain the following system of equations :    ∂ t P − ∇ · (PQ) = ∆P, where the tilde is dropped for simplicity. Throughout this paper, we are concerned with the qualitative behavior of the model for fixed values of χ and D. To simplify the presentation, we take χ = D = 1, resulting in the following system of equations : are real when p ≥ 0. This indicates that the principal part of (1) is hyperbolic in biologically relevant regimes. Such a feature enables researchers to adapt the analytic techniques for hyperbolic PDEs, such as symmetrization, to study the qualitative behavior of (2). Due to the competition between cellular aggregation and chemical degradation, as illustrated in (2), it is natural to ask: "Aggregation versus degradation, who wins?". It is well known in the mathematical analysis of chemotaxis models that when aggregation dominates over degradation, the formation of finite-time singularities is oftentimes detected in solutions, while the uniform distribution of cellular (e.g. VEC) population is observed when degradation prevails. In the one-dimensional case, it has been confirmed that the model (2) possesses global-in-time classical solutions for any finite value of χ > 0 and any classical initial data, which indicates that the destabilizing action of chemo-attraction indeed is dominated by diffusion and degradation, see the references in Section 1.2. However, the persistence of such a feature in the multi-dimensional spaces is unknown, which is the main question to be answered in this paper. It is worth mentioning that the dominance of chemoattraction by diffusion and/or degradation in multi-dimensional spaces has been partially verified for the related model: where χ > 0, for which the global-in-time classical solutions (cf. [4,46,11,23,52]) and generalized/weak solutions (cf. [41,24]) have been constructed, and the asymptotic stability of homogeneous states has been studied (cf. [49]), at least for small values of χ. But the system (2) has a very significant difference from the system (5). It can be shown that v-component in (5) has a positive lower bound if the initial value of v is positive and hence the singularity is naturally ruled out (see [11]), however the system (2) has no such nice property and the singularity become a real challenging issue. This paper aims to enrich the knowledge base by providing partial answer to the above question for the model (2) by studying the qualitative behavior of multi-dimensional solutions to the transformed equations (1).
To get a rough idea about the qualitative behavior of (1), let us first take a look at the scaling invariant property enjoyed by the model. Indeed, by a direct calculation, we can show that (1) holds its form under the scaling (p, q) → (p λ , q λ ) := λ 2 p(λ 2 t, λx), λq(λ 2 t, λx) .
Under this scaling, when the initial data are perturbed around the zero ground state, it holds that q λ 0 2 L 2 (R n ) = λ 2−n q 0 2 L 2 (R n ) , which reveals that norm-inflation (especially for the q-component) is only possible when n = 1.
Next, we note that the weak Lyapunov functional associated with (1) reads : under appropriate boundary conditions, wherep > 0 is a constant ground state and the "entropy expansion" is defined by It should be mentioned that the weak Lyapunov functional, (6), has been observed and extensively utilized in the literature for studying the qualitative behavior of (2), see e.g. [28,29,30,35,45,47]. Because of the scaling property of the q-component and the fact that the right hand side of the weak Lyapunov functional is zero only when n = 1, from the point of view of energy criticality we then see that the global well-posedness of (1) is sub-critical when n = 1, critical when n = 2, and super critical when n ≥ 3. Such an observation partially explains why the model is globally well-posed in one space dimension, as was observed in many previous works (see below), while the problem is still widely open in the multi-dimensional case.
1.2. Literature review. To put things into perspective, we first briefly survey the literature in connection with this work. Since the qualitative behavior of the transformed model, (1), differs in one-and multi-dimensional spaces, we list related references in two categories.
In the one-dimensional space, due to the energy sub-criticality, the following results concerning the global well-posedness and long-time behavior of large amplitude classical solutions to (1) have been established: • global well-posedness of Cauchy problem when ε = 0 [13], • global well-posedness of initial-boundary value problems when ε = 0 [10], • global well-posedness and long-time behavior of initial-boundary value problems when ε ≥ 0 [29,30,42,45], • global well-posedness and long-time behavior of Cauchy problem when ε ≥ 0 [28,35].
The results for the one-dimensional version of (1), especially those reported in [28,29,30,35,42,45], indicate that no matter how strong the chemotactic sensitivity is and how large the energy of initial data is, the population of VEC always distributes uniformly over space as time evolves. Hence, the aforementioned speculation that the density-dependent degradation rate acts effectively to prevent the population of VEC from aggregating is confirmed. From the point of view of biological applications, the results provided explicit formulas for the convergence rates of solutions toward equilibrium states in terms of the system parameters and initial and boundary data, which can bring useful information to biologists to quantify the time required for a cellular population to become uniformly distributed over space, at which an in vitro process of angiogenesis stops. There are also works investigating other properties of the one-dimensional version of the model, such as the formation of shock waves [43], asymptotic behavior of small amplitude solutions [51], and nonlinear stability of traveling waves [20,31,32,33,34]. Since they are not directly related to this paper, we shall not discuss them in detail.
In multi-dimensional spaces, due to the energy criticality and super-criticality, the qualitative behavior of the model has been investigated relatively little, comparing with the magnitude of research conducted in the one-dimensional case. The following results have been recently established: • local well-posedness and blowup criteria of classical solutions in R n when ε = 0 [9,27], • global well-posedness and long-time behavior of small classical solutions in R n when ε = 0 [14,27], • global well-posedness and long-time behavior of small classical solutions on bounded domains in R n when ε = 0 [30], wherep is a positive background state, • global well-posedness and long-time behavior of classical solutions in R 3 when ε ≥ 0 and (p 0 −p, q 0 ) L 2 (R 3 ) is small [39], • global well-posedness and long-time behavior of strong solutions in R n when ε ≥ 0 and (p 0 −p, q 0 ) H 1 (R n ) is small [44]. • Global generalized (weak) solutions on bounded domains in R 2 with Neumann boundary conditions [47], followed with a work addressing the eventual smoothness of solutions [48].

1.3.
Motivations. Now we would like to point out the facts that motivate the current work. First, we note that the results reported in [8,39,44] improved previous ones in [14,27,30] by relaxing the smallness assumptions on the initial data from the whole spectrum to lower frequencies. However, after a close inspection of the results reported in [8,39,44] we find that none of them gives a positive answer to the question of global well-posedness of solutions to (1) when initial data contain potentially large L 2 -norm of the perturbation. The first goal of this paper is to strengthen the previous results by further relaxing the constraints on the initial data (especially the smallness assumptions). We aim to establish the global wellposedness and long-time behavior results for initial data potentially having large L 2 energy. Specifically, we consider initial data having the following behavior : where η(p) = p ln p − p, andp is the equilibrium state associated with the specific problem under consideration. For simplicity, we shall call solutions with such type of initial conditions small entropic solutions. We remark that assuming the smallness of the first order expansion of the anti-logarithmic function allows p 0 −p H 2 to be arbitrarily large. Please see Remark 1. Second, due to biological consideration, we study (1) on bounded domains in R n with physical boundaries, and construct small entropic solutions to several initialboundary value problems for the model. In this case, another interesting question arises, which is related to the vanishing coefficient limit in a PDE system. It was mentioned in the beginning of the previous section that comparing with the diffusion coefficient, D, of VEC, the magnitude of the diffusion coefficient, ε, of VEGF can be negligible (cf. [26]). Hence, from the point of view of biological modeling, it is desirable to know whether the chemically non diffusible (ideal) model is a good approximation of the diffusible (realistic) one when ε is small. Mathematically, this is related to the question of under what types of boundary conditions solutions with ε > 0 converge to those with ε = 0, as ε → 0. Such a question also arises naturally from the numerical point of view, due to the smallness of ε can cause stiffness in the discrete equations, which can substantially increase the computational cost for simulating the diffusible model. Hence, knowing the conditions under which the two models are consistent helps provide useful information for both the modeling and numerical purposes.
Recently, it has been confirmed that in one space dimension, the chemically diffusible and non diffusible models are consistent when the spatial domain has no physical boundary [35] or on finite intervals (say [a, b]) subject to the Neumann-Dirichlet boundary conditions: ∂ x p| x=a,b = 0, q| x=a,b = 0 [29,45]. On the other hand, the Dirichlet-type boundary conditions: p| x=a,b =ū ≥ 0, q| x=a,b = 0, create boundary layers in the process of vanishing chemical diffusion limit [18,19,29,38].
Although the aforementioned results identified the boundary conditions under which the chemically non diffusible model is or is not a good approximation of the diffusible one, they are restricted to the one-dimensional case. A search in the database shows that the zero chemical diffusion limit of multi-dimensional solutions to (1) has not been studied. This is partially due to the geometrical complications of multi-dimensional bounded domains with physical boundaries, which make characterizing the boundary conditions under which the consistency between the diffusible and non diffusible models may or may not hold a challenging issue. To contribute to this contemporary body of knowledge, in this paper, we consider (1) on bounded domains with physical boundaries in R n (n > 1), prescribe different types of boundary conditions, and study the global well-posedness, long-time behavior and diffusion limit of small entropic solutions to the initial-boundary value problems.
1.4. Boundary conditions. The first set of boundary conditions is prescribed as follows : which we refer to as the Neumann-Stress-free boundary conditions. Here, Ω = (0, 1) n , n = 2, 3, n is the unit outward normal to ∂Ω, and τ i denotes the unit tangential vector on ∂Ω when n = 3. We would like to mention that the Neumann-Stress-free boundary conditions are naturally inspired by the Neumann (no-flux) boundary conditions for the original model (2). Indeed, the stress-free boundary conditions for the function q are consistent with the Neumann boundary condition for the original function v in (2), due to the transformation: q = ∇v/v. We will see later that in addition to the global well-posedness and long-time behavior of small entropic solutions, a direct outcome of prescribing the Neumann-Stress-free boundary conditions is the consistency between the chemically diffusible and non diffusible models in the limiting process as ε → 0. We further remark that the geometrical configuration of the rectangular domain, i.e. Ω = (0, 1) n , plays a crucial role in establishing the vanishing chemical diffusion limit result, due to the boundary conditions in (8), together with the equations in (1), yield (desirably) (∆q) · n| ∂Ω = 0, without which the result can not be obtained. To the best of our knowledge, the vanishing chemical diffusion limit of solutions to (1) on physical spatial domains with non-flat boundaries subject to the Neumann-Stress-free boundary conditions is still a significant mathematical challenge. Second, inspired by recent studies concerning the boundary layer phenomena associated with the one-dimensional version of (1) (cf. [18,19,29,38]), we prescribe the following Dirichlet-Stress-free boundary conditions: Unlike the case of Neumann-Stress-free boundary conditions, although we can still manage to prove the global well-posedness and long-time behavior of small entropic solutions to (1) associated with the Dirichlet-Stress-free boundary conditions, no diffusion limit result in this case can be established. This is mainly due to the mismatch between the boundary conditions for the chemically diffusible and non diffusible models. Therefore, we conjecture that boundary layers emerge from solutions to the IBVPs of (1) subject to the Dirichlet-Stress-free boundary conditions. Our numerical experiments, see Section 6, partially confirm such a speculation, and we leave the rigorous justification in a forthcoming paper.
2. Statement of results and ideas of proof.
2.1. Neumann-Stress-free boundary conditions. The first result of this paper addresses the global well-posedness and long-time behavior of classical solutions to (1) under the Neumann-Stress-free boundary conditions when ε > 0. In the sequel, Ω = (0, 1) n , n = 2, 3, unless otherwise specified.
In addition, the following decay estimate holds: where the positive constants α 1 , β 1 are independent of t and remain bounded as ε → 0, and β 1 > 0 for all ε ≥ 0. Remark 1. We remark that by assuming the smallness of the entropy expansion, i.e., the spatial integral of the first order expansion of the anti-logarithmic function, the usual L 2 energy of the initial datum, p 0 −p L 2 , can be potentially large, especially when the equilibrium statep > 0 is large. This can be seen from the following inequality: Indeed, let us consider the function It is straightforward to check that (p) 2 and F (w) > 0 for w ∈ [0,p 4 ), it holds that F (w) > 0 for w ∈ 0,p 4 . Therefore, F (w) ≥ 0 for all w ∈ [0, ∞). In the Appendix, we provide an explicit example of initial functions whose p-component can have arbitrarily small entropic energy, but potentially large L 2 energy.
The next result is concerned with the global well-posedness and long-time behavior of classical solutions to (1) under the Neumann-Stress-free boundary conditions when ε = 0.
In addition, the following decay estimate holds: where the positive constants α 2 , β 2 are independent of t.
The third theorem establishes the consistency and convergence rate between the chemically diffusible and non diffusible models under the Neumann-Stress-free boundary conditions. Theorem 2.3. Let (p ε , q ε ) and (p 0 , q 0 ) be the solutions obtained in Theorem 2.1 and Theorem 2.2, respectively, with the same initial data. Then for any t > 0, there hold that where the constants γ 1 and γ 2 are independent of t and remain bounded as ε → 0.

Define
where C Ω > 0 denotes a generic constant depending only on Ω. Then there exists a generic constant 2 > 0 depending only on Ω, such that for any 0 < ≤ 2 , if then there exists a unique solution to (1) & (9), such that In addition, the following decay estimate holds: where the positive constants α 3 , β 3 are independent of t and remain bounded as ε → 0, and β 3 → 0 as ε → 0.
We have several remarks in order.
Remark 2. In the theorems above, we assumed q 0 to be curl free, due to the consideration that the system of parabolic conservation laws (1), reduces to the transformed system (4) when q is curl free, and the latter one originates from the Keller-Segel model (2), through the transformation q = ∇ ln v. In Section 3.1 we show that under the initial curl free condition, the system (1) automatically produces curl free solutions (potential flows) under the prescribed boundary conditions, which allows one to recover the solutions to the original chemotaxis model.

Remark 3.
We remark that although the results reported in this paper are oriented around classical solutions to (1) for initial data in the H 3 space, the reader will see from the proofs that the energy estimates indeed can be closed at the H 2 level. Hence, one can also adopt the energy estimates derived in this paper to study the qualitative behavior of less regular (strong) solutions to the model. In addition, the energy method developed in this paper can be of independent interest, and can be adapted to study other PDE models with similar structures.
Remark 4. Another way to allow solutions to carry potentially large L 2 -energy of the zeroth frequency is to assume small oscillation in the initial data, i.e., assume ∇p 0 2 L 2 and ∇ · q 0 2 L 2 to be sufficiently small. This is particularly interesting in the case of the Cauchy problem, since it is well known that the zeroth frequency of a function can not be dominated by its first order derivative in the whole space. Indeed, the global well-posedness of classical solutions to the Cauchy problem of (1) in multi-dimensional spaces with initial data having only small oscillation is still an open problem in the field. We will give a positive answer to the question in a forthcoming paper. On the other hand, in the case of a fixed bounded domain, because of the Poincaré inequality, the smallness of ∇p 0 L 2 to be small. In this case, the proof will be in the same spirit of this paper and we shall not go through the technical details.
Remark 5. Theorem 2.2 shows that small entropic solutions to the chemically non diffusible model under the Neumann-Stress-free boundary conditions exist globally in time. Indeed, such solutions can be obtained from those of the chemically diffusible problem by taking ε → 0. On the other hand, from the proof of Theorem 2.4 we will see that the zero diffusion limit process can not be realized under the Dirichlet-Stress-free boundary conditions, due to the energy bounds derived under the Dirichlet-Stress-free boundary conditions are inversely proportional to ε. This leads to the question of whether the chemically non diffusible model is globally wellposed under the Dirichlet-Stress-free boundary conditions. However, we note that for the chemically non diffusible model, it makes no sense to prescribe the Dirichlet-Stress-free boundary conditions because otherwise the problem is over-determined. From a rigorous mathematical point of view, only the following boundary value problem makes sense: Nevertheless, we noted that one can not fully estimate the spatial derivatives of q by using classical inequalities (see e.g. [5]), even assuming the curl-free condition. The resolution might require a significant improvement of classical results in analysis. We leave the investigation for the future.

Numerical analysis.
We consider also the numerical approximation of (1) with the Dirichlet-Stress-free boundary conditions, using a finite element in space and backward Euler in time discretization. It is numerically advantageous to shift the p variable to have mean 0, and we apply this. For appropriate finite element spaces R h and Q h (see Section 5 for details), the scheme is defined as Algorithm 2.1. For a given time step size ∆t and initial conditions p 0 ∈ R h and Here, γ ≥ ε is a user chosen parameter to penalize the curl of the discrete variable q n h . Under a smallness assumption on the data, and a time step size restriction, we prove that the algorithm's solutions decay exponentially in time.
, with C being a generic constant independent of the mesh width, time step size, and ε. Then the solution to Algorithm 2.1 decays exponentially fast to 0. Extension of this result to a second order analogue of the scheme is possible, and discussed in a remark.
2.4. Ideas of proof. Now we briefly explain the difficulties of the problem and the main approaches to overcome the technical barriers. First of all, since the entropic energy of p 0 −p (i.e., in Orlicz space) is assumed to be small, the global well-posedness problem considered herein is different from many of the existing results concerning the global well-posedness of nonlinear PDE models for initial data having small energy in Sobolev spaces. Hence, some of the standard approaches for constructing mild solutions are not directly accessible for proving the results recorded in Theorem 2.1. We prove the theorem by combining a priori assumptions, energy estimates and standard continuation argument. The major technical difficulty consists in closing the energy estimates for the individual frequencies of the solution. We achieve our goal by applying various Gagliardo-Nirenberg and Poincaré type inequalities. The key idea is to minimize the application of the elementary inequality: 2ab ≤ a 2 + b 2 , in order to sustain the potential largeness of the individual (higher) frequencies and to avoid using the smallness of the whole spectrum of the solution.
Secondly, since we are also concerned with the zero chemical diffusion limit of the solution, an additional technical difficulty encountered in the proof of Theorem 2.1 is how to obtain the uniform ε-independent energy estimates of the solution. We realize such estimates by deriving a linear and inhomogeneous damping equation for the spatial derivatives of q, taking advantage of the additional boundary conditions obtained from the prescribed ones by exploring the geometry of the spatial domain and the structure of the equations, and again using the Gagliardo-Nirenberg and Poincaré type inequalities with minimal application of the elementary inequality. The results recorded in Theorems 2.2 and 2.3 and their proofs appear to be consequences of the proof of Theorem 2.1, with similar ideas but a considerable amount of elaboration. The proof of Theorem 2.4 also utilizes the similar ideas as Theorem 2.1, with the distinction that one is forced to invoke the classical elliptic theory to recover the energy estimates for the spatial derivatives of the solution through its temporal derivatives, due to the Dirichlet-Stress-free type boundary conditions. A consequence of the proof is that the energy estimates of the solution, especially those for the spatial derivatives, are inversely proportional to ε, which prevents one from getting the zero chemical diffusion limit result under the Dirichlet-Stress-free boundary conditions. For Theorem 2.5, proving that discrete solutions decay exponentially requires a different approach than the continuous case result: since the discrete function spaces contain globally continuous piecewise polynomials, some choices of test functions that were critical in the continuous proof are not viable choices in the discrete case.
The key pieces for the discrete proof are exploiting the extra positivity provided by the particular time derivative approximation (which is controlled with the time step size), using the discrete inverse inequality in combination with a time step size restriction, and a L 2 lower bound on the H 1 discrete solution norms provided by the Poincaré inequality. These first two tools are critical for control of the nonlinear growth terms, and the last is necessary for showing the solution decay (versus nongrowth). These key tools, delicately combined with Hölder's inequality, Sobolev embeddings, and Young's inequality, allowed for exponential decay to be proven for discrete solutions. The proof shows that the exponential decay is monotonic in the norm ( p 2 L 2 + q 2 L 2 ) 1/2 , and our numerical tests show the decay is not monotonic in either of p L 2 or q L 2 .
The rest part of this paper is organized as follows. In Section 3, we focus on the Neumann-Stress-free boundary conditions and study the global well-posedness, long-time behavior and diffusion limit of solutions to (1) under the conditions of Theorems 2.1-2.3. Section 4 is devoted to the study of (1) under the Dirichlet-Stress-free boundary conditions and we prove Theorem 2.4. Section 5 contains a detailed description of the numerical algorithm above (and discussion of higher order schemes), and provides numerical analysis that proves the scheme converges exponentially to the equilibrium state. In Section 6, we carry out numerical simulations for (1) under the Dirichlet-Stress-free boundary conditions. Here, we test the exponential decay results above, test the system behavior when the assumptions of Theorem 2.5 are violated, and investigate the formation of boundary layers subject to the Dirichlet-Stress-free boundary conditions. We note that due to the energy criticality in multi-dimensional spaces, the global well-posedness/finite-time blowup of classical solutions to (1) with arbitrarily large initial data is still elusive from the point of view of rigorous mathematical analysis. On the other hand, the boundary layer formation in multi-dimensional spaces is a very delicate problem and is highly sensitive to the geometry of the spatial domain under consideration. A rigorous mathematical analysis is challenging and involved, which is beyond the scope of the current paper. We leave the investigation of both topics for the future. The paper ends with concluding remarks in Section 7.
3. IBVP under Neumann-Stress-free boundary conditions. In this section, we study the qualitative behavior of (1) under the Neumann-Stress-free boundary conditions and complete the proofs for Theorems 2.1-2.3. Since the threedimensional geometry is more complicated than the two-dimensional one, throughout the rest part of this paper we only deal with the three-dimensional case, in order to simplify the presentation. Moreover, from the proofs we will see that the Poincaré inequality plays an important role in the energy estimates. Since the Poincaré inequality is independent of the spatial dimension, results in the threedimensional case can be automatically transferred to the two-dimensional case. For the convenience of the reader, we use (p, q) to denote the perturbation of the original functions around (p, 0), and restate the IBVP as where ε ≥ 0, Ω = (0, 1) 3 , n is the unit outward normal vector on ∂Ω, and τ i (i = 1, 2) denotes the unit tangential vector on ∂Ω.
Notation. Throughout this section and the rest part of this paper, we denote the usual L 2 norm of a function by · . Unless specified, we use C Ω to denote a generic constant which depends only on Ω, and use C to denote a generic constant which is independent of the unknown functions and t, but depends on the initial data.

3.1.
Observations about boundary conditions. We would like to emphasize that one of the purposes of this section is to investigate the limiting behavior of the solution to (1) when the VEGF diffusion coefficient, ε, tends to zero, under the Neumann-Stress-free boundary conditions. Hence, we first collect some information about the Neumann-Stress-free boundary conditions, which will help us obtain the uniform (independent of ε) a priori estimates of the solution, leading to the desired results concerning the vanishing diffusion limit and convergence rate. The following calculations are achieved by using the geometrical configuration of the spatial domain, Ω = (0, 1) 3 , and the equation for q. In addition, we denote q ≡ (q 1 , q 2 , q 3 ), and denote ω = (ω 1 , ω 2 , ω 3 ) ≡ ∇ × q.
3.1.1. Preservation of curl-free solution. First, we show that the initial curl-free condition, together with the equation for q, produces a curl-free solution. To see this, let us first consider the top surface: For any fixed point x 0 ∈ S T , the subsequent calculations are performed at x 0 . We note that the unit outward normal vector at x 0 is n = (0, 0, 1). Then according to the boundary conditions: q · n| ∂Ω = 0 and ω · τ i | ∂Ω = 0, we deduce that which implies immediately that and, due to the flatness of S T , Since which implies ∂ x2x3 q 1 = 0 = ∂ x1x3 q 2 . Hence, it holds that Moreover, it can be readily checked that the above calculations are also valid on the other surfaces of the domain. We omit the details for simplification. By taking the curl of the second equation of (15), we have By taking the L 2 inner product of the above equation with ω, and using (17) and (19), we have 1 2 Hence, the initial curl-free condition produces a curl-free solution as time evolves. We further notice that although the above calculations regrading the boundary conditions are restricted to flat domains, the conclusion is not. Indeed, we observe that since for any sufficiently regular vector field F ∈ R 3 the following identities hold: the equation for ω, (20), can be written as Hence, for a general domain with non-flat boundary, by taking the L 2 inner product of (21) with ∇ × q, we can show that Due to ω is perpendicular to the unit tangential vectors on the boundary of the domain, we know that ω is parallel to n on ∂Ω. Therefore, the integrand in the surface integral on the right hand side of the above equation vanishes, and we end up with 1 2 which produces the curl-free solution under the initial curl-free condition for general domains with non-flat boundaries. Because of the above observation, it suffices to deal with ω, in order to estimate the spatial derivatives of q. Moreover, as mentioned before, since ∆q = ∇(∇ · q) − ∇ × (∇ × q), under the curl-free condition we have ∆q = ∇(∇ · q).

Slip
Laplacian. Next, we show that (∆q)·n = 0 on the boundary ∂Ω. Again, to simplify the presentation, we only show the result on the top surface S T . The subsequent calculations are performed at any fixed point x 0 ∈ S T . Due to (16) and (18), we know that Since (∂ t q) · n| ∂Ω = 0 and ∇p · n| ∂Ω = 0, by using the equation of q, we get that (ε∆q) · n = ∂ t q − ∇p + ε∇(|q| 2 ) · n = 0.
Therefore, it holds that (∆q)·n = 0 on ∂Ω. Such a property will help us build higher order spatial regularity of the solution by directly differentiating the equations without revoking the elliptic regularity theory, and hence obtain the ε-independent estimate of the solution.
3.1.3. Estimation of spatial derivatives of q. Since q · n = 0 on ∂Ω, the following Poincaré type inequality holds: for some constant which is independent of q. Hence, for the curl-free solution, it holds that On the other hand, since from the classical Poincaré inequality we know that Moreover, since q · n = 0 on ∂Ω, from a classical result (cf. [5]) we know that Hence, for the curl-free solution, it holds that for some constant depending only on Ω and s. When s = 1, we deduce from (24) and (22) that When s = 2, by (25) and (23), we have Moreover, when s = 3, we have due to the Poincaré type inequality, the boundary condition for ∆q and (24). The inequalities, (22) and (25)-(27), will be frequently used in the proofs of our main results.
3.1.4. Estimation of spatial derivatives of p. Lastly, since ∇p · n = 0 on ∂Ω and p has zero mean, from classical elliptic regularity we know that (cf. [5]) for some constant depending only on Ω and s.
In the next subsection, we prove the global well-posedness and long-time behavior of classical solutions to (15) with ε > 0 when the initial data satisfy the assumptions of Theorem 2.1. The results are proved by combining local well-posedness, a priori estimates and continuation argument.

3.2.
Global well-posedness. First of all, we note that since the initial data belong to H 3 (Ω) and the spatial dimension n ≤ 3, and because of the quadratic nonlinearities in the system of equations, the local well-posedness of classical solutions to the initial-boundary value problem (15) for general initial data can be established by applying some of the classical approaches in the literature (e.g., regularization, Galerkin method, fixed point theory, et al; cf. [13,36,40,50]), and carrying out the similar energy estimates as in [27]. Moreover, it follows from the maximum principle that the solution satisfies p +p ≥ 0 within its lifespan. To simplify the presentation, we shall state the results without going through the technical details. for some T 0 ∈ (0, ∞).
Next, we shall derive a priori estimates for the local solution, under the conditions of Theorem 2.1, in order to extend it to be a global one. It should be mentioned that although our analysis is formal, it can be rigorously justified using routine procedures such as regularization, see e.g. [40] and the references therein. First, we observe that according to (10), for the IBVP (15), there hold that where Then it follows from Lemma 3.1 that there exists Next, we derive a priori estimates for the local solution within the time interval [0, T 1 ]. We begin with the estimate of the entropy associated with the system of equations.
3.2.1. Entropy estimate. Denote η(z) = z ln z − z, then we find from the first and second equations of (15) that where the integral on the right hand side can be estimated by using the Hölder inequality as To estimate the term on the right hand side of (33), we recall the Gagliardo-Nirenberg type inequalities for functions defined on a bounded domain in R 3 : which, after being applied to (33), imply where we have applied the Poincaré inequality (25). Hence, when we can update (32) as Hence, in view of (10), we see that This completes the entropy estimate.
Remark 6. We remark that in the two-dimensional case, in the place of (33), we can estimate the nonlinear term by using the 2D Gagliardo-Nirenberg inequality: f L 2 ∇f L 2 and the Poincaré inequality as follows: which is essentially the same as (33).
Next, we carry out regular L 2 -based energy estimates.
3.2.2. L 2 -estimate. By testing the equations in (15) with the targeting functions and applying (34), (22) and (25), we can show that Hence, when δM 1 satisfies (36), there holds that which yields Hence, in view of (10) we see that This completes the L 2 -estimate.

Remark 7.
We remark again that in the two-dimensional case, in the place of (38), we can estimate the nonlinear terms as follows: which is essentially the same as (38).
Next, we estimate the first order spatial derivatives of the solution.  where we applied (34), (26), (28), and the following Gagliardo-Nirenberg inequality: Hence, when we get from (42) that Integrating (45) with respect to time, it holds that Hence, in view of (10) we see that This completes the H 1 -estimate. Remark 9. We further remark that for general domains with non-flat boundaries, it can be readily checked that the a priori estimates presented in Sections 3.2.1-3.2.3 are still valid under the Neumann-Stress-free boundary conditions. However, it is not clear whether or not the energy estimates in the next subsection can be established for general domains, due to the estimates involve the boundary integral of the vector field: ∆p ∇(∇ · q), whose information on the boundaries of general domains is unknown. This is the main reason for which we are unable to obtain the zero chemical diffusion limit and convergence rate results for the model (1) on general domains with non-flat boundaries.
Next, we estimate the second order spatial derivatives of the solution.
3.2.4. H 2 -estimate. By taking the spatial gradient of the first equation and the spatial divergence of the second equation of (15), we get Testing the first equation of (48) with −∇∆p and the second one with −p∆(∇ · q), respectively, and using the boundary conditions: ∇p·n| ∂Ω = 0 and ∇(∇·q)·n| ∂Ω = 0, we get For the first term on the right hand side of (49), by applying the Sobolev embedding: H 2 → L ∞ , and previous ideas, we can show that Similarly, for the second term on the right hand side of (49), we can show that Hence, when there holds that Applying the Gronwall inequality to (51) and using (46), we have Hence, in view of (10) we see that In addition, by plugging (52) into (51), we have By integrating the above inequality with respect to t and using (46), we have (cf.
This completes the H 2 -estimate.

H 3 -estimate.
By using the similar ideas as in deriving the previous energy estimates, we can show that According to (26) and (53), we know that q(t) 2 . Hence, we update the above inequality as By applying the Gronwall inequality to the above inequality and using (46) and (54), we have In view of (30) we see that Moreover, by plugging (56) into (55), we have Integrating (58) with respect to t, we have where which is independent of t and ε. This completes the H 3 -estimate.
3.2.6. Further estimates for q. From previous estimates we see that the temporal integral of the spatial derivatives of q is inversely proportional to ε. This is not desirable, since we also aim to establish the global well-posedness and long-time behavior results for the chemically non diffusible model, i.e., (1) with ε = 0. In this section, we derive the ε-independent temporal integrability of ∇ · q 2 H 2 . Taking the divergence of the second equation of (15), then combining the result with the first equation, we find that ∂ t (∇ · q) +p∇ · q = ε∆(∇ · q) + ∂ t p − ε∆(|q| 2 ) − ∇ · (pq).
Taking the L 2 inner product of (60) with ∇ · q and noting that ∆q = ∇(∇ · q), we have For the first term on the right hand side of (61), we note that where we have used the second equation of (15). Then we update (61) as d dt For the second term on the right hand side of (63), we can show that For the third term on the right hand side of (63), by using (43), we can show that where we have used (26) and the following Gagliardo-Nirenberg type inequality: For the fourth term on the right hand side of (63), we can show that Hence, when we can update (63) as d dt Multiplying (39) by 2, then adding the result to (65), we find where Since the right hand side of (66) is uniformly integrable with respect to time, due to (40), integrating (66) with respect to time then yields, in particular, where the constant is independent of t and remains bounded as ε → 0, due to the presence of the termp 2 ∇ · q 2 on the left hand side of (66). In a similar fashion, we can show that where again the constant is independent of t and remains bounded as ε → 0. The proof is similar to the one presented above, we omit the technical details for simplicity. This completes the proof of the a priori estimates of the local solution, which will be utilized in the subsequent sections to establish the global well-posedness and long-time behavior results recorded in Theorem 2.1.

3.2.7.
Global well-posedness. We now prove the global well-posedness of classical solutions to (15) under the conditions of Theorem 2.1. For this purpose, we first collect the energy estimates from (37), (41), (47), (53) and (57): which, in particular, imply that By applying the local well-posedness, Lemma 3.1, we know that there existsT ∈ (0, ∞), such that there exists a unique classical solution to (15) on the time interval [T 1 , T 1 +T ], which satisfies In view of (31) and (69), we see that By repeating the arguments in Sections 3.2.1-3.2.5, we know that the solution satisfies the same estimates as in (67) for ∀ t ∈ [0, T 1 +T ], and in particular, the estimates in (68) are valid when T 1 is replaced by T 1 +T . Hence, starting from T 1 +T , by applying the local well-posedness again, we know that the solution can be extended byT , such that it satisfies the same uniform energy estimates within the time interval [0, T 1 + 2T ]. By repeating the procedure, we know that the solution indeed exists globally in time. Moreover, by repeating the arguments in Sections 3.2.5-3.2.6, we know that the solution satisfies the same energy estimates derived therein. This completes the proof of the global well-posedness result recorded in Theorem 2.1.

Exponential decay.
In this section, we derive the exponential decay rate of the solution towards the constant equilibrium state by coupling different energy estimates obtained in the previous sections. To illustrate the idea, we only show the decay estimate for the H 1 -norm of the solution. Let By multiplying (39) with 2C 1 , then adding the result to (45) and (66), we have where From the above expressions and Poincaré inequality we can tell that there exists a positive constant C 2 , depending only on ε, Ω,p, such that C 2 > 0 for all ε ≥ 0 and G 1 (t) ≥ C 2 K 1 (t). We update (71) as Hence, the exponential decay of p(t) 2 H 1 + q(t) 2 H 1 follows from (72) and the definition of K 1 (t). It is worth mentioning that the above arguments are valid for all ε ≥ 0. The decay of the higher order derivatives follows in a similar way, and we omit the technical details for simplicity. This completes the proof for Theorem 2.1.
3.4. Chemically non diffusible system. From the energy estimates derived in the previous sections we see that under the conditions of Theorem 2.1, the solution to (15) satisfies where the constant is independent of t and remains bounded as ε → 0. Hence, by essentially applying similar arguments as in the previous sections, we can show that under the conditions of Theorem 2.2, any solution to (12) satisfies for some constant which is independent of t, where (p, q) denotes the perturbation of the solution around (p, 0). In addition, a decay estimate similar to (71) can also be established for the solutions to (12) by repeating the arguments in Section 3.3.
To simplify the presentation, we omit the technical details here. This completes the proof for Theorem 2.2. Next, we quantify the convergence rate between the chemically diffusible and non diffusible models in terms of ε.
Step 1. Taking the L 2 inner products, we find For the first term on the right hand side of (74), by applying Young's inequality, we have where we used the decay estimates recorded in Theorems 2.1-2.2. For the second term on the right hand side of (74), by applying the weighted Young's inequality, we have where again we used the decay estimates recorded in Theorem 2.1. Plugging (75) and (76) into (74), we can show that where β 3 = min{β 1 /2, β 2 }. Applying the Gronwall inequality to (77), we have where the constant on the right hand side is independent of t and remains bounded as ε → 0. Next, we consider the convergence of the first order derivatives of the difference.
Step 2. Taking the L 2 inner products of the first two equations in (73) with the −∆ of the targeting functions, we deduce For the first term on the right hand side of (79), by applying Young's inequality, we have where the second term on the right hand side can be estimated as where we applied various Gagliardo-Nirenberg inequalities, Sobolev embedding and Poincaré inequality. So we can update (80) as For the second and third terms on the right hand side of (79), in a similar fashion, we can show that Plugging (82) and (83) into (79), we find We remark that at this point one can obtain the convergence rate of the first order spatial derivatives of the difference of the solutions, which is of O(ε), by invoking the decay estimates recorded in Theorems 2.1-2.2, applying the Gronwall inequality to (84) and using (54). However, the energy bound in front of the convergence rate will be an exponentially increasing function of t if one sticks to such a procedure, due to the lack of some uniformly integrable quantity in front of the third term on the right hand side of (84). Next, we shall improve the energy bound to be a uniform constant by utilizing the damping equation for ∇ ·q (c.f. (60)).
Step 3. Taking ∇· to the second equation in (73), then combining with the first equation, we have Taking the L 2 inner product of (85) with ∇ ·q, we have 1 2 For the first term on the right hand side of (86), by using the second equation of (73), we deduce By plugging (87) into (86), we have d dt For the first term on the right hand side of (88), by applying Young's inequality, we have where used similar arguments as those in (81). For the second integral on the right hand side of (88), by using similar arguments as those in (83), we can show that For the last term on the right hand side of (88), we have Plugging (89)-(91) into (88), we find d dt Next, we shall couple different energy estimates together to generate a desired one leading to the uniform (w.r.t. time) convergence rate for the first order spatial derivatives of the difference between the solutions.
Step. 4 Multiplying (77) by 2, then adding the result to (92), we have d dt (93) Multiplying (93) by 2, then adding the result to (84), we find where we applied the decay estimates recorded in Theorems 2.1 and 2.2, and By applying the Gronwall inequality to (94), we deduce where the constant on the most right is independent of t and ε. From the definition of H(t) we see that the first order spatial derivatives of the difference between the solutions converge uniformly (w.r.t. time) to zero as ε → 0 with the same rate (O(ε 2 )) as that for the zeroth frequency. Furthermore, by using the similar strategy and working with the higher order derivatives of the solutions, we can show that for some quantity H 2 (t), which is qualitatively equivalent to (p,q) 2 H 2 , and some positive constant β 4 independent of t and ε, the following holds By using the temporal integrability of ε ∆ 2 q ε (t) 2 , we deduce that which shows that the convergence rate of the second order spatial derivatives of the solution difference is slower than that of the zeroth and first frequencies. This completes the proof for Theorem 2.3.

4.
IBVP under Dirichlet-Stress-free boundary conditions. In this section, we prove the global well-posedness and long-time behavior results recorded in Theorem 2.4 for (1) under the Dirichlet-Stress-free boundary conditions. Let us recall the initial-boundary value problem: where (p, q) denotes the difference between the original solution and the constant state (p, 0) withp > 0. First of all, we note that because of the stress-free boundary conditions and the arguments in Section 3.1.1, again the initial curl-free condition is preserved by the solution component q.

4.1.
Global well-posedness. We prove the global well-posedness result recorded in Theorem 2.4 by using the same strategy as in Section 3.2. Again, the heart of the matter is the derivation of a priori energy estimates for the local classical solution, in order to extend it to be a global one. As in Section 3.2, in what follows we assume for some 0 < T < ∞: where σ, N 2 , N 3 > 0 are the same constants as those defined in Theorem 2.4. First, we observe that under the Dirichlet-Stress-free boundary conditions, the first two steps (entropy and L 2 estimates) in the proof of Theorem 2.1 are still valid for (95). Hence, we still have (37) and (40) for (95), i.e., The major feature that distinguishes the proof of Theorem 2.4 from that of Theorem 2.1 is the estimates of spatial derivatives of the solution. It was mentioned before that because of the boundary conditions in Theorem 2.4, one can not proceed to the higher order estimates by directly working on the spatial derivatives of the solution, as was done in the proof of Theorem 2.1. To resolve the issue, we turn to the estimates of the temporal derivatives of the solution and recover the spatial derivatives by invoking classical elliptic theory. Nevertheless, the reader will see that the energy bounds for the spatial derivatives of the solution depend inversely proportionally on ε, see e.g. (101) and (106). Hence, the process of diffusion limit can not be realized, and the formation of boundary layers is highly expected (see Section 6 for numerical experiments).
4.1.1. H 1 -estimate. Taking the L 2 inner products of the equations in (95) with the ∂ t of the targeting functions and applying Young's inequality, we have Note that by using (43) and (64), we can estimate the last two terms on the right hand side of (98) as Hence, when σ(N 3 ) 3 is smaller than some absolute constant (depending on Ω), there holds that Plugging (99) into (98), we get Integrating (100) with respect to time, we see that where we have used (97). In view of (13) we see that This completes the proof of the H 1 -estimate of the solution. It should be stressed that the energy bound for the first order spatial derivatives of the solution is inversely proportional to ε, which is the major difference between the Neumann-Stress-free and Dirichlet-Stress-free boundary conditions.

H 2 -estimate.
Taking ∂ t to the equations of (95), we get Computing the L 2 inner products with the temporal derivatives of targeting functions, we get 1 2 For the first term on the right hand side of (102), we can show that Similarly, for the second term, we can show that Hence, when σN 2 is smaller than some absolute constant, there holds that Integrating (103) with respect to t and using (101), we have According to the first equation of (95), we know that where the third term on the right hand side can be estimated as follows: Hence, when σN 2 is smaller than some absolute constant, we have which implies ∆q . Hence, when σN 2 is smaller than some absolute constant, there holds that ε ∆q ≤ 2 ∂ t q + 2 ∇p , which implies, by (104), Then, we deduce Hence, in view of (13) we see that In addition, after integrating (103) with respect to time, we have which, together with the equations in (95), implies that for some constant which is independent of t, but is inversely proportional to ε. In a similar fashion we can establish the uniform-in-time H 3 -estimate of the solution, for which we omit the technical details to simplify the presentation. This completes the proof of the a priori estimates of the solution. The global well-posedness then follows from the a priori estimates and similar arguments as in Section 3.2.7.

4.2.
Exponential decay. Again, to simplify the presentation, we only present the proof for the decay estimate of the H 1 -norm of the solution. Similar to the proof in Section 3.3, by coupling (39) (still valid under the Dirichlet-Stress-free boundary conditions) and (100) together, we have where Since C 3 ≥ 8, we can update (107) as where We note that according to the Poincaré inequality, On the other hand, by definition, Then it holds that by which we can update (108) as Since C 4 > 0 is independent of t, it follows that we can show by using (109) that from which and the Poincaré inequality the exponential decay of p(t) 2 follows. It should be stressed that since C 3 → ∞ as ε → 0, it follows that C 4 → 0 as ε → 0. Hence, the decaying of the solution slows down as ε → 0 and becomes invalid when ε = 0, which is in contrast to the decay estimates for the Neumann-Stress-free boundary conditions (cf. Section 3.3). This completes the proof of Theorem 2.4.

5.
Numerical discretization and analysis. We consider now numerical discretizations and associated analysis, for the system (1) with Dirichlet-Stress-free boundary conditions (9) in Ω ⊂ R 2 . We take for simplicity Ω to be convex,p = 1, and rename p = p −p to be a shifted concentration (so that p| ∂Ω = 0), which produces the system with boundary conditions Hence convergence of p to 0 in time corresponds to convergence of the unshifted p to the equilibrium state. We prove that solutions to a backward Euler in time, finite element in space, discretization of (111)-(112) decay exponentially, and the decay is monotonic in the energy norm p 2 + q 2 1/2 . The proof of this result is fairly technical, and although it avoids the use of complicated test functions (e.g. no natural logs in test functions), it uses various discrete Sobolev inequalities and the inverse inequality. It also takes full advantage of the extra positivity coming from the backward Euler time derivative approximation. We give numerical results in the next section which illustrate this result, and show that in different norms (e.g. p ) the decay need not be monotonic. The weak formulation of (111)-(112) reads: Find (p, q) ∈ (R, Q) × (0, T ] satisfying (p t , r) + (pq, ∇r) − (∇ · q, r) + (∇p, ∇r) = 0, ∀r ∈ R, Here, γ ≥ ε > 0 denotes a user-specified penalty parameter. In the continuous case, ∇ × q = 0, and so the choice of γ has no effect on the solution. However, in discretizations, curl-free initial conditions will not produce curl-free solutions at later times (due to discretization error in standard choices of finite elements), and γ ε > 0 will provide for a penalization of non-zero curl. The penalization parameter γ was chosen to be γ = 10, 000 in all of our tests; this choice provided for ∇ × q h ≈ 10 −5 q h , whereas if γ = ε, then ∇ × q h ≈ 10 −1 q h .
Recall that the Poincaré inequality holds in R and Q: there exists a constant C depending only on Ω satisfying r ≤ C ∇r , for every r ∈ R, χ ∈ Q.
We have the following Sobolev bounds: there exists a constant C, depending only on Ω, satisfying r L 4 ≤ C r 1/2 ∇r 1/2 , Consider a regular, conforming triangulation τ h , with h the maximum element diameter, and define finite element spaces as follows: We will assume the mesh is sufficiently regular so that the inverse inequality holds in R h and Q h : Moreover, due to the discrete Agmon inequality, we get the existence of a constant C dependent only on Ω such that [6] We consider a conforming finite element spatial discretization with a backward Euler temporal discretization. For a given, regular, conforming triangulation τ h (Ω) (with maximum element diameter h), define the finite element spaces Q h ⊂ Q and R h ⊂ R by where P 2 (τ h ) denotes the space of globally continuous, element-wise quadratic functions on τ h . The discrete algorithm we use is defined as follows.

5.2.
Discretization and analysis. The numerical algorithm we consider is finite element in space and backward Euler in time, and is defined as follows. A second order in time algorithm is also discussed at the end of the subsection. Throughout our analysis, C will denote a generic constant, possibly changing at each occurance, that is independent of h, ∆t, and .
Algorithm 5.1. For a given time step size ∆t and initial conditions p 0 ∈ R h and Under a smallness assumption on the data, and a time step size restriction, we prove now that solutions to Algorithm 5.1 decay exponentially fast in time.
Theorem 5.1. Assume that ε ≤ 1, the initial condition is sufficient small, q 0 h 2 + p 0 h 2 ≤ Cε, and the time step size is chosen to satisfy 0 < ∆t ≤ Then the solution to Algorithm 5.1 decays exponentially fast to zero: where λ > 0 is independent of ε, ∆t, and h. Since this shows the solution decays in the sense of the squared norms of the solution, then the assumptions will also hold at the next time level (since these involve also the squared norms of the solution), and moreover with the same ∆t. Choose r h = p n h and χ h = q n h in Algorithm 5.1. Apply the polarization identity on the time derivative terms and add the equations together to obtain We consider the first term on the right hand side, and add and subtract p n−1 h to the first component of the inner product, and write the term as We will now majorize the first resulting term of (114) by appling Hölder's inequality to followed by a Sobolev inequality for L 4 , which yields where we applied the inverse inequality in the last step. Finally, apply Young's inequality to the right hand side of (115) to get We will bound the second right hand side term of (114) using similar arguments in Similarly, for the second right hand side term in (113), we obtain the bound Combining (113) with (117)-(118), and using that ∇q n h 2 ≤ ∇ · q n h 2 + ∇ × q n h 2 due to the convex domain assumption, we find that 1 2∆t To majorize the first right hand side term, we first apply Young's inequality, and then add and subtract the gradient terms at the next time level, which yields after Cauchy-Schwarz, where in the last step the inverse inequality was applied. In a similar fashion, the second right hand side term is bounded as Using this in (119) yields the bound 1 2∆t With the assumption that the time step size satisfies ∆t ≤ Ch 2 p n−1 h 2 + q n−1 h 2 , we can drop the third and fourth terms from the left hand side, and the bound reduces to 1 2∆t Since ε ≤ 1 is assumed, we reduce further to obtain 1 2∆t Finally using that q n−1 h 2 + p n−1 h 2 ≤ Cε, we get the estimate which after applying the Poincaré inequality implies there exists λ > 0 such that This completes the proof of Theorem 5.1.

Remark 10.
A similar analysis can be applied to the BDF2 analogue of Algorithm 5.1, which can prove that this second order algorithm will also yield solutions that convergence exponentially to zero under similar restriction of the initial data and a similar CFL condition.
6. Numerical experiments. We now present results of several numerical experiments, to illustrate our theory above and to give some insight into the behavior of the chemotaxis system studied herein. We test the system decay for small and large initial data, and the general system behavior on test problems with square and circular domains. All tests in this section were run using Algorithm 5.1; we also ran these same tests using a second order BDF2 analogue of Algorithm 5.1 and found similar results (omitted).
6.1. Numerical Test 1: Exponential decay in time. We consider here a numerical test to illustrate Theorem 2.1; that is, to show that the solution to Algorithm 5.1 decays exponentially in the natural energy norm p 2 + q 2 1/2 . We take ε = 1.0 × 10 −j for j = 0, .., 4, and consider Ω = (0, 1) 2 . The initial conditions are chosen to be and we run Algorithm 5.1 with k = 2 and ∆t = 0.01 to an end time of T = 1, using a 24 × 24 uniform triangular mesh that was further refined near the boundary and provided 16,131 total spatial degrees of freedom at each time step. Plots of solution norms versus time are given in Figure 1, and we observe a monotonic exponential decay in the norm p 2 + q 2 , as expected from Theorem 2.1. Also as expected, the decay is slower with smaller ε. Interestingly, we observe that the decay in the norms p(t) and q(t) is not monotonic, as early in the simulations these norms both slightly increase. 6.2. Numerical Test 2: Blowup for large initial data. The theorems above that prove decay for the chemotaxis system and discretizations rely on small initial data, and so we test here numerically the case of large initial data. We consider the same test problem and discretization parameters as above in Numerical Test 1, but fix ε = 0.01 and change p 0 to be p 0 = −P * sin(πx) sin(πy).
Tests were run with varying P * > 0. We give the values of the energy norms at T = 1 in Table 1, and observe that for P * = 57 the T = 1 solution is O(10 13 ) (we call this blowup) while the P * < 56 solutions remain 'small'. Thus, the value of 57 appears critical for system blowup; several more tests done using P * > 57 were performed, and all exhibited blowup. Since there was a time step restriction for the decay theorem above, we also looked for the critical P * critical by rerunning this test using various time step sizes, results for which are shown in Table 2. Here we observe some dependence of the critical P * on ∆t for larger values, however the P * critical appears to converge as the time step size gets small. From these tests, it is unclear whether the blowup occurs due to numerical instability, or instability of the underlying PDE with large initial data. We leave for future work an in depth numerical study of blowup for this PDE and its discretizations. Our next two numerical tests are motivated by the previous results (cf. [18,19,29]) on the one-dimensional version of (1), where the formation of boundary layer is numerically observed and rigorously justified on a finite interval, I = (a, b), subject to the Dirichlet boundary conditions: p| x=a, x=b =p (> 0), q| x=a, x=b = 0. For the multi-dimensional system under the Dirichlet-Stress-free boundary conditions, (9), although we are currently unable to prove any analytical result, it is worth testing the problem numerically to gain insights into the complex structure of the problem. 6.3. Numerical Test 3: Solution behavior and formation of boundary layer. For our next test, we use the same test problem as Numerical Test 1, to explore the possible existence of boundary layers as well as the general behavior of the solution.
Contour plots of the solutions p and |q| at various times with ε = 0.001 are shown in Figure 2, as well as centerline (x = 0.5) plots of p and q 2 . We observe that p has a rapid initial decay on 0 < t < 0.4, becomes positive (except at the boundary) by t = 0.4, and continues decaying until the end time. The contour plots show that |q| decays less rapidly than p, and suggest the presence of a boundary layer near the walls but not in the corners. The centerline plot for q 2 shows the immediate formation of a boundary layer at the beginning of the simulation, which then slowly decays with time. In Figure 3, a plot of the |q(0.5, y, 0.2)| is given for varying ε, and we observe the boundary layer sharpens as ε decreases. 6.4. Numerical Test 4: A circular domain. For our last test, we take Ω to be the unit circle, and choose p 0 = − sin(π(x 2 + y 2 )) and q 0 = −∇ cos(π(x 2 + y 2 )). The mesh was created by first creating a Delaunay triangulation of the domain, and then refining it near the boundary, which provided for 18,375 total degrees of freedom. The evolution of solutions are displayed in Figure 4 for ε = 0.01, again as contour plots of p and |q|, as well as centerline plots of p and q 2 . We observe that q holds roughly the same pattern as it uniformly decays, while the p variable experiences significant initial transformation and growth before decaying.
We again consider the decay of solutions, and boundary layer behavior, for varying ε. The decay of p and q with time is shown in Figure 5, and we observe that      circle is done with a convex polygon with 100 corners, which could have an effect on the boundary layer if corners play some significant role in this phenomenon.  Figure 5. Shown above are plots of ( p(t) 2 + q(t) 2 ) versus time t (left), and q 2 (0, y, 0.1), for varying ε, for numerical experiment 4.

7.
Conclusion and looking ahead. We have studied the qualitative behavior of a system of parabolic conservation laws, (1), on the unit square or cube subject to various types of boundary conditions. In the case of the Neumann-Stress-free type boundary conditions (8), for naturally prepared initial data with small entropic energy, it is shown that classical solutions to the initial-boundary value problem are globally well-posed for any non-negative value of the chemical diffusion coefficient, and the solutions converge to constant equilibrium states determined by initial and boundary data exponentially rapidly as time goes to infinity, see Theorems 2.1-2.2. Moreover, it is shown that under the Neumann-Stress-free boundary conditions, the chemically diffusible model is consistent with the non diffusible model in the process of zero chemical diffusion limit, see Theorem 2.3.
For the Dirichlet-Stress-free type boundary conditions (9), it is shown that classical solutions to the initial-boundary value problem with small entropic energy are globally well-posed for any positive value of the chemical diffusion coefficient, and the solutions converge to constant equilibrium states determined by initial and boundary data exponentially rapidly as time goes to infinity, see Theorem 2.4. However, the process of zero chemical diffusion limit can not be realized under the Dirichlet-Stress-free boundary conditions, due to the mismatch between the boundary conditions for the chemically diffusible and non diffusible models.
Numerical analysis for the chemotaxis system with the Dirichlet-Stress-free type boundary conditions was given in Theorem 5.1 that complements Theorem 2.4, as it proves that solutions to the discrete problem possess similar decay properties as the continuous problem. This result was illustrated with experiments in Section 6, and confirmed a monotonic (in the energy norm) exponential convergence to the equilibrium state. Extensions of these results to higher order discretizations proved elusive, although numerical results suggest they hold; thus, analysis for the higher order schemes is left as an open problem.
Our numerical simulations exhibit consistency with the rigorous results concerning the decay of the solutions when the initial data are small in certain energy space. More interestingly, they show that boundary layers emerge under the Dirichlet-Stress-free boundary conditions as the chemical diffusion coefficient tends to zero for general initial data, and also that numerical solutions with large initial data may blowup in finite time (however, from our tests it is unclear if the blowup is due to the numerical discretization or instability of the PDE). As was shown in previous works, in the 1D case the solutions always tend to constant equilibrium states no matter how large the initial data are (due to the energy sub-criticality in 1D). On the other hand, as we have seen in this paper that the global well-posedness of multidimensional solutions for arbitrarily large initial data is still unknown because of the energy criticality and super-criticality. For the boundary layer formation, similar to the previous results reported in [18,29] regarding the formation of boundary layers in the one-dimensional version of (1), the singular limit occurs in the q-component, while the p-component behaves regularly in the limiting process. However, a rigorous justification of the boundary layer phenomenon on general multi-dimensional bounded domains is a significant mathematical challenge, due to geometrical complications. We leave the investigation for the future.
Finally, we would like to mention that analyzing the global behavior of (1) in the whole space (Cauchy problem) is considerably more challenging than the bounded domain case, due to the low frequency part of a function can not be dominated by its high frequencies in the whole space. We note that the proof of the main results of this paper rely heavily on the Poincaré inequality on bounded domains. Without such a powerful tool, the energy estimates for the low frequency part of the solution differs significantly from those constructed in this paper. We leave the detailed analysis to a forthcoming paper.