Analysis of a Cahn--Hilliard system with non-zero Dirichlet conditions modeling tumor growth with chemotaxis

We consider a diffuse interface model for tumor growth consisting of a Cahn--Hilliard equation with source terms coupled to a reaction-diffusion equation, which models a tumor growing in the presence of a nutrient species and surrounded by healthy tissue. The well-posedness of the system equipped with Neumann boundary conditions was found to require regular potentials with quadratic growth. In this work, Dirichlet boundary conditions are considered, and we establish the well-posedness of the system for regular potentials with polynomial growth of order less than six and also existence and uniqueness of weak solutions for singular potentials. New difficulties are encountered due to the growth assumptions for the potential, but for regular potentials, we retain the continuous dependence on initial and boundary data for the chemical potential and for the order parameter in strong norms as established in the previous work. Furthermore, we deduce the well-posedness of a variant of the model with quasi-static nutrient by rigorously passing to the limit where the ratio of nutrient diffusion time-scale to the tumor doubling time-scale is small.


Introduction
We consider the following system of equations describing a two component mixture of tumor tissue and healthy (surrounding) tissue in the presence of a chemical species acting as nutrient for the tumor, in Ω × (0, T ), ϕ(x, 0) = ϕ 0 (x), σ(x, 0) = σ 0 (x) in Ω. (1.1e) Here T > 0 denotes a fixed time and Ω ⊂ R 3 denotes a bounded domain with Lipschitz boundary Γ. The system (1.1) describes a diffuse interface model for tumor growth via a Cahn-Hilliard equation coupled with a reaction-diffusion equation for a nutrient species, whose concentration we denote as σ. The order parameter which distinguishes the two components is denoted by ϕ, where the region {ϕ = 1} corresponds to the tumor phase and the region {ϕ = −1} corresponds to the healthy tissue phase.
Here Ψ denotes a potential with two equal minima at ±1, ε > 0 is a parameter relating to the interfacial thickness, γ > 0 denotes the surface tension, λ p , λ a , and λ c are nonnegative constants denoting the rate of proliferation, apoptosis, and consumption of the nutrient, respectively. The positive mobility D(ϕ) corresponds to the diffusivity of the nutrient, and h(ϕ) is an interpolation function such that h(−1) = 0 and h(1) = 1, with the simplest example being h(ϕ) = 1 2 (1 + ϕ). Here, κ ≥ 0 is a constant such that, for a constant mobility D(ϕ) = D 0 , the ratio κ D 0 represents the ratio of nutrient diffusion time-scale and the tumor doubling time-scale. Finally, χ and η are non-negative constants representing chemotaxis (movement of tumor cells towards regions of high nutrients) and active transport (establishment of persistent nutrient concentration gradient in the vicinity of the tumor interface), respectively. We refer the reader to [9] for more details regarding the derivation of the model and the modeling of the chemotaxis and active transport mechanisms.
We point out that the well-posedness of (1.1) with κ = 1, h(⋅) ∈ C 0 (R) ∩ L ∞ (R), homogeneous Neumann boundary conditions for µ and ϕ, and Robin boundary condition for σ has been studied by the authors in [8]. In contrast, here we consider Dirichlet boundary conditions for the following reasons.
• In [8], the class of admissible potential Ψ is restricted to potentials with at most quadratic growth. Here, we are able to extend our class of potentials to have polynomial growth of order up to (but not including) 6 in dimension d = 3 (see (2.1) below).
It is possible to consider more general boundary conditions ϕ ∞ for ϕ, but in this work we restrict to the case ϕ ∞ = −1 which allows for the physical interpretation that the tumor region is enclosed by the healthy tissue. Using the Yosida approximation and due to a priori estimates which does not depend on the Yosida parameter, we can extend our analysis to potentials of the form Ψ(y) =β(y) + Λ(y), (1.3) whereβ ∶ R → [0, ∞] is a convex, proper (i.e., not identically ∞), lower semicontinuous function and Λ ∶ R → R is a C 1,1 perturbation with at most quadratic growth. It is possible thatβ does not possess a classical derivative, and so Ψ ′ may not be well-defined. But we can consider the subdifferential ∂β as a notion of generalized derivative.
Using the notation β(x) ∶= ∂β(x), we denote β 0 (x) to be the element of the set β(x) such that f for x ∈ D(β).
The well-posedness theory of (1.1) will allow us to deduce the existence and uniqueness of weak solutions to in Ω, (1.5e) where ψ ∈ β(ϕ) = ∂β(ϕ) denotes a selection, as β(ϕ) can be multivalued. In particular, if we consider Ψ to be a potential of double-obstacle typê in Ω × (0, T ), (1.8a) The source term h(ϕ)(σ − µ) appearing in (1.8a) and (1.8c) is motivated by linear phenomenological constitutive laws for chemical reactions (see [14] for more details), and is different compared to our choice of source terms in (1.1). The well-posedness of the system has been established in [5,7] for large classes of potentials Ψ and nonlinearities h. Another class of models that describes tumor growth uses a Cahn-Hilliard-Darcy system where v denote a mixture velocity, p denotes the pressure, and S is a mass exchange term. The existence of strong solutions in 2D and 3D have been studied in [16] for the case S = 0. For the case where S ≠ 0 is prescribed, existence of global weak solutions and unique local strong solutions in both 2D and 3D can be found in [15]. We also mention the work of [2] which treats a related system known as the Cahn-Hilliard-Brinkman system. The structure of this paper is as follows. In Section 2, we state the assumptions and the results for (1.1), (1.2) and (1.5). In Section 3 we perform a Galerkin procedure to deduce existence of weak solutions to (1.1), and show further regularity and continuous dependence on initial and boundary data. In Section 4 we pass to the limit κ → 0 in (1.1) and deduce the existence result for (1.2), and then we prove further regularity properties and continuous dependence on initial and boundary data. The regular weak solutions to (1.1) satisfies an energy inequality and this will allow us to employ the Yosida approximation to prove the existence of weak solutions to (1.5), which is done in Section 5.
For convenience, we use the notation L p ∶= L p (Ω) and W k,p ∶= W k,p (Ω) for any p ∈ [1, ∞], k > 0 to denote the standard Lebesgue spaces and Sobolev spaces equipped with the norms ⋅ L p and ⋅ W k,p . In the case p = 2 we use H k ∶= W k,2 with the norm ⋅ H k , and the notation H 1 0 and H −1 to denote the spaces H 1 0 (Ω) and its dual H −1 (Ω). We recall the Poincaré inequality for H 1 0 : There exists a constant C p > 0 that depends only on Ω such that We also state the Sobolev embedding H 1 ⊂ L r , r ∈ [1, 6], for dimension d = 3: There exists a constant C s > 0 depending on Ω and r such that We will also use the following Gronwall inequality in integral form (see [8,Lemma 3.1] for a proof): Let α, β, u and v be real-valued functions defined on [0, T ]. Assume that α is integrable, β is non-negative and continuous, u is continuous, v is non-negative and integrable. If u and v satisfy the integral inequality then it holds that (1.11)

Main results
In this section we state the main results on existence, regularity, uniqueness and continuous dependence first for regular potentials and then for singular potentials. We also state the results for the quasi-static limit κ → 0 for both cases. The results are stated for dimension d = 3, but similar results also holds for d = 1, 2.
(A1) λ p , λ a , λ c , χ and η are fixed non-negative constants, while κ, γ and ε are fixed positive constants. (A2) The functions D, h belong to the space C 0 (R) ∩ L ∞ (R), and there exist positive for some constant k 1 > 0, and exponent p ∈ (1, 5). (A4) The initial and boundary data satisfy We point out that by (2.1) the potential Ψ has at most polynomial growth of order up to (but not including) 6. By the Sobolev embedding H 1 ⊂ L 6 in three dimensions, a consequence of (2.2) is that Ψ(ϕ 0 ) ∈ L 1 . Definition 2.1. We call a triplet of functions (ϕ, µ, σ) a weak solution to (1.1) if, such that, for all ζ, λ, ξ ∈ H 1 0 and for a.e. t ∈ (0, T ), where ⟨⋅, ⋅⟩ denotes the duality pairing between H 1 0 and its dual space H −1 .
for some positive constant C independent of κ, ϕ, µ and σ.
We make the following additional assumptions to prove continuous dependence on initial and boundary data.
Assumption 2.2. In addition to Assumption 2.1, we assume that for some positive constant k 3 and for all s 1 , s 2 ∈ R.
Here, we point out that, the constant C does not depend on κ, and by Theorem 2.2, the criterion ϕ ∈ L 2 (0, T ; H 3 ) is satisfied if Γ is a C 3 boundary. Moreover, Theorem 2.3 provides continuous dependence for the chemical potential µ in L 2 (0, T ; L 2 ) and for the order parameter ϕ in L ∞ (0, T ; L 2 ). This is in contrast to the classical continuous dependence for ϕ in L ∞ (0, T ; H −1 ) one expects for the Cahn-Hilliard equation, compare [7, Theorem 2] and Theorem 2.7 below.
Before we give the result concerning the quasi-static limit κ → 0, we introduce the definition of weak solutions to (1.2).
We point out that we are only able to show H 2 -regularity for ϕ, in contrast with the H 3regularity for the regular potentials. This is due to the fact that the H 3 -regularity estimate is not independent of the Yosida parameter. Before we state the result on continuous dependence, we introduce the inverse Dirichlet-Laplacian operator That is, N (f ) is the unique solution to the Dirichlet-Laplacian problem with right-hand side f . Note that and we can define a norm on H −1 as From this we deduce that (2.14) In particular, ⋅ * and ⋅ H −1 are equivalent norms. Moreover, we obtain from (2.11) the relation for f ∈ H 1 (0, T ; H −1 ) and a.e. t ∈ (0, T ). This holds from the fact that ∂ t f ∈ L 2 (0, T ; H −1 ) and thus N (∂ t f ) is well-defined and satisfies Multiplying the above equality with ζ ∈ C ∞ c (0, T ) and integrating in time, we obtain the following chain of equalities (here (⋅, ⋅) denotes the L 2 scalar product) Theorem 2.7 (Partial continuous dependence and full uniqueness). Under (C1) and (C2) of Assumption 2.2, for any two weak solution quadruples obtained from Theorem 2.6 with boundary conditions µ ∞, In particular, if σ ∞,1 = σ ∞,2 , ϕ 1,0 = ϕ 2,0 and σ 1,0 = σ 2,0 hold, then we have ϕ 1 = ϕ 2 , We note that continuous dependence on the initial and boundary data can only be shown for ϕ and σ. In particular, we do not have control over the differences µ 1 − µ 2 and ψ 1 − ψ 2 . In the proof we have to apply the inverse Dirichlet-Laplacian operator (2.11) and thus it is necessary that µ ∞,1 = µ ∞,2 . Partial continuous dependence are often observed in the case β is multivalued (see for instance [11,Remark 3

.1] and [10, Remark 2.3]).
In the same spirit as Theorem 2.4, we will pass to the limit κ → 0 in (1.5) to deduce the existence of weak solutions for the quasi-static model of (1.5) with κ = 0.
Since the above result is proved by combining the proofs of Theorems 2.4, 2.5 and 2.7, we will omit the proof of Theorem 2.8.

Regular potentials 3.1 Galerkin approximation
We prove Theorem 2.1 via a Galerkin procedure. Let us consider the set of eigenfunctions of the Dirichlet Laplacian {w j } j∈N which are chosen such that they form an orthonormal basis of L 2 and an orthogonal basis of H 1 0 . Let W k = span{w 1 , . . . , w k } ⊂ H 1 0 denote the finite dimensional subspace spanned by the first k eigenfunctions. We introduce the Galerkin ansatz and we require that ϕ k , µ k and σ k satisfy for all 1 ≤ j ≤ k. For convenience we introduce the following matrices where δ ij is the Kronecker delta. Note that M is the identity matrix precisely due to the orthonormality of {w j } j∈N in L 2 . Furthermore, we use the notation , so that when we substitute the Galerkin ansatz into (3.2), we obtain the following nonlinear initial value problem in vector form Note that ψ k , h k , m k h,σ , s k D,σ , M h and S D depend nonlinearly on the solution. We can also express (3.3) as an initial value problem just in terms of α k and τ k by substituting (3.3b) into (3.3a). Moreover, by the continuity of h(⋅) and D(⋅), we observe that the right-hand side of (3.3a) and (3.3c) depends continuously on α k and τ k . Thus, by the Cauchy-Peano theorem we infer the existence of local solutions Next, we derive some a priori estimates to deduce that α k , β k and τ k can be extended to the interval [0, T ] and to allow us to pass to the limit. In the following, the various constants C may vary line to line, but they do not depend on k or on κ.
Then, by Young's inequality and the Poincaré inequality (1.10) applied to f = µ k − µ ∞ and f = ω k , we obtain where we have used Similarly, (3.10) Then, substituting (3.7), (3.8) and (3.10) into (3.6) leads to and applying the Gronwall inequality (1.11) with ∀s ∈ (0, T ], for some positive constant C independent of k and κ. Then, by the Poincaré inequality, (3.9) and (2.3), we find that for some positive constant C independent of k and κ. This a priori estimate in turn guarantees that the solution {ϕ k , µ k , σ k } to (3.3) exists on the interval [0, T ], and thus t k = T for each k ∈ N.
We now provide some a priori estimates on the time derivatives. Let Π k denote the orthogonal projection onto W k = span{w 1 , . . . , w k }. Then, for any ζ ∈ L 2 (0, T ; H 1 0 ), we see that where {ζ k j } 1≤j≤k ⊂ R k are the coefficients such that Π k ζ = ∑ k j=1 ζ k j w j . Thus, from (3.2a), we find that and so (3.12) Integrating in time from 0 to T and using (3.11), we have for some positive constant C independent of k and κ. In a similar manner, from (3.2c) we find that for some positive constant C independent of k and κ.
Fix j and consider δ(t) ∈ C ∞ c (0, T ). Then, δ(t)w j ∈ L 2 (0, T ; H 1 0 ). We multiple (3.2) with δ(t) and integrate in time from 0 to T , leading to Together with the weak convergence σ k ⇀ σ in L 2 (0, T ; L 2 ), we obtain by the product of weak-strong convergence The terms involving D(⋅) can be treated in a similar fashion. From the above compactness results, we have For any q < 6, choose r = 3 2 − 3 q < 1. Then, the Sobolev embedding H r ⊂ L q yields that C([0, T ]; H r ) ⊂ L q (0, T ; L q ) for q < 6. Thus, ϕ k → ϕ strongly in L q (0, T ; L q ) ≅ L q (Ω × (0, T )) for q < 6.
We here also state the version of the generalized Lebesgue theorem used above. Then,

Energy inequality
By the a.e. convergence ϕ k → ϕ in Ω × (0, T ), the non-negativity and continuity of Ψ and Fatou's lemma, we have that for a.e. s ∈ (0, T ], Taking the supremum in s ∈ (0, T ] leads to Moreover, by the following convergences and the weak lower semi-continuity of the Sobolev norms, we find that passing to the limit k → ∞ in (3.11) leads to (2.5).

Further regularity
Suppose that Γ is a C 2 boundary. From (2.4b), we have that which can be seen as the weak formulation of Recall that Ψ ′ satisfies the growth condition Ψ ′ (y) ≤ k 1 (1 + y p ), p ∈ (1, 5), (3.20) and to prove the regularity assertion, we employ a bootstrapping argument. By the Sobolev embedding H 1 ⊂ L 6 we observe that the term µ + χσ in the right-hand side of (3.19a) belongs to L 6 . The bootstrapping argument is as follows. For j ∈ N, we define a sequence of positive numbers {l j } j∈N such that Then, it can be shown that for p < 5 we have and so {l j } j∈N is a strictly increasing sequence. Then, it holds that where the factor 6 6−(5−p)l 1 is strictly greater than 1. This implies that the sequence {l j } j∈N tends to infinity. Furthermore, by the Gagliardo-Nirenburg inequality [6, Theorem 10.1, p. 27], it holds that and so First step. By the growth assumption (3.20), we have H 1 , and due to the boundedness of ϕ in L ∞ (0, T ; H 1 ), the right-hand side of (3.19a) belongs to L 2 (0, T ; L l 1 ) and elliptic regularity theory ([13, Theorem 2.5.1.1] or [12,Theorem 8.13]) yields that ϕ ∈ L 2 (0, T ; W 2,l 1 ).
Remark 3.2. The restriction of r ∈ [1, 6] for the first regularity assertion of Theorem 2.2 is due to the fact that H 3 -regularity requires a C 3 boundary. If Γ is only a C 2 boundary, then at best we can only deduce ϕ ∈ L 2 (0, T ; W 2,6 ) even if the right-hand side of (3.19a) belongs to L 2 (0, T ; H 1 ).

.1 Existence of weak solutions
Let κ ∈ (0, 1] be fixed and (ϕ κ , µ κ , σ κ ) be a weak solution to (1.1) satisfying the hypothesis of Theorem 2.4. Then, we see that the following energy inequality (2.5) is satisfied for where we have used that κ ≤ 1 to obtain that the right-hand side does not depend on κ.

Continuous dependence
The continuous dependence result can be easily obtain by setting κ = 0 in Section 3.4. Given two weak solution triples satisfying the hypothesis of Theorem 2.5, let ϕ, µ, and σ denote their differences, respectively, which satisfy for all ζ, λ, ξ ∈ H 1 0 (Ω). As in Section 3.4, let µ ∞ ∶= µ ∞,1 −µ ∞,2 and σ ∞ ∶= σ ∞,1 −σ ∞,2 denote the difference of boundary data, and we substitute ζ = γεϕ, ξ = Z(σ − σ ∞ ), λ = µ − µ ∞ and λ = Yϕ for positive constants Y, Z yet to be determined. This is equivalent to setting κ = 0 in (3.25). Then, J 1 vanishes and the estimation of J 2 , J 3 , J 4 and J 5 do not depend on κ. One would have to adjust the value for Z accordingly to account for the absence of the first term on the right-hand side of (3.26a), but the same arguments will lead to the assertion of Theorem 2.5.

Maximal monotone operators and the Yosida approximation
In this section, we will briefly review the basic concepts regarding the Yosida approximation of maximal monotone operators. For more details related to maximal monotone operators, subdifferentials and the Yosida approximation, we refer the reader to [17,Chapter III], [3,Chapter III], [22,Chapter 32] and [20, p. 161].
Let T ∶ X → 2 X * be a possibly multivalued mapping with effective domain D(T ) ∶= {x ∈ X ∶ T x ≠ ∅}. We denote the range of T as R(T ) ∶= {y ∈ T x ∶ x ∈ D(T )} and the graph of T as G(T ) ∶= {[x, y] ∈ X × X * ∶ x ∈ D(T ), y ∈ T x}. The mapping T is said to be monotone if ⟨f − g, x − y⟩ ≥ 0 ∀x, y ∈ D(T ), f ∈ T x, g ∈ T y, and T is said to be maximal monotone if and only if It is known that the subdifferential of a convex, proper, lower semicontinuous function is a maximal monotone mapping ([17, Proposition 2.13, p. 124], [18,Theorem 4]). In the case where X is a Hilbert space H, we identify H with its dual H * and let T ∶ H → 2 H be a maximal monotone multivalued mapping with [0, 0] ∈ G(T ), i.e., 0 ∈ T (0). In the Hilbert space setting, we also have an equivalent formulation for the monotone mappings. The mapping T ∶ H → 2 H is monotone if and only if Moreover, for any λ > 0, [17,Corollary 2.9,p. 120] gives that R(I + λT ) = H, where I is the identity operator, i.e., I + λT is surjective. Thus, for any fixed u ∈ H, there exists We define the resolvent J λ ∶ H → D(T ) and the Yosida approximant T λ ∶ H → H of T as respectively. In particular, x = J λ u and f = T λ u. Moreover, we have u = J λ u + λT λ u for all λ > 0 and T λ u = f ∈ T x = T (J λ u).
Then the monotonicity of T and the fact that T λ u ∈ T x = T (J λ u) yield that i.e., the Yosida approximant T λ is a Lipschitz continuous operator with constant 1 λ . By (5.1), the monotonicity of T implies that i.e., J λ is Lipschitz continuous with constant 1.