Asymptotic Behavior of Solutions of a Free Boundary Problem Modelling the Growth of Tumors with Stokes Equations

We study a free boundary problem modelling the growth of non-necrotic tumors with fluid-like tissues. The fluid velocity satisfies Stokes equations with a source determined by the proliferation rate of tumor cells which depends on the concentration of nutrients, subject to a boundary condition with stress tensor effected by surface tension. It is easy to prove that this problem has a unique radially symmetric stationary solution. By using a functional approach, we prove that there exists a threshold value $\gamma_*>0$ for the surface tension coefficient $\gamma$, such that in the case $\gamma>\gamma_*$ this radially symmetric stationary solution is asymptotically stable under small non-radial perturbations, whereas in the opposite case it is unstable.

Note that the general situation can be easily reduced into this special situation by using the rescaling σ → σ/σ, p → p/ν and γ → γ/ν.
Tumor growth modelling and analysis has attracted considerable attention during the past more than ten years. Most tumor models assume that the tumor tissue has the structure of a porous medium for which Darcy's law applies (see, e.g., [2], [3] and the references cited therein). For such tumor models, many interesting results of rigorous analysis have been obtained, for which we refer the interested reader to see [4]- [8], [18]- [22], [26], [27] and the references cited therein. The tumor whose tissue does not have the structure of a porous medium but instead is more like a fluid was recently considered by Franks et al in the literatures [12]- [15], where some new models were proposed to mimic the early stages of the growth of ductal carcinoma in the breast. A basic feature of a ductal carcinoma in the breast in early stages is that it is confined to the duct of a mammary gland, which consists of epithelial cells, a meshwork of proteins, and extracellular fluid. In modelling, this leads to the replacement of the Darcy's law used in porous medium structured tumor models by the Stokes equations. See [12]- [15] for details. The models of Franks et al [12]- [15] have been concisely reformulated by Friedman in [16] (see also [17]).
The problem (1.1)-(1.9) above is a simplification of the tumor model proposed in [16]. The simplifications are made in two aspects. First, the model in [16] contains a system of nonlinear hyperbolic conservation laws with source terms describing the movements and interchanges of three different species of cells: proliferating cells, quiescent cells and necrotic cells. In this paper we only consider one species of proliferating cells, so that no hyperbolic conservation laws appear in (1.1)- (1.9). Second, in [16] the equation for σ is of the following evolutionary type: where f is as given in (1.10), but in this paper the stationary form (1.1) is considered. All these simplifications are made for the purpose to make the analysis simpler. If either one of the above two aspects of simplifications are not made, then the model will be much more complicated to analyze, and new mathematical techniques have to be employed. We leave it for future work.
In [16] Friedman established local wellposedness in Hölder spaces of his model. Meanwhile, he proved that in the special case that the tumor contains only one species of cells (i.e., the tumor contains only proliferating cells), there exists a unique radially symmetric stationary solution. Based on these results, a number of interesting questions are raised in [16] (see also [17]), one of which is as follows: Is this radially symmetric stationary solution asymptotically stable under non-radial perturbations? A heuristic result toward an answer to this question was obtained by Friedman and Hu in [19], where they proved that this radially symmetric stationary solution is linearly asymptotically stable for small µ/γ, i.e., there exists a threshold value (µ/γ) * such that if we denote by (σ s , v s , p s , Ω s ) this stationary solution, then in the case µ/γ < (µ/γ) * the trivial solution of the linearization at (σ s , v s , p s , Ω s ) of the original problem is asymptotically stable. Moreover, they also proved that in the case µ/γ > (µ/γ) * the radially symmetric stationary solution is unstable. We also refer the interested reader to see [18] for the study of existence of non-radial stationary solutions.
The purpose of this paper is to prove that, at least for the simplified model (1.1)-(1.9), the answer to the above question is yes for large γ but no for small γ. To this end, we shall use a functional approach inherited from the references [7], [10], [26] and [27] to study the problem (1.1)-(1.9), namely, we shall first reduce the problem (1.1)-(1.9) into an evolution equation containing merely the function ρ describing the free boundary ∂Ω(t), which can be considered as a differential equation in certain Banach space. We shall prove that this differential equation is of the parabolic type. Next we use the geometric theory for parabolic differential equations in Banach spaces (see [1] and [24]) to study the stability of the stationary solution. Since our discussion does not depend on the specific linear forms of the equations (1.1) and (1.2), throughout this paper we shall not consider the specific forms of f and g given by (1.10), but instead assume that they are general smooth functions satisfying the following assumptions: To give a precise statement of our main result, let us first introduce some notations.
Given a bounded domain Ω ⊆ R 3 and two numbers m ∈ N and θ ∈ (0, 1), we denote by h m+θ (Ω) the so-called little Hölder space on Ω of index m + θ, which is, by definition, the closure of C ∞ (Ω) in the usual Hölder space C m+θ (Ω). Similarly, given a smooth hypersurface Γ in R 3 , we denote by h m+θ (Γ) the closure of C ∞ (Γ) in C m+θ (Γ).
It can be easily shown (see Theorem A in Appendix A) that under Assumptions (A1)-(A3), the problem (1.1)-(1.8) has a unique radially symmetric stationary solution. Later on we use the same notation (σ s , v s , p s , Ω s ) as before to denote this radially symmetric stationary solution of the problem (1.1)-(1.8). Note that this means that there exists R s > 0 such that where r = |x| and v s represents a scalar function. Clearly, a coordinate translation of a solution of (1.1)-(1.8) is still a solution of it. Thus, for any the stationary solution obtained by the coordinate translation x → x + x 0 of the stationary solution (σ s , v s , p s , Ω s ). Given ρ ∈ C 1 (∂Ω s ) with ρ C 1 (∂Ωs) sufficiently small, we denote by Ω ρ the domain enclosed by the hypersurface r = R s + ρ(ω), where ω ∈ ∂Ω s . It is obvious that for Since we shall only be concerned with small perturbations of the stationary solution (σ s , v s , p s , Ω s ), it is natural to assume that the domains Ω(t) and Ω 0 in (1.1)-(1.9) are small perturbations of Ω s . It follows that there exist functions ρ(t) (= ρ(ω, t)) and ρ 0 (= ρ 0 (ω)) on ∂Ω s such that Ω(t) = Ω ρ(t) and Ω 0 = Ω ρ 0 . Using these notations, the initial condition (1.9) can be rewritten as follows: The solution (σ, v, p, Ω) of the problem (1.1)-(1.9) will be correspondingly rewritten as (σ, v, p, ρ), and the radially symmetric stationary solution (σ s , v s , p s , Ω s ) will be re-denoted as (σ s , v s , p s , 0).
The main result of this paper is the following theorem: Theorem 1.1 Assume that Assumptions (A1)-(A3) hold. For given m ∈ N, m ≥ 3, and 0 < θ < 1, we have the following assertion: There exists a positive threshold value γ * such that for any γ > γ * , the radially symmetric stationary solution (σ s , v s , p s , 0) is asymptotically stable in the following sense: There exists constant ε > 0 such that for any ρ 0 ∈ h m+θ (∂Ω s ) satisfying ρ 0 C m+θ (∂Ωs) < ε, the problem (1.1)-(1.9) has a unique solution (σ, v, p, ρ) for all t ≥ 0, and there are positive constants ω, K independent of the initial data and a point x 0 ∈ R 3 uniquely determined by the initial data, such that the following holds for all t ≥ 0: (1.14) For γ < γ * the stationary solution (σ s , v s , p s , 0) is unstable. 2 It is interesting to compare this result with the corresponding result for the porous medium structured tumor model obtained by Cui and Escher in [8], where it is proved that, for the porous medium structured tumor model, there exists a threshold value for the surface tension coefficient γ, which we denote asγ * , such that the unique radially symmetric stationary solution is asymptotically stable if γ >γ * , but unstable if γ <γ * . We shall show that γ * >γ * . This implies that radially symmetric stationary solution is more stable for a tumor whose tissue has a porous medium structure than a tumor whose tissue is more like a fluid. See Lemma 3.5 for the proof of the assertion that γ * >γ * .
The structure of the rest part is as follows. In Section 2 we first convert the problem into an equivalent initial-boundary value problem on a fixed domain by using the so-called Hanzawa transformation, and next further reduce it into a scalar equation containing the single function ρ, which can be regarded as a differential equation in the Banach space h m−1+θ (S 2 ). We shall also prove that this equation is of the parabolic type. In Section 3 we study the linearization of (1.1)-(1.8) at the radially symmetric stationary solution, and study the spectrum of the linearized operator. In the last section we give the proof of Theorem 1.1.

Reduction of the problem
In this section we reduce the problem (1.1)-(1.9) into a differential equation in a Banach space. For simplicity of the notation, later on we always assume that R s = 1. Note that this assumption is reasonable because the case R s = 1 can be easily reduced into this case after a rescaling. It follows that Let m and θ be as in Theorem 1.1. For u ∈ h m+θ (B 3 ), we denote by tr S 2 (u) the trace of u on S 2 , i.e., tr S 2 (u) = u| S 2 . We know that tr S 2 (u) ∈ h m+θ (S 2 ) and the operator tr is linear, bounded and surjective. Let Π be a bounded right inverse of it, i.e., Π ∈ L(h m+θ (S 2 ), h m+θ (B 3 )) and tr S 2 Π(u)) = u for any u ∈ h m+θ (S 2 ).
) be an extension operator, i.e., E has the property that Here BU C m+θ (R 3 ) denotes the space of all C m functions u on R 3 such that u itself and all its partial derivatives of order≤ m are bounded and uniformly θ-th order Hölder continuous in R 3 . We denote . Hence there exists a constant C 0 > 0 such that Given ρ ∈ O m+θ δ (S 2 ), we define the Hanzawa transformation Φ ρ : R 3 → R 3 as follows: Using (2.1) and (2.2) we can easily verify that Φ ρ is a h m+θ diffeomorphism from This implies that x ∈ Γ ρ if and only if there exists ω ∈ S 2 such that x = [1 + ρ(ω)]ω. Thus, in the polar coordinates (r, ω) of R 3 , where r = |x| and ω = x/|x|, the hypersurface Γ ρ has the following equation: r = 1 + ρ(ω).
Since our purpose is to study asymptotical stability of the radially symmetric stationary solution, later on we always assume that the initial domain Ω 0 is contained in a small neighborhood of Let ρ be as above, and let Φ i ρ be the i-th component of Φ ρ , i = 1, 2, 3. We denote where ∇ ω represents the orthogonal projection of the gradient ∇ x onto the tangent space T ω (S 2 ) 1) . Here and hereafter, for a matrix A we use the notation A ij to denote the element of . We now introduce four partial differential operators A(ρ), B(ρ), B(ρ)· and B(ρ)⊗ on R 3 as follows: Here and hereafter we use the convention that repeated indices represent summations with respect to these indices, and ∂ j = ∂/∂x j , j = 1, 2, 3. Obviously, The definitions (2.4)-(2.7) can be respectively briefly rewritten as follows: As in [16] we introduce the following vector functions: Let n and κ be respectively the unit outward normal and the mean curvature of Γ ρ (see (1.5)). We denote A direct computation shows that We remind the reader to notice that, obviously, (2.10) Finally, for σ, v and p as in (1.1)-(1.8), we denote We also denote w ρ j = w j • Φ ρ , j = 1, 2, 3. Using these notations, we see easily that the Hanzawa transformation transforms the equations (1.1)-(1.5), (1.7) and (1.8) into the following equations, respectively: In what follows we rewrite (1.6) into an explicit equation expressed with the function ρ = ρ(ω, t).
It follows that the normal velocity of Γ ρ (t) is as follows (see [10]): Hence (1.6) can be rewritten as follows: we see that after the Hanzawa transformation, this equation has the following form: Finally, we rewrite (1.13) as follows: In summary, we have the following preliminary result: In the sequel we further reduce the problem (2.11)-(2.19) into a scalar equation containing the single unknown ρ. The idea is to first solve the system of equations (2.11)-(2.17) to get σ, v and p as functionals of ρ, and next substitute v obtained in this way into the equation (2.18).

2
Next, for given ρ ∈ O m+θ δ (S 2 ) we consider the following boundary value problem: Lemma 2.3 Let δ be sufficiently small and let ρ ∈ O m+θ δ (S 2 ) be given. A necessary and sufficient condition for (2.22)-(2.26) to have a solution is that ϕ, g and h satisfy the following relations: If this condition is satisfied, then (2.12)-(2.26) has a unique solution Proof: Integrating by parts and employing the divergence theorem we see that for any v, w ∈ (C 2 (Ω ρ )) 3 and p ∈ C 1 (Ω ρ )) there holds the following integral identity (cf. [11]): (2.29) Here, for two matrix A and B we use the notation A · B to denote the inner product of A and B, i.e., A · B = A ij B ij . By making the Hanzawa transformation, we see that for any v, w ∈ (C 2 (B 3 )) 3 and p ∈ C 1 (B 3 ) there holds (2.30) Besides, clearly ∇ ⊗ w j + (∇ ⊗ w j ) T = 0, ∇ · w j = 0, j = 1, 2, 3, which yield, after the Hanzawa transformation, that Hence, if ( v, p) is a solution of (2.22)-(2.26), then by takeing w = w ρ j (j = 1, 2, 3) in (2.30), we see that (2.27) holds. Similarly we have (2.28). This proves the necessity of (2.27) and (2.28).
Next we assume that the conditions (2.27) and (2.28) are satisfied, and proceed to prove that there exists a unique solution to the problem (2.22)-(2.26). We first prove uniqueness of the solution. Let ( v 1 , p 1 ) and ( v 2 , p 2 ) be two solutions of (2.22)-(2.26). Then v = v 1 − v 2 and p = p 1 − p 2 satisfy the corresponding homogeneous equations. Thus, by letting w = v in (2.30), we get This combined with (2.25) and (2.26) yields, by the Korn inequality [25]), that v = 0. From this it follows immediately that also p = 0. Hence the solution is unique if it exists.
From these relations we can easily show that if ρ = 0 then ζ = 0. By continuity (a small perturbation of a nonsingular matrix is still nonsingular), this implies that if δ is sufficiently small then for any ρ ∈ O m+θ δ (S 2 ) we also have ζ = 0. Hence our claim is true. It follows that ( v, p) is a solution of (2.22)-(2.26). This completes the proof of Lemma 2.3.
Proof: Let notations be as in the proof of Lemma 2.3. We denote by I 1 , I 2 and I 3 the natural embedding operators from h m−k−1+θ (B 3 ), (h m−k−2+θ (B 3 )) 3 and (h m−k−1+θ (S 2 )) 3 , respectively, into Y, and by J the projection operator from X onto (h m−k+θ (B 3 )) 3 . Then by letting we immediately see that the desired assertion follows.
We summarize: , p, ρ) be a solution of the problem (2.11)- (2.19). Then ρ is a solution of the initial value problem (2.39). Conversely, if ρ is a solution of the initial value problem (2.39), then by letting σ = R(ρ) and ( v, p) = P L(ρ) −1 (g( σ), 1 3 In the sequel we prove that if δ is sufficiently small then for any ρ ∈ O m+θ δ (S 2 ), DQ(ρ) is a infinitesimal generator of an analytic semigroup in h m−1+θ (S 2 ) with domain h m+θ (S 2 ), so that the differential equation in (2.39) is of parabolic type. Here and in what follows, the notation D· represents Fréchet derivatives of smooth operators from h m+θ (S 2 ) to h m−1+θ (S 2 ).
We first note that the mean curvature operator K has the following expression K(ρ) = K 1 (ρ)ρ + K 0 (ρ), (2.41) where, for each ρ, K 1 (ρ) is a second-order linear elliptic partial differential operator on S 2 with coefficients being functions of ρ and its first-order derivatives, and K 0 is a first-order nonlinear partial differential operator on S 2 , so that ). (2.42) (see [7] and [10]). Substituting (2.41) into (2.38) we see that where, for each ρ, Q 2 (ρ) is a bilinear operator, Q 1 (ρ) is a linear operator, and Q 0 is a nonlinear operator; they are respectively defined as follows: We note that, by (2.21), (2.37), (2.42) and Lemma 2.4 we have , h m−1+θ (S 2 ))), where BL(·× ·, ·) denotes the Banach space of all bilinear operators with respect the corresponding spaces.
Given two Banach spaces E 0 and E 1 such that E 1 is continuously and densely embedded into E 0 , we denote by H(E 1 , E 0 ) the subset of all linear operators A ∈ L(E 1 , E 0 ) such that −A generates a strongly continuous analytic semigroup on E 0 .
Besides, a simple computation shows that L(0)(∇ψ η , 0, 0)= 0, 0, T(∇ψ η , 0)n, It follows that It can be easily seen that B 0 ∈ L(h m+θ (S 2 ), h m+θ (S 2 )). Furthermore, minor changes to the proof of Lemma A.2 in [9] show that also B 1 ∈ L(h m+θ (S 2 ), h m+θ (S 2 )) (see Lemma B.1 and Corollary B.2 in Appendix B for details). Hence (2.49) follows. This completes the proof of Lemma 2.6. 2 (2.50) By this corollary we see that, at least in a small neighborhood of the origin, the differential equation (2.39) is of the parabolic type in the sense of Amann [1] and Lunardi [24], so that the geometric theory for parabolic differential equations in Banach spaces presented in these literatures can be applied to (2.39). In the following sections we shall use this theory to prove Theorem 1.1.

Linearization
In this section we compute the spectrum of the operator DQ(0). Note that since h m+θ (S 2 ) is compactly embedded into h m−1+θ (S 2 ), by Lemma 2.6 we see that the spectrum of the operator DQ(0) consists of all eigenvalues.
In summary, we have proved the following result: where α 0 and α l (γ) defined in (3.49) and (3.48), respectively. 2 As usual, for a given closed linear operator B in a Banach space X, we denote by ρ(B) and σ(B) the resolvent set and the spectrum of B, respectively. The set of all eigenvalues of B is denoted by σ p (B). As mentioned in the beginning of this section, we have σ(DQ(0)) = σ p (DQ(0)). Hence, from Lemma 3.2 we immediately obtain the following result: Moreover, the multiplicity of the eigenvalue 0 is 3. 2 The next result shows some useful properties of α 0 and γ l (l ≥ 2): Lemma 3. 4 We have the following assertions: (i) α 0 < 0.
From this fact the assertion (iii) immediately follows. The proof is complete. 2 By virtue of the the assertion (ii) of the above lemma, we introduce γ * = max l≥2 γ l .