Regularity theorems for a biological network formulation model in two space dimensions

We present several regularity results for a biological network formulation model originally introduced by D. Cai and D. Hu in {\it Phys. Rev. Lett.}, {\bf 111}(2013), 138701. A consequence of these result is that a stationary weak solution must be a classical one in two space dimensions. Our mathematical analysis is based upon the weakly monotone function theory and Hardy space methods.


Introduction
Let Ω be a bounded domain in R N and T a positive number. Set Ω T = Ω × (0, T ). We study the behavior of weak solutions of the system (1.2) coupled with the initial boundary conditions p(x, t) = 0, m(x, t) = 0, (x, t) ∈ Σ T ≡ ∂Ω × (0, T ), (1.3) m(x, 0) = m 0 (x), x ∈ Ω (1.4) for given function S(x) and physical parameters D, E, γ with properties: (H1) S(x) ∈ L q 0 (Ω), q 0 > N 2 ; and (H2) D, E ∈ (0, ∞), γ ∈ ( 1 2 , ∞). This system was originally derived in ( [12], [13]) as the formal gradient flow of the continuous version of a cost functional describing formation of biological transportation networks on discrete graphs. In this context, the scalar function p = p(x, t) is the pressure due to Darcy's law, while the vectorvalued function m = m(x, t) is the conductance vector. The function S(x) is the time-independent source term. Equation (1.1) can be interpreted as Kirchhoff's law for the flux u ≡ −(I + m ⊗ m)∇p, and the cost is proportional to |(u·∇p)|+c|m| 2γ , where c is a constant [10]. Values of the parameters D, E, and γ are determined by the particular physical applications one has in mind. To give an example, we have γ = 1 2 for blood vessel systems. We would like to refer the reader to [10] for more discussions in this regard.
In general nonlinear problems do not possess classical solutions. A suitable notion of a weak solution must be obtained for (1.1)- (1.4). It turns out [9] that we can introduce the following definition.
A result in [9] asserts that (1.1) -(1.4) has a weak solution provided that, in addition to assuming S(x) ∈ L 2 (Ω) and (H2), we also have Finite time extinction or break-down of solutions in the spatially one-dimensional setting for certain ranges of the relaxation exponent γ was carefully studied in [10]. Further modeling analysis and numerical results can be found in [1]. We also mention that the question of existence in the case where γ = 1 2 is addressed in [10]. In this case the term |m| 2(γ−1) m is not continuous at m = 0. It must be replaced by the following function However, the general regularity theory remains fundamentally incomplete. In particular, it is not known whether or not weak solutions develop singularities in high space dimensions. Recently, Jian-Guo Liu and the author [16] obtained a partial regularity theorem for (1.1)-(1.4). It states that the parabolic Hausdorff dimension of the set of singular points can not exceed N , provided that N ≤ 3.
In this paper, we continue to study the regularity properties of weak solutions. We focus our attention on the case where N = 2. Our main result is: Main Theorem. Let (H1) and (H2) be satisfied, and let (m, p) be a stationary weak solution to (1.1)-(1.3), i.e., the functions m and p are independent of time. Assume that the space dimension N is 2. Then (m, p) is locally a classical solution, provided that S(x) is locally Hölder continuous.
If we apply the proof of the partial regularity theorem in [16] to the situation considered here, we can only conclude that the set of singular points is countable. Even though our estimates are interior ones, we do not foresee any major difficulty in extending our results to the boundary. We encourage the interested reader to try that.
We begin by studying the time-dependent problem. A key observation is that the function p can be decomposed into two pieces in a small neighborhood: The first piece is weakly monotone [17]. A result in [17] asserts that a weakly monotone function in W 1,N loc (Ω) is locally continuous. The proof of this result becomes applicable to our case if N = 2. The second piece is bounded due to a result in [16]. The combination of two is enough to yield the local continuity of p in the space variables. This result is then used to prove that |m| β ∈ L 2 (0, T ; W 1,2 loc (Ω)), |m| β |∇p| 2 ∈ L 1 (0, T ; L 1 loc (Ω)) for each β > 0. Thus in the two-dimensional stationary case, |m| β can be viewed as a BMO function for each β > 0 (see Section 2 for definition and other relevant information). This combined with a result in [14] asserts that (|m| 2 + 1) σ is an A q -weight for each q > 1, σ > 0. Equipped with this, we are able to establish that (m · ∇p) 2 , |∇p| 2 both lie in the local Hardy space, from whence follows the local continuity of m.
To describe the mathematical difficulty involved in our problem, first notice the term (m · ∇p)∇p in (1.2), which represents a cubic nonlinearity. Currently, there is little work done on this type of nonlinearities. Second, the elliptic coefficients in (1.1) satisfy where ·, · denotes the inner product in R N . Since m is not bounded a priori, the largest eigenvalue λ l and the smallest eigenvalue λ s of the coefficient matrix may not satisfy λ l ≤ cλ s .
Here and in what follows the letter c denotes a generic positive number. Thus existing results for degenerate and/or singular elliptic equations such as these in [11] are not applicable. A condition in [2] seems to be satisfied by our elliptic coefficients, but the results there cannot be used to improve the regularity of the terms |∇p| 2 , (m · ∇p) 2 , neither are we able to employ a weighted version of Gehring's lemma [15] in our analysis. Instead, we are motivated by an idea from [3]. See Proposition 3.3 for details.
The rest of the paper is organized as follows. In Section 2, we collect some relevant results about maximal functions, BMO functions, Hardy spaces, and A q -weights. Various regularity results are presented in Section 3. Our main theorem is a consequence of these results.

Preliminary results
In this section, we review some relevant results about maximal functions, BMO functions, Hardy spaces, and A q weights.
Consider the Hardy-Littlewood maximal function M f of a given measurable function f , which is defined by denotes the ball in R N with center at x and radius r, and |B r (y)| is its Lebesgue measure. Also, if no confusion arises, we always suppress the dependence of a function on its independent variables. We refer the reader to [20] for the full story of maximal functions. Here we only mention the well known inequality where · q is the norm in L q (R N ). Note that (2.2) fails when q = 1. Define a class T of normalized test functions on R N by Define the "grand maximal function" f * of a distribution on R N by Here we write φ r for the function r −N φ(r −1 y). Note the similarity between this and (2.1). In , and the Hardy space norm is defined by There is an alternative definition to this that is equivalent and simpler. Specifically, if φ is any C ∞ function on R N with compact support and R N φdx > 0, then f lies in For the purpose of applications to boundary value problems for PDE, we need a local version of the Hardy space.
Let Ω be an open set in R N . We say that a distribution f on Ω lies in H 1 loc (Ω) if for each compact set K 0 ⊂ Ω there is an ε 0 > 0 so that , and the converse is true when f ≥ 0.
This lemma can be found in [20].
It turns out that the inequality We have the set inclusions The following result is contained in [14].
This lemma has played a key role in the proof of our main result. It is well known that there are many useful results for L q spaces when 1 < q < ∞ that fail when q = 1 or q = ∞. This happens, for example, when one is faced with the equation ∆u = g and one wants to have L q estimates for ∇ 2 u in terms of g. Hardy spaces provide alternatives to L q when q = 1 for which there are counterparts to the familiar estimates for 1 < q < ∞. In general BM O is the right substitute for L ∞ . We refer the reader to [20] for more detailed information in this regard. Here we only cite the following result from [20].
. Then u is locally a continuous function on Ω when N = 2.

Main results
In this section, we first develop a couple of regularity results for (1.1)-(1.4) in two space dimensions. Our main theorem is then established a consequence of these results.
The reason our first two propositions in this section are proved for the time-dependent case is not only for the purpose of generality but also due to a private communication to us by P. Markowich stating that numerical experiments for (1.1)-(1.4) in the generality considered here show no signs of singular behavior in solutions. This suggests possible existence of a classical solution to the problem. Unfortunately, the method we have developed here relies on Lemma 2.4, and there is currently no suitable parabolic version of this lemma. As a result, our proof of the main theorem can not be extended to the time-dependent case. Thus it remains to be seen that the numerical evidence mentioned earlier can be verified analytically.
Even though the elliptic coefficients in (3.1) may not be bounded above, we can easily infer from the proof of Lemma 2.3 in [21] that this problem has a unique solution p 1 = p 1 (x, t) in the following sense: (D4) p 1 ∈ W 1,2 0 (B r (y)), m · ∇p 1 ∈ L 2 (B r (y)); (D5) for each ξ ∈ W 1,2 0 (B r (y)) with m · ∇ξ ∈ L 2 (B r (y)) one has Furthermore, we are in a position to assert from (H1) and Proposition 2.1 in [16] that there is a positive number c = c(N, q 0 ) such that Here and in what follows sup (resp. inf) means ess sup (resp. ess inf). Consequently, p 0 ≡ p − p 1 is the unique solution of the boundary value problem −div [(I + m ⊗ m)∇p 0 ] = 0 in B r (y), (3.5) p 0 (x, t) = p(x, t) on ∂B r (y) (3.6) in the sense of (D4)-(D5) with an obvious modification to the boundary condition. That is, we can decompose p(x, t) into the sum of p 0 (x, t) and p 1 (x, t) on B r (y), or equivalently, Set k l = sup ∂Br(y) p. Then (p 0 − k l ) + ∈ W 1,2 0 (B r (y)) with m · ∇(p 0 − k l ) + ∈ L 2 (B r (y)). Thus we can use it as a test function in (3.5), thereby obtaining p 0 ≤ k l a.e. on B r (y).
In fact, we can further conclude that p 0 is a weakly monotone function [17], i.e., (3.8) sup for each sub-domain Ω ′ of B r (y). By Morrey's inequality on spheres as formulated by Gehring [7], we obtain .
Of course, the above inequality can also be established via an elementary calculus argument. Keeping the preceding estimates in mind, we compute ω r (y) ≡ osc Br(y) p = sup = sup ≤ osc Br(y) p 0 + osc Br(y) p 1 ≤c r ∂Br(y) (3.10) Remember that q 0 > 1. Square both sides of the above inequality, divide through the resulting inequality by r, and integrate to obtain (3.11) Note that ω r (y) is a decreasing function of r. We deduce from the preceding inequality that Obviously, p 0 − p is a legitimate test function for (3.5). Upon using it, we derive This finishes the proof.
The continuity of p in the space variables enables us to derive a local version of Proposition 2.1 in [16]. As we shall see, the key difference is that here β does not have to be small. In fact, we shall establish that |m| β ∈ L 2 (0, T ; W 1,2 loc (Ω)), |m| β |∇p| 2 ∈ L 1 (0, T ; L 1 loc (Ω)) for each β > 0.
Proof. Let K > 0, β > 0 be given and v be defined as in (3.15). For L > K, define Set v L = θ L (|m| 2 ). Fix y ∈ Ω, pick r ∈ (0, dist(y, ∂Ω)), and select a smooth cutoff function ξ : Then the function v β L mξ 2 is a legitimate test function for (1.2). Upon using it, we arrive at 1 2 where ε > 0. In the derivation of the third term above, we have used the fact that (3.19) ∇v L = 0 on the set where |m| 2 > L 2 or |m| 2 < K 2 .
Also observe that ∇m = ∇⊗m, and we still have ∇ 1 2 |m| 2 = ∇mm. Set p y,r (t) = − Br(y) p(x, t)dx. Note that m ⊗ m∇p = (m · ∇p)m. Keep this in mind and use v β L (p − p y,r (t))ξ 2 as a test function in (1.1) to deduce where ε > 0 is given as before. By virtue of (3.19), we have that v β−2 L |m| 2 |∇v L | 2 = v β−1 L |∇v L | 2 . Also, v L ≥ K 2 and max Br(y) (p − p y,r (t)) 2 ≤ σ 2 (r). Choose ε suitably small, multiply through the above inequality by 2E 2 , add the resulting inequality to (3.18), and thereby obtain d dt Br(y) Here we have used the fact that p ∈ L ∞ (Ω T ). This is due to Proposition 2.1 in [16]. In view of Proposition 3.1, lim r→0 σ(r) = 0. We can choose r sufficiently small so that the first term on the right-hand in the above inequality can be absorbed into the similar term on the left-hand side there. Integrating the resulting inequality with respect to t and then taking L → ∞ yield the desired result. The proof is complete.
The core of our development is the following proposition, whose proof is inspired by an argument in [3], based upon important contributions due to Müller [18]. Also see Proposition 2.1 in [5], which has become known as the div-curl lemma. Proof. Note that this proposition would be a trivial consequence of Proposition 2.1 in [5] if we had (m · ∇p)m ∈ L 2 (Ω) 2 due to (1.1). Let K 0 ⊂ Ω be compact. Pick 0 < ε 0 < dist(K 0 , ∂Ω). Fix φ ∈ T with R 2 φdx > 0. For each r ∈ (0, ε 0 ) and each y ∈ K 0 we use 1 where χ Ω is the indicator function of Ω and M (|S(x)| s χ Ω )(y) is the value of the maximal function M (|S(x)| s χ Ω ) at y. The fourth term in (3.22) is the most difficult one to handle. Fix σ > 1. By Proposition 3.2, we have |m| 2σ ∈ W 1,2 loc (Ω). Pick a C ∞ function ξ on R 2 satisfying Obviously, we have |m| 2σ ξ ∈ W 1,2 (R 2 ). We can easily check that |∇(|m| 2 ξ + 1) σ | ∈ L 2 (R 2 ). An application of Poincaré's inequality indicates that ( We estimate from Poincaré's inequality that |∇p| s 1 (|m| 2 ξ + 1) Integrate the above inequality over K 0 and keep in mind the inequality (2.2) and the fact that the exponents 2 s 1 , q 0 s on the right-hand side of the preceding inequality are both bigger than 1 to derive The last step is due to Proposition 3.2. This implies the desired result. The proof is complete.
We are ready to prove the main theorem.
(More general results of this nature can be found in [4].) It immediately follows from (H2) that the right-hand side of (3.30) has the same integrability as |∇p|. Consequently, we can appeal to a local version of the Calderon-Zygmund inequality to derive m ∈ W 2,q loc (Ω) N for each q > 1.
Thus m ∈ C 1,δ loc (Ω) N for some δ ∈ (0, 1). Now we are in a position to invoke the classical Schauder estimates for (1.1). To be specific, we can derive from a local version of Theorem 6.13 in [8] that p ∈ C 2,δ loc (Ω). It can easily be inferred from (H2) that the last term in (3.30) is locally Hölder continuous. Use the Schauder estimates for (3.30) to get m ∈ C 2,δ 0 loc (Ω) N for some δ 0 ∈ (0, 1). The proof is complete.