The Impact of Time Delay in a Tumor Model

TL;DR

This paper analyzes how time delay influences the stability and size of tumors in a mathematical model, revealing that increased delay can enlarge tumors and that stability depends on tumor aggressiveness.

Contribution

It introduces a tumor growth model with time delay, proving stability conditions and showing how delay affects tumor size and stability thresholds.

Findings

Existence of a unique stable stationary tumor solution for any tumor aggressiveness.

Stability depends on a critical aggressiveness parameter, with instability for higher values.

Time delay increases tumor size and has a greater impact with higher tumor aggressiveness.

Abstract

In this paper we consider a free boundary tumor growth model with a time delay in cell proliferation and study how time delay affects the stability and the size of the tumor. The model is a coupled system of an elliptic equation, a parabolic equation and an ordinary differential equation. It incorporates the cell location under the presence of time delay, with the tumor boundary as a free boundary. A parameter $\mu$ in the model is proportional to the "aggressiveness" of the tumor. It is proved that there exists a unique classical radially symmetric stationary solution $(\sigma_*, p_*, R_*)$ which is stable for any $\mu > 0$ with respect to all radially symmetric perturbations (c.f. \cite{delay1}). However, under non-radially symmetric perturbations, we prove that there exists a critical value $\mu_*$, such that if $\mu<\mu_*$ then the stationary solution $(\sigma_*, p_*, R_*)$ is linearly stable; whereas if $\mu>\mu_*$ the stationary solution is unstable. It is actually unrealistic to expect the problem to be stable for large tumor aggressiveness parameter, therefore our result is more reasonable. Furthermore, we established that adding the time delay in the model would result in a larger stationary tumor, and if the tumor aggressiveness parameter is larger, then the time delay would have a greater impact on the size of the tumor.

Equations

$$\lambda {\sigma }_t=\Delta \sigma -\sigma \text{in the tumor region }\Omega (t),$$

$$S=\mu [\sigma (\xi (t-\tau ;x,t),t-\tau )-\tilde{\sigma }],$$

$$\left\{\begin{aligned} &\frac{\mathrm{d}\xi}{\mathrm{d}s}=\vec{V}(\xi,s),\quad t-\tau\leq s\leq t,\\ &\xi\big{|}_{s=t}=x.\end{aligned}\right.$$

$$-\Delta p=\mu [\sigma (\xi (t-\tau ;x,t),t-\tau )-\tilde{\sigma }].$$

$$V_n=-\nabla p\cdot n=-\frac{\partial p}{\partial n}\text{on }\partial \Omega (t).$$

$$\sigma =1\text{ on }\partial \Omega (t),$$

$$p=\kappa \text{ on }\partial \Omega (t),$$

$$\Omega (t)={\Omega }_0-\tau \le t\le 0,$$

$$\sigma (x,t)={\sigma }_0(x)\text{ in }{\Omega }_0,-\tau \le t\le 0.$$

$$\lambda {\sigma }_t-\Delta \sigma +\sigma =0,x\in \Omega (t),t>0,$$

$$-\Delta p=\mu [\sigma (\xi (t-\tau ;x,t),t-\tau )-\tilde{\sigma }],x\in \Omega (t),t>0,$$

$$\displaystyle\left\{\begin{array}[]{lr}\displaystyle\frac{\mathrm{d}\xi}{\mathrm{d}s}=-\nabla p(\xi,s),\hskip 20.00003ptt-\tau\leq s\leq t,\\ \xi=x,\hskip 20.00003pt\hskip 20.00003pt\hskip 20.00003pt\hskip 20.00003pts=t,\\ \end{array}\right.$$

$$\sigma =1,x\in \partial \Omega (t),t>0,$$

$$p=\kappa ,x\in \partial \Omega (t),t>0,$$

$$V_n=-\frac{\partial p}{\partial n},x\in \partial \Omega (t),t>0,$$

$$\Omega (t)={\Omega }_0,-\tau \le t\le 0,$$

$$\sigma (x,t)={\sigma }_0(x),x\in {\Omega }_0,-\tau \le t\le 0.$$

$$\partial \Omega (t):r=R_{\ast }+\epsilon {\rho }_0(\theta ),$$

$$\partial \Omega (t):r=R_{\ast }+\epsilon \rho (\theta ,t)+O({\epsilon }^2),$$

$$\sigma (r,\theta ,t)={\sigma }_{\ast }(r)+\epsilon w(r,\theta ,t)+O({\epsilon }^2),$$

$$p(r,\theta ,t)=p_{\ast }(r)+\epsilon q(r,\theta ,t)+O({\epsilon }^2).$$

$$\begin{array}[]{rcl}|\rho(\theta,t)-(a_{1}\cos(\theta)+b_{1}\sin(\theta))|\leq Ce^{-\delta t},\hskip 20.00003ptt>0,\end{array}$$

$$\lambda {\sigma }_t(r,t)-\Delta \sigma (r,t)+\sigma (r,t)=0,\text{in }{B}_{R(t)},t>0,$$

$$-\Delta p(r,t)=\mu [\sigma (\xi (t-\tau ;r,t),t-\tau )-\tilde{\sigma }],\text{in }{B}_{R(t)},t>0,$$

$$\displaystyle\left\{\begin{array}[]{lr}\displaystyle\frac{\mathrm{d}\xi}{\mathrm{d}s}=-\frac{\partial p}{\partial r}(\xi,s),\hskip 20.00003ptt-\tau\leq s\leq t,\\ \xi=r,\hskip 20.00003pt\hskip 20.00003pt\hskip 20.00003pt\hskip 20.00003pts=t,\\ \end{array}\right.$$

$$\frac{\partial \sigma }{\partial r}(0,t)=0,\sigma (R(t),t)=1,$$

$$\frac{\partial p}{\partial r}(0,t)=0,p(R(t),t)=\frac{1}{R(t)},$$

$$\frac{\mathrm{d}R}{\mathrm{d}t}=-\frac{\partial p}{\partial r}(R(t),t),$$

$$R(t)=R_0,-\tau \le t\le 0,$$

$$\sigma (r,t)={\sigma }_0(r),0\le r\le R_0,-\tau \le t\le 0,$$

$$R^{\prime }(t)$$

$$\displaystyle=\frac{\mu}{R(t)}\Big{[}\int_{0}^{R(t)}\sigma(\xi(t-\tau;r,t),t-\tau)r\mathrm{d}r-\int_{0}^{R(t)}\widetilde{\sigma}r\mathrm{d}r\Big{]}.$$

$$\left\{\begin{aligned} &\frac{\mathrm{d}}{\mathrm{d}s}\Big{(}\frac{\partial\xi}{\partial r}\Big{)}=-\frac{\partial^{2}p}{\partial r^{2}}(\xi,s)\frac{\partial\xi}{\partial r},\quad t-\tau\leq s\leq t,\\ &\frac{\partial\xi}{\partial r}\Big{|}_{s=t}=1,\end{aligned}\right.$$

$$\mathrm{d}r=\mathrm{d}r^{\prime}\exp\Big{\{}\int_{t-\tau}^{t}-\frac{\partial^{2}p}{\partial r^{2}}(\xi(s;r,t),s)\mathrm{d}s\Big{\}}=(1+O(\tau))\mathrm{d}r^{\prime}.$$

$$\left\{\begin{aligned} &\frac{\partial\xi}{\partial s}=-\frac{\partial p}{\partial r}(\xi,s),\\ &\xi\big{|}_{s=t}=0.\end{aligned}\right.$$

$$\left\{\begin{aligned} &\frac{\partial\xi}{\partial s}=-\frac{\partial p}{\partial r}(\xi,s),\\ &\xi\big{|}_{s=t}=R(t).\end{aligned}\right.$$

$$r=r^{\prime }+{\int }_{t-\tau }^t-\frac{\partial p(\xi (c;r,t),c)}{\partial r}\mathrm{d}c=r^{\prime }+O(\tau ),$$

$$R^{\prime }(t)$$

$$\displaystyle=\frac{\mu}{R(t)}\Big{[}\int_{0}^{R(t-\tau)}\sigma(r,t-\tau)r\mathrm{d}r-\int_{0}^{R(t)}\widetilde{\sigma}r\mathrm{d}r+O(\tau)\Big{]}.$$

$$-\Delta {\sigma }_{\ast }(r)+{\sigma }_{\ast }(r)=0,r<R_{\ast },$$

Full text

1. Introduction

Over the last few decades, an increasing number of PDE models describing solid tumor growth in forms of free boundary problems have been proposed and studied. All these models provided a better and deeper understanding of the tumor growth. The basic reaction-diffusion tumor model was studied in Greenspan [28, 29], Cui and Escher [9, 8], Escher and Matioc [13], Bazaliy and Friedman [3, 2], Friedman and Hu[22, 21], Friedman and Reitich [26, 27]. Furthermore, the basic model can be extended to more sophisticated ones by adding different factors. For example, Byrne and Chaplain [5], Cui [7], Cui and Friedman [10], Wang [37], Wu and Zhou [40, 39], Xu et. al. [43, 41, 42] analyzed the tumor growth under the effect of inhibitor; Friedman, Hu and Kao [18, 19, 20, 24] considered a multiscale tumor model by adding cell cycle; and Friedman et. al. [15, 25, 30] added the effect of angiogenesis. See also [1, 6, 17, 23, 34, 35, 36, 32, 33, 38] for other extensions to a variety of different tumor models.

One biological meaningful extension of the basic tumor model is to add the effect of time delay $\tau$. In real life, time delays can arise everywhere, since every process, whether it is long or short, would consume time. Time delays can represent gestation times, incubation periods, transport delays, or can simply lump complicated biological processes together, accounting for the time required for these processes to complete. The basic tumor model can be viewed as an approximation of model with time delay, since time delay $\tau$ is rather small compared with the time range $[0,T]$ we consider. However, compared with the basic model, model with time delay are more accurate and consistent with real life.

Here we propose a tumor growth model with time delay in cell proliferation. The time delay is reflected between the time at which a cell commences mitosis and the time at which the daughter cells are produced (it takes approximately 24 hours). In this model, oxygen and glucose are viewed as nutrients, with its concentration $\sigma$ satisfying the reaction-diffusion equation

[TABLE]

where $-\sigma$ is the nutrients consumed by the tumor. Since the diffusion rate of oxygen or glucose (e.g. $\sim\text{1 }\text{min}^{-1}$) is much faster than the rate of cell proliferation (e.g. $\sim\text{1 }\text{day}^{-1}$), $\lambda$ is very small and can sometimes be set to be 0 (quasi-steady state approximation).

By conservation of mass, cell proliferation rate $S=\mathrm{div}\vec{V}$, where $\vec{V}$ denotes the velocity field of cell movement within the tumor. Due to the presence of time delay, the tumor grows at a rate which is related to the nutrient concentration when it starts mitosis. For a simple approximation, we assume a linear relationship between the cell proliferation rate and the nutrient concentration:

[TABLE]

(it is also a first order Taylor expansion for fully nonlinear model), where $\widetilde{\sigma}$ is a threshold concentration, $\mu$ is a parameter expressing the “intensity” of tumor expansion, and $\xi(s;x,t)$ represents the cell location at time $s$ as cells are moving with the velocity field $\vec{V}$. The function $\xi(s;x,t)$ satisfies the ODE

[TABLE]

In another word, $\xi$ tracks the path of the cell currently located at $x$. In this problem, (1.3) describes how cells, including both interior cells and cells on the boundary, move due to the presence of time delay. Adding time delay to the basic tumor growth model makes our problem more reasonable and yet more challenging.

Furthermore, if the tumor is assumed to be of porous medium type where Darcy’s law (i.e., $\vec{V}=-\nabla p$, where $p$ is the pressure, here we consider extracellular matrix as “porous medium” in which cell moves) can be used, then

[TABLE]

Assuming the velocity field is continuous up to the boundary, we obtain the normal velocity of the moving boundary, namely,

[TABLE]

In addition, assume $\sigma$ and $p$ satisfy the boundary conditions:

[TABLE]

where $\kappa$ is the mean curvature. Equation (1.6) represents nutrient supply at the boundary and equation (1.7) represents cell-to-cell adhesiveness.

Finally, it remains to prescribe initial conditions. Instead of defining the initial conditions at time 0, for this time-delay problem, we are required to supply the initial conditions on an interval $[-\tau,0]$. For simplicity we assume initial data are time independent on the interval $[-\tau,0]$:

[TABLE]

Now the problem is reduced to mainly finding two unknown functions $\sigma$ and $p$, together with the unknown tumor region $\Omega(t)$:

[TABLE]

We shall reformulate the radially symmetric case in the next section. The radially symmetric model with time delay was studied in [4, 11, 16, 43, 41], it will be justified in the next section that the models in these papers are first-order approximations of our model in radially symmetric case. In [43], Xu, Zhou, and Bai rigorously proved that the stationary solution is always stable with respect to all radially symmetric perturbations in the case $\lambda=0$. In reality, however, we cannot ensure that perturbation is strictly radially symmetric, thus it is more natural and reasonable to ask the stability under non-radial perturbations. In this paper, we shall consider the linear stability of the unique radially symmetric stationary solution $(\sigma_{*},p_{*},R_{*})$ with respect to non-radial perturbations in the quasi-steady state case $\lambda=0$; the existence of such a solution is guaranteed:

Theorem 1.1**.**

The system $($1.10$)$–$($1.17$)$ admits a unique radially symmetric classical stationary solution $(\sigma_{*},p_{*},R_{*})$.

Next assume the initial conditions are perturbed as follows:

[TABLE]

Substituting

[TABLE]

into $($1.10$)$–$($1.17$)$ and collecting the $\varepsilon$-order terms, we obtain a linearized system around the unique radial solution $(\sigma_{*},p_{*},R_{*})$, and establish the following results:

Theorem 1.2**.**

There exists a critical value $\mu_{*}>0$ such that for any $\mu<\mu_{*}$, the radially symmetric stationary solution $(\sigma_{*},p_{*},R_{*})$ is linearly stable in the sense

[TABLE]

for some constants $a_{1}$, $b_{1}$ and $\delta>0$. If $\mu>\mu_{*}$, this stationary solution is linearly unstable.

Remark 1.1**.**

The system $($1.10$)$–$($1.19$)$ is invariant under coordinate translations, that is the reason why we exclude $a_{1}\cos(\theta)+b_{1}\sin(\theta)$ in $($1.22$)$.

Remark 1.2**.**

In this paper, we consider only 2 space dimensional case; the 3 space dimensional case can be considered in a similar manner without any difficulties, but the computations will be much more involved.

Remark 1.3**.**

In contrast to the result in [43], where stability holds for all $\mu$ with respect to radially symmetric perturbations, our $\mu_{*}$ is finite, and instability stems from $n=2$ mode when $\mu>\mu_{*}$. Recall that $\mu$ represents the tumor aggressiveness, and it is a biologically reasonable result that larger tumor aggressiveness induces instability.

Here are a couple of interesting results of the impact of time delay on the size and stability of the stationary tumor:

Theorem 1.3**.**

Adding the time delay would result in a larger stationary tumor as compared to the system without delay. The biger the tumor proliferation intensity $\mu$ is, the greater impact that time delay has on the size of the stationary tumor.

Theorem 1.4**.**

Adding the time delay to the system would not alter the critical value $\mu_{*}$ for which the stability of the stationary solution changes.

The outline of this paper is as follows. In Section 2, we reformulate the radially symmetric case. In Section 3, we establish the existence and uniqueness of the radially symmetric stationary solution $(\sigma_{*},p_{*},R_{*})$. In Section 4, we introduce the linearization of the system at $(\sigma_{*},p_{*},R_{*})$ and carry out the details of our lengthy proofs of Theorems 1.2, 1.3 and 1.4.

2. Radially Symmetric Case

In radially symmetric case, the system $($1.10$)$–$($1.19$)$ becomes

[TABLE]

where $B_{R(t)}$ denotes the disk centered at 0 with radius $R(t)$. By integrating $($2.2$)$ over $B_{R(t)}$ and using $($2.8$)$, we obtain (recall that the space dimension is $2$)

[TABLE]

We shall make a substitution $r^{\prime}=\xi(t-\tau,r,t)$ in the above integration, and derive some properties of $\xi(s,r,t)$ that will be needed in the substitution. Taking another derivative with respect to $r$ on both sides of (2.3), we have

[TABLE]

from which we find that $r^{\prime}=\xi(t-\tau;r,t)$ satisfies

[TABLE]

Furthermore, the domain of integration $r\in(0,R(t))$ becomes $r^{\prime}\in(0,R(t-\tau))$ after changing variable. This is justified in three steps as follows.

From \frac{\partial\xi}{\partial r}=\exp\big{\{}\int_{t}^{s}-\frac{\partial^{2}p}{\partial r^{2}}(\xi,c)\mathrm{d}c\big{\}}>0, we find that $r^{\prime}=\xi(t-\tau;r,t)$ is a monotone increasing function of $r$.

• 2.

If $r=0$, then $r^{\prime}=\xi(t-\tau;0,t)=0$. As a matter of fact, $\xi\equiv 0$ is the unique solution to the ODE (here it is clear that $\frac{\partial p}{\partial r}(0,s)=0$)

[TABLE]

• 3.

Now we claim $\xi(t-\tau;R(t),t)=R(t-\tau)$. Indeed, both $\xi(s,R(t),t)$ and $R(s)$ satisfy the same ODE

[TABLE]

By the uniqueness of the ODE solution, we derive $\xi(s;R(t),t)=R(s)$. Letting $s=t-\tau$, we conclude $\xi(t-\tau;R(t),t)=R(t-\tau)$.

Thus we conclude that the domain of integration under the change of variable $r^{\prime}=\xi(t-\tau,r,t)$ becomes $(0,R(t-\tau))$. Integrating $($2.5$)$ over the interval $[t-\tau,t]$, we get

[TABLE]

and substitute it into $($2.11$)$ to obtain

[TABLE]

The time delay $\tau$ in our model is actually very small; it is therefore reasonable to drop the $O(\tau)$ terms from $($2.12$)$ in radially symmetric case, with this approximation our model coincides with those in [4, 11, 16, 43, 41].

3. Radially Symmetric Stationary Solution

In this section we prove the existence and uniqueness of the radially symmetric stationary solution $(\sigma_{*}(r),p_{*}(r),$ $R_{*})$ to the system $($1.10$)$–$($1.17$)$. After setting all $t$-derivative terms to be [math], a stationary solution $(\sigma_{*}(r),p_{*}(r),R_{*})$ satisfies

[TABLE]

Lemma 3.1**.**

For sufficiently small $\tau$, there exists a unique classical solution $(\sigma_{*},p_{*},R_{*})$ to the problem $($3.1$)$-$($3.7$)$.

Proof.

To begin with, we introduce a change of variables

[TABLE]

Solving the ODE $($3.5$)$ and substituting in $($3.2$)$ and $($3.7$)$, we obtain a new system in the fixed domain $\{\widehat{r}<1\}$. After dropping the $\;\widehat{}\;$ in the above variables, the new PDE system takes the following form:

[TABLE]

We first solve $($3.11$)$ explicitly as

[TABLE]

Take $R_{\min}$ and $R_{\max}$ to be determined later. For any $R_{\min}\leq R_{*}\leq R_{\max}$, it is clear that $\sigma$ is uniquely determined by $($3.14$)$. Substituting $($3.14$)$ into $($3.12$)$, we shall prove that $p$ is also uniquely determined when $R_{*}$ is bounded by using contraction mapping principle.

Note that $\xi(s;r,0)$ ($-\tau\leq s\leq 0$) might be out of the region $[0,1]$ by following the ODE $($3.10$)$ even if $0\leq\xi(0,s,0)=r<1$ (i.e., started within the unit disk); however $\xi(s;r,0)$ should be very close to $r$ if $\tau$ is small enough. It is natural to assume $\xi(s;r,0)$ locates within $[0,2]$, so we take

[TABLE]

For each $p\in\mathcal{P}$, we first solve $\xi$ from $($3.10$)$, and substitute it into $($3.12$)$, hence we shall obtain a unique solution $\overline{p}\in W^{2,\infty}[0,1]$ from the following system:

[TABLE]

It follows from integrating $($3.15$)$ that

[TABLE]

[TABLE]

where $\sigma_{\max}=\max\limits_{0\leq r\leq 2}\sigma(r;R_{*})$. Note that $\overline{p}$ derived from $($3.15$)$ is only defined for $r\leq 1$, we shall extend $\overline{p}$ to a bigger region. Define an extension of $\overline{p}$ as

[TABLE]

It is easy to check $\|\widetilde{p}\|\in W^{2,\infty}[0,2]$, and $\|\widetilde{p}\|_{W^{2,\infty}[0,2]}\leq 2\|\overline{p}\|_{W^{2,\infty}[0,1]}$. Combining with $($3.16$)$ and $($3.17$)$, we have

[TABLE]

Define the mapping $\mathcal{L}:p\rightarrow\widetilde{p}$. If we choose $M\geq M_{1}$, then by $($3.19$)$, $\widetilde{p}\in\mathcal{P}$. Thus $\mathcal{L}$ maps $\mathcal{P}$ to itself. In the next step, we shall prove that $\mathcal{L}$ is a contraction.

Let $p_{1},p_{2}\in\mathcal{P}$, we solve $\xi_{1},\xi_{2}$ from the following two systems:

[TABLE]

[TABLE]

Integrating $($3.22$)$ and $($3.25$)$ in $s$ and making a subtraction, we obtain the following estimate

[TABLE]

for all $-\tau\leq s\leq 0$ and $0\leq r\leq 1$, hence

[TABLE]

We then substitute $\xi_{1},\xi_{2}$ into $($3.15$)$ and solve for $\overline{p}_{1}$ and $\overline{p}_{2}$, respectively. From $($3.15$)$, $\overline{p}_{1}-\overline{p}_{2}$ satisfies $(\overline{p}_{1}-\overline{p}_{2})(1)=0$ and

[TABLE]

Using $($3.26$)$, we have

[TABLE]

and similarly,

[TABLE]

[TABLE]

where M_{2}=\frac{\mu}{2}\|\frac{\partial\sigma}{\partial r}\|_{L^{\infty}[0,2]}\Big{(}1+\frac{M\tau}{(R_{\min})^{3}-M\tau}\Big{)}, M_{3}=\frac{\mu}{4}\|\frac{\partial\sigma}{\partial r}\|_{L^{\infty}[0,2]}\Big{(}1+\frac{M\tau}{(R_{\min})^{3}-M\tau}\Big{)} and M_{4}=\frac{3\mu}{2}\|\frac{\partial\sigma}{\partial r}\|_{L^{\infty}[0,2]}\Big{(}1+\frac{M\tau}{(R_{\min})^{3}-M\tau}\Big{)} are independent of $r$. It is clear that $M_{4}>M_{2}>M_{3}$, thus

[TABLE]

$\overline{p}_{1}$ and $\overline{p}_{2}$ are extended in the same way as in $($3.18$)$, hence

[TABLE]

It is easy to derive $\widetilde{p}_{1}-\widetilde{p}_{2}\in W^{2,\infty}[0,2]$ and $\|\widetilde{p}_{1}-\widetilde{p}_{2}\|_{W^{2,\infty}[0,2]}\leq 2\|\overline{p}_{1}-\overline{p}_{2}\|_{W^{2,\infty}[0,1]}$. Combining with $($3.27$)$, we have

[TABLE]

We thus obtain a contraction mapping $\mathcal{L}$ by taking $\tau$ small so that $2M_{4}\tau<1$.

Now for any particular $R_{*}\in[R_{\min},R_{\max}]$, $\sigma$ and $p$ are uniquely determined, it remains to show that there exists a unique solution $R_{*}$ satisfying $($3.13$)$. Substituting $($3.14$)$ into $($3.13$)$, we find that it is equivalent to show that there exists a unique solution to the equation:

[TABLE]

In order to prove the above statement, we set

[TABLE]

then F(R,0)=\int_{0}^{1}\Big{[}\frac{I_{0}(rR)}{I_{0}(R)}-\widetilde{\sigma}\Big{]}r\mathrm{d}r=P_{0}(R)-\frac{\widetilde{\sigma}}{2}, where $P_{0}(R)=\frac{I_{1}(R)}{RI_{0}(R)}$. From [12] (pg.61) and [22] (the equations (2.19) (2.21) and (2.26)), it is known that $P_{0}(R)$ is decreasing in $R$ and $0<P_{0}(R)\leq\frac{1}{2}$. Since $0<\widetilde{\sigma}<1$, there exists a unique solution to the equation $F(R,0)=0$, denoted by $R_{S}$. In addition, from the fact that $F(R,0)$ is monotone decreasing in $R$, we have

[TABLE]

Next we take derivative of $F(R,\tau)$ with respect to $R$ and expand the partial derivative in $\tau$ to get

[TABLE]

When $\tau$ is small enough, the signs of $\frac{\partial F(R,\tau)}{\partial R}$ and $\frac{\partial F(R,0)}{\partial R}$ should be the same, thus $F(R,\tau)$ is also monotone decreasing in $R$. In addition, from the fact that $F(R,\tau)$ is continuous in $\tau$, we have

[TABLE]

Hence there exists a unique solution $R_{*}$ satisfying $F(R_{*},\tau)=0$, i.e., equation $($3.13$)$, when $\tau$ is small enough; furthermore we have $R_{\min}=\frac{1}{2}R_{S}<R_{*}<\frac{3}{2}R_{S}=R_{\max}$. The proof is complete. ∎

Remark 3.1**.**

Since $\partial p_{*}/\partial r=0$ on the boundary for the radially symmetric stationary solution, we have $\xi(s;R_{*},0)\equiv R_{*}$ and $\xi(s;r,0)$ will stay within the unit disk if initially $\xi(0;r,0)=r<R_{*}$. From $($3.2$)$ it is clear that $\frac{\partial p}{\partial r}$ is not identically [math] for $0<r<R_{*}$, therefore $\xi(s;r,0)$ will not be a constant for $-\tau\leq s\leq 0$.

For our stationary solution, the free boundary does not move in time. But the velocity field inside the tumor domain is not zero, and movements are necessary to replace dead cells with new daughter cells to reach an equilibrium. Because of the time delay, such replacement requires a time $\tau$ for the mitosis to complete and for the daughter cells to move into the right place; and that is incorporated into the equation $($3.5$)$. In that sense, the delay-time derivative cannot be set to zero even for our stationary solution and our solution differs from the classical definition of stationary solution where time derivatives are all zero.

Throughout the paper, we denote the corresponding radially symmetric stationary solution by $(\sigma_{*},p_{*},R_{*})$.

4. Linear Stability

In this section, we consider the linearized problem of system $($1.10$)$–$($1.17$)$ with $\lambda=0$ and shall determine a critical value $\mu_{*}$ such that $(\sigma_{*},p_{*},R_{*})$ is linearly stable in the interval $0<\mu<\mu_{*}$ and linearly unstable for $\mu>\mu_{*}$. We shall also discuss the impact of time delay $\tau$ on the stability and the size of tumor.

We begin by making some small non-radially symmetric perturbations on the initial conditions (Note that the perturbations are made in a time interval $[-\tau,0]$ instead of an initial time due to the presence of time delay, and we assume for simplicity that the perturbation is uniform on the interval $[-\tau,0])$:

[TABLE]

To linearize $($1.10$)$–$($1.17$)$, we let

[TABLE]

Since we are considering a domain which is a small perturbation of a disk, we shall express $\xi(s;r,\theta,t)$ in equation $($1.14$)$ in polar coordinates $(\xi_{1}(s;r,\theta,t),\xi_{2}(s;r,\theta,t))$, where $\xi_{1}$ represents radius, and $\xi_{2}$ represents angle. Thus, the vector $\xi$ is expressed in the form $\xi=\xi_{1}\vec{e}_{1}(\xi)$, where $\vec{e}_{1}(\xi)=\cos(\xi_{2})\vec{i}+\sin(\xi_{2})\vec{j}$ and $\vec{e}_{2}(\xi)=-\sin(\xi_{2})\vec{i}+\cos(\xi_{2})\vec{j}$ are the two basis vectors in polar coordinates. We then expand $\xi_{1},\xi_{2}$ in $\varepsilon$ as

[TABLE]

Accordingly, $-\nabla$ is also expressed in polar coordinates, i.e., $-\nabla=-\vec{e}_{1}\frac{\partial}{\partial r}-\frac{1}{r}\vec{e}_{2}\frac{\partial}{\partial\theta}.$ Since $\frac{\mathrm{d}\vec{e}_{1}(\xi)}{\mathrm{d}s}=\vec{e}_{2}(\xi)\frac{\mathrm{d}\xi_{2}}{\mathrm{d}s},$ equation $($1.14$)$ is equivalent to

[TABLE]

from which we obtain two sets of ODEs in polar coordinates:

[TABLE]

[TABLE]

Substituting $($4.3$)$ and $($4.4$)$ into the above ODEs, and dropping higher order terms, we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that the equation for $\xi_{10}$ is the same as the equation for $\xi$ in radially symmetric case (i.e., $($4.5$)$ and $($2.5$)$ are the same), thus $\xi_{10}$ is independent of $\theta$; and from $($4.7$)$ we can easily derive $\xi_{20}\equiv\theta$.

Substituting $($4$)$ and $($4.5$)$–$($4.8$)$ into $($1.10$)$–$($1.17$)$, using also the mean-curvature formula in the 2-dimensional case for the curve $r=\rho$:

[TABLE]

and collecting only the linear terms in $\varepsilon$, we obtain the linearized system in $B_{R_{*}}$ ($B_{R_{*}}$ denotes the disk centered at 0 with radius $R_{*}$), namely,

[TABLE]

where the equations for $\xi_{10}$ and $\xi_{11}$ are listed in $($4.5$)$ and $($4.6$)$, respectively. Since $\xi_{21}$ does not appear explicitly in $($4.9$)$–$($4.13$)$, it is not needed.

In what follows, we seek solutions of the form

[TABLE]

Noting that in a similar manner, we can also seek solutions of the form

[TABLE]

Using the relation $\Delta=\partial_{rr}+\frac{1}{r}\partial_{r}+\frac{1}{r^{2}}\partial_{\theta\theta}$ in $($4.9$)$–$($4.13$)$, we obtain the following system in $B_{R_{*}}$:

[TABLE]

[TABLE]

[TABLE]

where the steady state solution $(\sigma_{*},p_{*},R_{*})$ satisfies $($3.1$)$–$($3.7$)$, $\xi_{10}$ satisfies $($4.5$)$ and $\varphi_{n}$ satisfies the following equation:

[TABLE]

4.0. Properties of Bessel Functions

In the sequel, we shall use modified Bessel functions $I_{n}(\xi)$ for $n\geq 0$. For convenience, we collect some properties of these functions here.

Recall that the modified Bessel function $I_{n}(\xi)$ satisfies the differential equations

[TABLE]

and is given by

[TABLE]

By [14, 22, 27], $I_{n}(\xi)$ satisfies

[TABLE]

These properties of $I_{n}(\xi)$ are needed in the subsequent discussions.

4.1. Expansion in $\tau$.

It is impossible to solve the system $($3.1$)$–$($3.7$)$, $($4.14$)$–$($4.19$)$ explicitly. However, we would like to study the impact of $\tau$ on this system. Since the time delay $\tau$ is actually very small, we look for the expansion in $\tau$ for the system $($3.1$)$–$($3.7$)$, $($4.14$)$–$($4.19$)$. Let us denote

[TABLE]

In order to compute the expansion in $\tau$ for $($4.14$)$–$($4.19$)$, we need to compute $\frac{\partial\sigma_{*}}{\partial r}$ in $($4.15$)$ $($4.16$)$, and $\frac{\partial^{2}p_{*}}{\partial r^{2}}$ in $($4.18$)$ $($4.19$)$. To do that, we expand system $($3.1$)$ — $($3.7$)$ in $\tau$. It follows from $($3.1$)$ and $($3.6$)$ that

[TABLE]

and therefore,

[TABLE]

To find $\frac{\partial^{2}p_{*}}{\partial r^{2}}$, we start with $($3.2$)$ and $($3.5$)$. We first integrate equation $($3.5$)$ over the interval $(-\tau,0)$ to obtain

[TABLE]

i.e.,

[TABLE]

We then substitute the above expression for $\xi(-\tau;r,0)$ into $($3.2$)$, since

[TABLE]

we derive the equations for $p_{*}^{0}$ and $p_{*}^{1}$,

[TABLE]

The boundary condition $p_{*}(R_{*})=\frac{1}{R_{*}}$ is expanded as follows:

[TABLE]

Thus, we have

[TABLE]

Next we expand the integral equation $($3.7$)$ using $($4.31$)$:

[TABLE]

By $($4.29$)$, the first part of $($4.36$)$ is integrated explicitly as

[TABLE]

Combining $($4.36$)$ and $($4.37$)$, we derive

[TABLE]

which leads to a set of two equations,

[TABLE]

These two equations determine $R_{*}^{0}$ and $R_{*}^{1}$, respectively.

Similarly, $w_{n}^{0}$ and $w_{n}^{1}$ satisfy the same equation $($4.14$)$. Expanding $($4.15$)$ we find

[TABLE]

which gives

[TABLE]

The next step is to expand $($4.16$)$ and $($4.19$)$ in $\tau$. Noting that

[TABLE]

and using $($4.5$)$,

[TABLE]

we deduce,

[TABLE]

Applying $($4.42$)$–$($4.44$)$ into $($4.16$)$, we derive equations for $q_{n}^{0}$ and $q_{n}^{1}$, respectively,

[TABLE]

To get the boundary condition for $q_{n}^{0}$ and $q_{n}^{1}$, we write $($4.17$)$ as

[TABLE]

hence

[TABLE]

Finally from $($4.18$)$ we have

[TABLE]

which implies

[TABLE]

We now group all the zeroth-order terms and the first-order terms in $\tau$, respectively, leading to two separate systems.

4.2. zeroth-order terms in $\tau$

Collecting the zeroth-order terms from $($4.29$)$, $($4.32$)$, $($4.34$)$, $($4.38$)$, $($4.40$)$, $($4.45$)$, $($4.47$)$, and $($4.49$)$, we obtain the following system in $B_{R_{*}^{0}}$,

[TABLE]

We first solve $p_{*}^{0}(r)$ and $w_{n}^{0}(r,t)$ from $($4.53$)$ and $($4.54$)$ as

[TABLE]

from which we compute the following terms needed in the subsequent computation,

[TABLE]

in deriving $($4.60$)$ we also made use of $($4.52$)$. To find $q_{n}^{0}$, let $\eta_{n}^{0}=q_{n}^{0}+\mu w_{n}^{0}$. Combining $($4.54$)$ and $($4.55$)$, we find that $\eta_{n}^{0}$ satisfies

[TABLE]

and its solution is given in the form

[TABLE]

thus,

[TABLE]

where $C_{1}(t)$ is determined by the boundary condition $($4.55$)$. Using also $($4.58$)$, we get

[TABLE]

In order to calculate $\frac{\partial q_{n}^{0}}{\partial t}(R_{*}^{0},t)$ in $($4.56$)$, we use $($4.23$)$, $($4.58$)$, $($4.62$)$, and $($4.63$)$ to obtain

[TABLE]

Taking another derivative with respect to $r$, we have

[TABLE]

which will be needed in the subsequent calculations. Now substituting $($4.60$)$ and $($4.64$)$ into $($4.56$)$, we derive

[TABLE]

which integrates to

[TABLE]

We shall discuss the long-time behavior of $\rho_{n}^{0}(t)$ based on $($4.66$)$. As will be seen, the analysis is different for $n=0$, $n=1$, and $n\geq 2$.

Lemma 4.1**.**

For $n=0$ and any $\mu>0$, there exists $\delta>0$ such that $|\rho_{0}^{0}(t)|\leq|\rho_{0}^{0}(0)|e^{-\delta t}$, for all $t>0$.

Proof.

When $n=0$, $\frac{n(n^{2}-1)}{(R_{*}^{0})^{3}}=0$, $($4.66$)$ becomes

[TABLE]

It suffices to show

[TABLE]

This inequality is equivalent to (3.22) in [31], which has been established already. ∎

Remark 4.1**.**

$n=0$* represents radially-symmetric perturbations. Indeed, in this case*

[TABLE]

When $\tau$ is small, we do not expect the first-order to have a major contribution, and the above result is just another indication that the stability discussed in [43] is valid for all $\mu$.

Lemma 4.2**.**

For $n=1$ and any $\mu>0$, we have $\rho_{1}^{0}(t)=\rho_{1}^{0}(0)$, for all $t>0$.

Proof.

When $n=1$, $\frac{n(n^{2}-1)}{(R_{*}^{0})^{3}}=0$, by $($4.66$)$, we have

[TABLE]

since $I_{0}(x)-I_{2}(x)=\frac{2}{x}I_{1}(x)$ by $($4.25$)$. ∎

Lemma 4.3**.**

For $n\geq 2$,

[TABLE]

The proof of this lemma can be found in [31] (Lemma 3.3). For $n\geq 2$, we define $\mu_{n}^{0}$ to be the solution to

[TABLE]

that is,

[TABLE]

Lemma 4.3 implies that $\mu^{0}_{n}>0$. We then have the following lemma.

Lemma 4.4**.**

For $n\geq 2$, $\mu_{n}^{0}<\mu_{n+1}^{0}$.

Proof.

By $($4.69$)$, we only need to establish the inequality

[TABLE]

Using Lemma 4.3, it suffices to show

[TABLE]

The above inequality has been established in [31]. The proof is complete. ∎

Since Lemmas 4.1 and 4.2 are valid for all $\mu$, we define $\mu_{0}^{0}=\mu^{0}_{1}=\infty$. And set

[TABLE]

Then by Lemma 4.4,

[TABLE]

Combining Lemmas 4.3, 4.4 and equation $($4.72$)$, it is easy to derive the following result.

Lemma 4.5**.**

For $n\geq 2$ and $\mu<\mu_{*}$, there exists $\delta>0$ such that

[TABLE]

where $\delta$ is independent of $n$.

Proof.

Since $\mu<\mu_{*}$, there exists $\delta_{1}>0$ independent of $n$ such that

[TABLE]

is valid for $n$ sufficiently large, i.e., $n>n_{0}$. On the other hand, for each $n\in[2,n_{0}]$, there exists a corresponding $\delta_{n}>0$ such that

[TABLE]

Choosing $0<\delta<\min\{\delta_{1},\delta_{2},\cdots,\delta_{n_{0}}\}$, we thus have

[TABLE]

holds for all $n\geq 2$. Therefore, it follows from $($4.66$)$ that

[TABLE]

Remark 4.2**.**

Lemma 4.5 indicates that when the time delay $\tau$ is small enough, and tumor proliferation intensity $\mu$ is smaller than a critical value (i.e., $\mu<\mu_{*}$), then the stationary solution $(\sigma_{*},p_{*},R_{*})$ is linearly stable even under non-radially symmetric perturbations. However, in contrast to the result in [43], we showed that the system is unstable with respect to perturbation when $\mu>\mu_{*}$. As indicated earlier, the instability comes from $n=2$ mode, which does not contradict the result in [43].

4.3. Sign of $R_{*}^{1}$

In 4.1, we have derived the equation for $R_{*}^{1}$. Since $R_{*}=R_{*}^{0}+\tau R_{*}^{1}+O(\tau^{2})$, we would like to know how the time delay $\tau$ affects the size of the tumor $R_{*}$, thus we are interested in the sign of $R_{*}^{1}$.

Theorem 4.1**.**

$R_{*}^{1}>0$, and $R_{*}^{1}$ is monotone increasing in $\mu$.

Proof.

Substituting $($4.29$)$, $($4.30$)$, and $($4.57$)$ into $($4.39$)$, recalling also the equality $\frac{\widetilde{\sigma}}{2}=\frac{I_{1}(R_{*}^{0})}{R_{*}^{0}I_{0}(R_{*}^{0})}$ from $($4.52$)$, we obtain,

[TABLE]

It follows from $($4.24$)$ that

[TABLE]

and by applying $($4.22$)$ and $($4.25$)$ we have

[TABLE]

Using the above equations, we shall write the integral in $($4.74$)$ explicitly as

[TABLE]

It then follows from $($4.74$)$ and $($4.75$)$ that

[TABLE]

Let

[TABLE]

[TABLE]

To determine the sign of $R_{*}^{1}$, we need to determine the signs of $A$ and $B$.

By $($4.26$)$ and $($4.27$)$, it follows

[TABLE]

[TABLE]

hence

[TABLE]

[TABLE]

From $($4.76$)$

[TABLE]

so we can directly derive $R_{*}^{1}>0$ by using $($4.77$)$ and $($4.78$)$. Furthermore, it is easy to tell that $R_{*}^{1}$ is monotone increasing in $\mu$.∎

Remark 4.3**.**

Since $R_{*}^{1}>0$, adding the time delay to the system would result in a larger stationary tumor. It is pretty reasonable because compared with models without time delay, there is more time for the tumor to grow in models with time delay. Theorem 4.1 also indicates that the biger the tumor proliferation intensity $\mu$ is, the greater impact that time delay has on the size of the tumor.

4.4. first-order terms in $\tau$

In what follows, we are going to tackle the system involving all the first-order terms in $\tau$, and we are more interested in the impact of time delay $\tau$ on our system. We now collect first-order equations and their respective boundary conditions from $($4.30$)$, $($4.33$)$, $($4.35$)$, $($4.39$)$, $($4.41$)$, $($4.46$)$, $($4.48$)$ and $($4.50$)$:

[TABLE]

[TABLE]

In $($4.89$)$, $p_{*}^{0}$ and $q_{n}^{0}$ are already computed, we only need to compute $\frac{\partial^{2}p_{*}^{1}}{\partial r^{2}}(R_{*}^{0})$ and $\frac{\partial q_{n}^{1}}{\partial r}(R_{*}^{0},t)$. Integrating $($4.82$)$ over $(0,r)$ with $r\mathrm{d}r$, we obtain

[TABLE]

We then substitute the expressions of $\sigma_{*}^{0}$ from $($4.51$)$, $p_{*}^{0}$ from $($4.57$)$, and $\sigma_{*}^{1}$ from $($4.80$)$ into the above equality to derive

[TABLE]

Since, by $($4.22$)$ and $($4.23$)$,

[TABLE]

[TABLE]

[TABLE]

the integral in $($4.90$)$ evaluates to

[TABLE]

Using $($4.22$)$, $($4.23$)$, and $($4.25$)$, taking another derivative and evaluating at $R_{*}^{0}$, also recalling the equality $\frac{\widetilde{\sigma}}{2}=\frac{I_{1}(R_{*}^{0})}{R_{*}^{0}I_{0}(R_{*}^{0})}$ from $($4.52$)$, we derive

[TABLE]

We completed computation of $\frac{\partial^{2}p_{*}^{1}}{\partial r^{2}}(R_{*}^{0})$. We now proceed a long and tedious journey to compute $\frac{\partial q_{n}^{1}}{\partial r}(R_{*}^{0},t)$. From $($4.85$)$, $w_{n}^{1}(r,t)$ can be solved in the form

[TABLE]

Substituting it into the boundary condition in $($4.85$)$, using also $($4.22$)$, $($4.23$)$, $($4.51$)$, $($4.59$)$, and $($4.80$)$, we derive

[TABLE]

Thus $C_{3}(t)$ is uniquely determined, and

[TABLE]

As in the computation of $q_{n}^{0}$ and $w_{n}^{0}$, we let $\eta_{n}^{1}=q_{n}^{1}+\mu w_{n}^{1}$. Combining $($4.85$)$ and $($4.88$)$, we find that $\eta_{n}^{1}$ satisfies

[TABLE]

with the boundary condition

[TABLE]

For simplicity, let us denote the differential operator by $L_{n}:=-\partial_{rr}-\frac{1}{r}\partial_{r}+\frac{n^{2}}{r^{2}}$, and rewrite the solution $\eta_{n}^{1}$ to $($4.94$)$ and $($4.95$)$ as $\eta_{n}^{1}=u_{n}^{(1)}+u_{n}^{(2)}+u_{n}^{(3)}+u_{n}^{(4)}$, where $u_{n}^{(1)},u_{n}^{(2)},u_{n}^{(3)}$ and $u_{n}^{(4)}$ satisfy the following equations, respectively.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We start by analyzing $u_{n}^{(1)}$. Substituting $($4.51$)$ and $($4.62$)$ into $($4.96$)$, recalling also $($4.23$)$, $($4.59$)$, and $($4.63$)$, we derive the explicit form of $($4.96$)$, namely,

[TABLE]

Recalling the definition of Bessel function $I_{n}(r)$ in $($4.21$)$, we have $\lim\limits_{r\rightarrow 0}\frac{I_{1}(r)}{r}=\frac{1}{2}$ and $\lim\limits_{r\rightarrow 0}\frac{I_{n}(r)}{r}=0$ for $n\geq 2$, thus the right hand side of $($4.100$)$ is less than $Q(n)\rho_{n}^{0}(t)$ when $0\leq r<R_{*}^{0}$. Here $Q(n)$ is a polynomial function of $n$.

Since $\rho_{n}^{0}(t)$ has different behaviors under $n\geq 2$, $n=0$ and $n=1$, we divide the following procedures into three cases: (i) $n\geq 2$; (ii) $n=0$; and (iii) $n=1$.

**Case 1: When $n\geq 2$

**

For $n\geq 2$, we introduce the following lemma to estimate $u_{n}^{(1)}$.

Lemma 4.6**.**

Consider the elliptic problem

[TABLE]

where $n\geq 2$. If $b(\cdot,t)\in L^{2}(B_{R})$, then this problem admits a unique solution $w$ in $H^{2}(B_{R})$ with estimates

[TABLE]

[TABLE]

where the constant $C$ in $($4.103$)$ and $($4.104$)$ is independent of $n$.

Proof.

Let us consider the approximate equation to $($4.101$)$ in $\varepsilon<r<R$ with zero boundary values on $x=R$ and $x=\varepsilon$, where $\varepsilon>0$ is arbitrarily small. We denote by $w_{\varepsilon}$ the corresponding classical solution. Multiplying $($4.101$)$ by $\frac{w_{\varepsilon}}{r^{2}}$ and integrating over $B_{R}\backslash B_{\varepsilon}$, we obtain:

[TABLE]

from which it follows that

[TABLE]

The equation $($4.101$)$, together with the fact that $b(\cdot,t)\in L^{2}(B_{R})$, implies $\frac{\partial^{2}w_{\varepsilon}}{\partial r^{2}}\in L^{2}(B_{R}\backslash B_{\varepsilon})$. Therefore $w_{\varepsilon}\in H^{2}(B_{R}\backslash B_{\varepsilon})$, and

[TABLE]

where $C$ is independent of $n$. Letting $\varepsilon\rightarrow 0$, we obtain a solution $w$ to $($4.101$)$ and $($4.102$)$ with estimate $($4.103$)$. The uniqueness of solution $w$ in $H^{2}(B_{R})$ follows by taking $b=0$ and using $($4.103$)$.

Next, since $H^{2}(B_{R}\backslash B_{R/2})\hookrightarrow C^{1+1/2}(B_{R}\backslash B_{R/2})$, we have

[TABLE]

which immediately implies $($4.104$)$. ∎

Applying Lemma 4.6 on $($4.100$)$, we obtain when $n\geq 2$,

[TABLE]

[TABLE]

From Lemma 4.5, we know that when $n\geq 2$ and $\mu<\mu_{*}$, there exists a constant $\delta>0$ such that $|\rho_{n}^{0}(t)|\leq|\rho_{n}^{0}(0)|e^{-\delta n^{3}t}\text{ for all }t>0.$ It follows that

[TABLE]

[TABLE]

Similarly, we can derive the same estimates for $u_{n}^{(2)}$ and $u_{n}^{(3)}$, namely,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In order to take summation with respect to $n$ in the perturbed Fourier series, we restrict $($4.107$)$ — $($4.112$)$ to $t\geq t_{0}$ for some small positive constant $t_{0}$. With these estimates in hand, the stability for the case $\mu<\mu_{*}$ will be determined by $u_{n}^{(4)}$. The solution $u_{n}^{(4)}$ to $($4.99$)$ is clearly given in the form,

[TABLE]

where $C_{4}(t)$ is determined by the boundary condition in $($4.99$)$. Combining the boundary condition from $($4.88$)$ and $($4.93$)$, we have

[TABLE]

where $H$ and $\widetilde{H}$ are functions of $R_{*}^{0}$ and $R_{*}^{1}$. Now let us combine $u_{n}^{(1)},u_{n}^{(2)},u_{n}^{(3)}$ and $u_{n}^{(4)}$ together,

[TABLE]

Note that we only need to evaluate $\frac{\partial q_{n}^{1}}{\partial r}(R_{*}^{0},t)$, thus

[TABLE]

So far, we have calculated all the expressions needed in $($4.89$)$. To get the equation for $\rho_{n}^{1}(t)$, we substitute $\frac{\partial^{2}p_{*}^{0}}{\partial r^{2}}$ from $($4.60$)$, $\frac{\partial^{3}p_{*}^{0}}{\partial r^{3}}$ from $($4.61$)$, $\frac{\partial^{2}p_{*}^{1}}{\partial r^{2}}$ from $($4.91$)$, $\frac{\partial^{2}q_{n}^{0}}{\partial r^{2}}$ from $($4.65$)$, and $\frac{\partial q_{n}^{1}}{\partial r}$ from $($4.114$)$, into $($4.89$)$, and get for $n\geq 2$,

[TABLE]

thus,

[TABLE]

To further analyze $($4.115$)$, we introduce the following lemma.

Lemma 4.7**.**

Suppose $f(t)$ satisfies

[TABLE]

with $d_{1}\neq d_{2}>0$ and the initial value $|f(0)|$ bounded, then we have

[TABLE]

where $d=\min\{d_{1},d_{2}\}$.

Proof.

$($4.116$)$ is equivalent to

[TABLE]

and thus

[TABLE]

Integrating $s$ from 0 to $t$, we derive for any $t>0$,

[TABLE]

[TABLE]

If $d_{1}>d_{2}>0$, $e^{-d_{1}t}<e^{-d_{2}t}$, the above equation implies $|f(t)|\leq Ce^{-d_{2}t}$; if $d_{1}<d_{2}$, the above equation implies $|f(t)|\leq Ce^{-d_{1}t}$. ∎

From Lemma 4.5, we have for $\mu<\mu_{*}$,

[TABLE]

Therefore we apply Lemma 4.7 on $($4.115$)$, and get $|\rho_{n}^{1}(t)|\leq Ce^{-\delta n^{3}t}$, i.e., $|\rho_{n}^{1}(t)|$ is exponentially decreasing.

**Case 2: When $n=0$

**For $n=0$, the estimates (4.71) and (4.72) follow from the standard $L^{2}$ and Schauder theory for elliptic equations. By Lemma 4.1, there exists $\delta>0$ such that

[TABLE]

Following similar procedures as in case 1, we can derive

[TABLE]

Again, from Lemma 4.1,

[TABLE]

Using Lemma 4.7, we have $|\rho_{0}^{1}(t)|\leq Ce^{-\delta t}$ , in other words, $|\rho_{0}^{1}(t)|$ is also exponentially decreasing.

Case 3: When $n=1$

Theorem 4.2**.**

For $n=1$ and any $\mu>0$, we have $\rho_{1}^{1}(t)=\rho_{1}^{1}(0)$, for all $t>0$.

Proof.

For $n=1$, from $($4.58$)$, $($4.62$)$, and $($4.63$)$, we have

[TABLE]

Differentiating $q_{1}^{0}(r,t)$ with respect to $r$ twice and evaluating at $r=R_{*}^{0}$, using also $($4.22$)$, we obtain

[TABLE]

We already derived the formula of $w_{1}^{1}(r,t)$ in $($4.93$)$. Using $($4.25$)$ to simplify, we have

[TABLE]

Next let us find the expression for $q_{1}^{1}(r,t)$. Substituting $($4.119$)$, $($4.120$)$, and $($4.122$)$ into $($4.88$)$, noting that $\frac{\partial w_{1}^{0}}{\partial t}=-\frac{I_{1}(r)}{I_{0}(R_{*}^{0})}\frac{\partial\rho_{1}^{0}(t)}{\partial t}=0$ by Lemma 4.2, and using also $($4.22$)$, $($4.23$)$, $($4.25$)$, $($4.51$)$, $($4.52$)$, and $($4.57$)$, we derive the equation for $q_{1}^{1}(r,t)$

[TABLE]

with boundary condition

[TABLE]

To solve ODE $($4.123$)$ and $($4.124$)$, we separate the solution into particular solutions $u_{1}^{(1)}$, $u_{1}^{(2)}$, $u_{1}^{(3)}$ and general solution $u_{1}^{(4)}$, where $u_{1}^{(1)}$, $\cdots$, $u_{1}^{(4)}$ satisfy the following equations, respectively,

[TABLE]

From the properties of Bessel function $($4.22$)$–$($4.25$)$, the functions $y_{1}(r)=\frac{r}{2}(1-2I_{2}(r))$, $y_{2}(r)=\frac{r}{4}(-I_{1}^{2}(r)+I_{0}(r)I_{2}(r))$ and $y_{3}(r)=-I_{1}(r)$ satisfy the following equations

[TABLE]

[TABLE]

[TABLE]

Thus, we obtain

[TABLE]

In addition, the general solution $u_{1}^{(4)}$ takes the form of

[TABLE]

Combining $($4.129$)$, $($4.130$)$, $($4.131$)$, and $($4.132$)$, we find that the solution to $($4.123$)$ should be

[TABLE]

where $C_{5}(t)$ is determined by the boundary condition in $($4.128$)$, namely,

[TABLE]

which simplifies to

[TABLE]

In order to compute $\rho_{1}^{1}(t)$ from $($4.89$)$, we need the derivative value $\frac{\partial q_{1}^{1}}{\partial r}(R_{*}^{0},t)$. We combine $($4.132$)$ and $($4.133$)$ to obtain

[TABLE]

Applying $($4.129$)$, $($4.130$)$, $($4.131$)$, and $($4.134$)$, we then derive

[TABLE]

Finally, it remains to substitute $\frac{\partial^{2}p_{*}^{0}}{\partial r^{2}}$ from $($4.60$)$, $\frac{\partial^{3}p_{*}^{0}}{\partial r^{3}}$ from $($4.61$)$, $\frac{\partial^{2}p_{*}^{1}}{\partial r^{2}}$ from $($4.91$)$, $\frac{\partial^{2}q_{1}^{0}}{\partial r^{2}}$ from $($4.121$)$, and $\frac{\partial q_{1}^{1}}{\partial r}$ from $($4.135$)$, into equation $($4.89$)$, and collect terms containing $\rho_{1}^{0}(t)$ and $\rho_{1}^{1}(t)$,

[TABLE]

Recalling we have formula for $R_{*}^{1}$ in $($4.76$)$,

[TABLE]

substituting it into $($4.136$)$, we deduce

[TABLE]

Hence $\rho_{1}^{1}(t)=\rho_{1}^{1}(0)$, which completes the proof. ∎

Remark 4.4**.**

We just established, for $n=1$, after ignoring $O(\tau^{2})$ terms,

[TABLE]

is a transform of the origin. Thus $n=1$ mode would not affect the stability.

Acknowledgment. We are grateful to the referees for a very careful reading of our manuscript, and for helpful suggestions.