Rigorous derivation of a Hele-Shaw type model and its non-symmetric traveling wave solution
Yu Feng Qingyou He Jian-Guo Liu Zhennan Zhou


]Rigorous derivation of a Hele-Shaw type model and its non-symmetric traveling wave solution




Yu Feng: Beijing International Center for Mathematical Research, Peking University, No. 5 Yiheyuan Road Haidian District, Beijing, P.R.China 100871
fengyu@bicmr.pku.edu.cn
Qingyou He: Sorbonne Université, CNRS, Laboratoire de Biologie Computationnelle et Quantitative, F-75005 Paris, France
qyhe.cnu.math@qq.com
Jian-Guo Liu: Department of Mathematics and Department of
  Physics, Duke University, Durham, NC 27707, USA
jian-guo.liu@duke.edu
Zhennan Zhou: Institute for Theoretical Sciences, Westlake University, Hangzhou, Zhejiang Province, 310030, China
zhouzhennan@westlake.edu.cn


[2020]35R35, 76D27, 92C10, 70K50





[
    Yu Feng, Qingyou He, Jian-Guo Liu, Zhennan Zhou
    April 29, 2024
===================================================




In this paper, we consider a Hele-Shaw model that describes tumor growth subject to nutrient supply. This model was recently studied in <cit.> via asymptotic analysis.  Our contributions are twofold: Firstly, we provide a rigorous derivation of this Hele-Shaw model by taking the incompressible limit of the porous medium reaction-diffusion equation, which solidifies the mathematical foundations of the model. Secondly, from a bifurcation theory perspective, we prove the existence of non-symmetric traveling wave solutions to the model, which reflect the intrinsic boundary instability in tumor growth dynamics.












§ INTRODUCTION


Tumor boundary instability, characterized by the formation and evolution of finger-like protrusions, has been a significant area of concern in oncology research <cit.>. This phenomenon, where tumors develop irregular boundaries, is believed to facilitate more efficient invasion of surrounding healthy tissue. Consequently, malignant tumors tend to exhibit more irregular borders compared to benign ones. In this paper, we investigate this boundary instability induced by nutrient consumption and supply by studying a Hele-Shaw type model.

To begin with, we let x∈Ω(t)⊆ℝ^n denote the region occupied by the tumor tissue at time t, p(x,t) be the pressure inside the tumor, and c(x,t) represents the nutrient concentrate that supports tumor growth. The model we consider is given by




-p    = G_0 c,       in  Ω(t),


p   =0,       on  ∂Ω(t),

-Δc+λc   =0,       in  Ω(t),

-Δc+ c   =c_B,       in  ℝ^n∖Ω(t),


where the parameters G_0, λ, c_B>0. This model was first proposed in <cit.>, and it can be interpreted as follows. Here, G_0 is the growth parameter, and G_0c is the growth rate function. Inside the tumor, nutrients are consumed by the tumor cells at a rate of λ, whereas outside the tumor, nutrients are supplied by the vascular network in the healthy region, with the supply rate proportional to the concentration difference c_B-c. Additionally, the tumor region Ω(t) evolves over time, and the boundary evolution is characterized via Darcy's law 

    v|_∂Ω=-∇ p|_∂Ω.


Recently, the authors in <cit.> studied the boundary instability of (<ref>) via an asymptotic analysis approach, which complements the current understanding of this model <cit.>. They introduced a small perturbation around the symmetric solutions with amplitude ϵ and a profile chosen from a set of basis functions indexed by the frequency l∈ℕ, which reduced the evolution of the boundary perturbation to the dynamics of the perturbation amplitude. The authors in <cit.> derived the relative rate equation of the amplitude evolution ϵ^-1dϵ/dt=E(,l) and characterized the boundary instability by determining its sign. When E(,l) takes a positive value, the finger-like structures grow; otherwise, the boundary degenerates to the symmetric one. The main result in <cit.> interprets that the nutrient consumption rate λ can trigger the boundary instability in (<ref>). Specifically, the boundary remains stable for any perturbation frequency if λ≤ 1. However, for λ>1, a threshold value exists for the perturbation frequency, below which boundary instability occurs, although higher frequencies remain stable.

In this paper, we contribute to this model in two aspects. Firstly, we rigorously derive the Hele-Shaw model (<ref>) by investigating the incompressible limit of a porous medium equation (PME) type cell density model, which poses new technical challenges. Secondly, we establish the existence of non-symmetric traveling wave solutions to (<ref>) by relating such solutions to the non-trivial branch of solutions of a nonlinear functional equation. In the following, we begin with a detailed literature review of incompressible limit studies on tumor growth models, situating our results within the current literature. 

It is well known that PME type equations possess a limit as m→∞, which is shown to be a Hele-Shaw problem. Perthame et al. first generalized related studies to tumor growth models in the seminal work <cit.>, which facilitates numerous impressive works in this direction  <cit.>. The Hele-Shaw asymptotic limit of the tumor growth model <cit.> was initially studied. 

For the tumor growth model with Brinkman's pressure law governing the motion,  authors in <cit.> established an optimal uniform decay rate of the density and the pressure in m, the Hele-Shaw (incompressible) limits of the two-species case were proved in <cit.>.  The Hele-Shaw limit of the PME with the non-monotonic (and nonlocal) reaction terms through the approach of the obstacle problem was completed in <cit.>. The existence of the weak solution and the free boundary limits of a tissue growth model with autophagy or necrotic core were obtained in <cit.> respectively.  In addition, the inconpressible limits for the chemotaxis (even with growth term) were shown in <cit.>. In addition, the singular limit of the PME with a drift was discussed in <cit.>. Recently, the convergence of free boundaries in the incompressible limit of tumor growth models was considered in <cit.>. In this paper, the incompressible limit can be derived for the nutrient model both in the parabolic and elliptic regimes, our approach is inspired by the recent methodology proposed by <cit.>.  In the absence of contact inhibition as in <cit.>, there is no natural bound for the pressure function in the system, and the source term provides a pure growth effect.  Hence, some additional analysis techniques are required. 

Regarding the boundary instability study to (<ref>), as mentioned before, the authors in <cit.> derived the amplitude evolution equation, E(,l), of a sequence of basis functions indexed by the frequency l∈ℕ, and determined the boundary instability by its sign.  In this paper, we further reveal the intrinsic boundary instability of (<ref>) by proving the existence of non-symmetric traveling wave solutions in a prescribed domain. The proof is based on relating such traveling wave solutions to the non-trivial solutions of a nonlinear functional equation F(ξ,)=0, where ξ is a function that describes the boundary profile, and  is the nutrient consumption rate parameter.
In fact, each solution of F(ξ,)=0 corresponds to a traveling wave solution to (<ref>). 
In particular, the symmetric traveling wave solutions of (<ref>) are associated with the trivial solution branch (0,). And, proving the existence of non-symmetric solutions reduces to finding the non-trivial solutions (ξ,) of F(ξ,)=0. This is achieved by investigating the Fréchet derivative of the nonlinear map  F(ξ,) and we conclude by the celebrated Crandall-Rabinowitz theorem.




We emphasize that the model (<ref>) shares some similarity with the free boundary models proposed by Greenspan  <cit.> and further developed by Friedman et al. in <cit.>, but there are essentially differences. As discussed above, the model (<ref>) is derived from the incompressible limit of PME-type equations. Thus, the pressure, as the limit of power functions of density, has to remain nonnegative and vanish on the tumor boundary (see <cit.>). However, the pressures in the models originated from <cit.> allow negative values, and their boundary values rely on its curvature <cit.>. 
In light of such differences, the boundary instability of (<ref>) can not be implied by the previous works <cit.>, and calls for an independent study.

Finally, it is worth mentioning in some recent works <cit.>, the authors studied the boundary behavior of a variant model close to (<ref>), whereas the boundary in that model is always stable in the sense of the Wasserstein distance. 


The paper is organized as follows. In Section <ref>, we provide a rigorous derivation of the Hele-Shaw model by taking the incompressible limit of a density model of the porous medium equation type. Then, in Section <ref>, we review the asymptotic analysis results to this Hele-Shaw model in <cit.>. Finally, in Section <ref>, we conclude with proving the existence of non-symmetric traveling wave solution in this model via a Crandall-Rabinowitz argument.



§ INCOMPRESSIBLE LIMIT DERIVATION


This section is devoted to the derivation of model (<ref>) by taking the incompressible limit of a cell density model in the PME type. We begin with introducing this cell density model in the following subsection.


 §.§ A cell density model of PME type

We let ρ_m(x,t) to denote the tumor cell density, start from the initial state ρ_m,0(x), and D(t) be the supporting set of ρ_m(x,t) that is 

    D(t)={x|ρ_m(x,t)>0}.

Physically, it presents the tumoral region at time t. We assume the pressure inside the tumor follows the constitutive relationship of P_m(ρ_m)=m/m-1ρ_m^m-1 and the tumor cell moving velocity is governed by the Darcy's law v=-∇ P_m. Thus, in particular, the boundary expansion is characterized by v|_∂ D=-∇ P_m|_∂ D and D(t) remains finite for any t<∞ provided D(0) is compact. On the other hand, we employ c_m(x,t) to denote the nutrient concentration and assume the cells' production rate is proportional to it. With the above assumption, the evolution of ρ_m satisfies the following porous medium equation (PME) with source term:

    ∂_tρ_m-∇·(ρ_m∇ P_m)=G_0 c_mρ_m,   t⩾ 0,    G_0>0.

Regarding the nutrient concentration c_m, it is governed by the following reaction-diffusion equation in general:

    τ∂_t c_m-Δ c_m +Ψ(ρ_m,c_m)=0,

where the parameter τ≥ 0 characterizes the nutrient change time scale, and the binary function Ψ(ρ_m,c_m) describes the overall effects of the nutrient supply outside the tumor and the nutrient consumption inside the tumor. Considering the timescale parameter τ≪ 1 (see, e.g. <cit.>), we drop it to get the elliptical form of nutrient equation,

    -Δ c_m +Ψ(ρ_m, c_m)=0.

Then, we focus on the so-called in vivo regime in <cit.>, in which the nutrients are provided by vessels of the healthy tissue surrounding the tumor. Mathematically, we assume the nutrients is consumed at a rate of >1 in the tumor cell saturated region S(t)={x|ρ_m(x,t)≥ 1}. While, outside S(t), the nutrient supply is determined by the concentration difference from the background, c_B - c_m, with the rate of (1-ρ_m)_+. Thus, the binary function takes the form of

    Ψ(ρ_m,c_m)=λρ_m c_m·χ_S-(1-ρ_m)_+(c_B-c_m)·χ_S^c,

where f_+=max{f,0}. Or equivalently,



-Δc_m + λρ_m c_m= 0,       for  x∈S(t),

-Δc_m =(1-ρ_m)_+(c_B-c_m),       for  x∈ℝ^n∖S(t).



By now, we have finished introducing the cell density model ((<ref>) coupled with (<ref>), equivalently, (<ref>)). However, considering the rigorous derivation of (<ref>) is technical, we postpone the proof to Section <ref> and Section <ref> and instead establish a formal derivation in the next section to light up the way.



 §.§ A formal derivation

In this subsection we derive (<ref>) formally. Firstly, by multiplying mρ_m^m-1 on the both sides of equation (<ref>) one gets the equation for pressure

    ∂_t P_m=|∇ P_m|^2+(m-1) P_m(Δ P_m+G_0 c_m).

Then, by sending m→∞ and let (ρ_∞, P_∞, c_∞) to denote the limit density, pressure, and nutrient, respectfully. One can formally obtain the so-called complementarity condition:

    P_∞(Δ P_∞+G_0 c_∞)=0.

At the same time, the limit density ρ_∞ satisfies the following equation in the distributional sense:

    ∂/∂ tρ_∞-∇·(ρ_∞∇ p_∞)=G_0 ρ_∞ c_∞,

with P_∞ compels ρ_∞ only take value in the range of [0,1] for any initial date ρ_∞,0∈[0,1]. And P_∞ belongs to the Hele-Shaw monotone graph:

    P_∞(ρ_∞)={[ 0,      0⩽ρ_∞<1,; [0,∞),    ρ_∞=1.;                  ].

Observe that, in general, the supporting set of ρ_∞ is larger than that of P_∞. However, a transparent regime called "patch solutions" exists for a large class of initial data, in which the two sets coincide with each other, denoted as Ω(t). At the same time, the limit density evolves in the form of ρ_∞=χ_Ω(t), where χ_A presents the characteristic function of set A. In this specific regime, (<ref>) reduce to



-Δc_∞+ λc_∞= 0,       for  x∈Ω(t),

-Δc_∞=c_B-c_∞,       for  x∈ℝ^n∖Ω(t).


While, (<ref>) and (<ref>) together yields



-p_∞   = G_0 c_∞,       in  Ω(t),

p_∞   =0,       on  ∂Ω(t).


Thus, by dropping the subscripts and replacing P by p in (<ref>) and (<ref>) one gets (<ref>). We emphasize that before taking the incompressible limit, the pressure function (<ref>) and the nutrient functions (<ref>) are coupled strongly to each other. However, fortunately, by taking the incompressible limit and further restricting ourselves to the 
patch solution regime, the two equations decoupled automatically.



 §.§ Regularity and stiff limit


In this subsection, we turn to the rigorous derivation of (<ref>). Recall the density model writes,

    ∂_tρ_m=Δρ_m^m+G_0ρ_m c_m,   
    τ∂_tc_m=Δ c_m-Ψ(ρ_m,c_m),     
    
    c_m→ c_B,  as   |x|→∞.

where Ψ(ρ_m,c_m) is given by (<ref>), and the pressure 
P_m:=m/m-1ρ_m^m-1 satisfies 

    ∂_t P_m=(m-1)P_m(Δ P_m+G_0c_m)+|∇ P_m|^2.

For concision, in the later section we let χ:=χ_S and χ:=χ_S^c. In order to being corresponding to the following Hele-Shaw problem, we consider τ=0 in the rest of this section.
We first give some regularity estimates, then we obtain the stiff limit through the weak solution.

Assume that m>2 and the initial data (ρ_m,0,P_m,0) satisfy

    ρ_m,0_L^1(ℝ^n)≤ C,      Δρ_m,0^m+G_0ρ_m,0 c_m,0_L^1(ℝ^n)≤ C,         ∇ρ_m,0_L^1(ℝ^n)≤ C,
        P_m,0(x)≤ C_0R_0^2h(x/R_0),       ∇ P_m,0_L^2(ℝ^n)≤ C
for some R_0,C_0>0 and h(x):=1/2n(1-|x|^2)_+. Then, it holds for 0≤ t≤ T that

    0≤ c_m≤ c_B, (x,t)∈ Q_T,   supp(P_m(t))⊂ B_R_T(0), 0≤ t≤ T,


    P_m_L^∞(Q_T)+ρ_m_L^∞(Q_T)≤ C(T),


    sup_0≤ t≤ T[c_m-c_B_L^2(ℝ^n)+∇ c_m_L^2(ℝ^n)+∂_t c_m_L^1(ℝ^n)]≤ C(T),


    sup_0≤ t≤ T[∂_tρ_m_L^1(ℝ^n)+∇ρ_m_L^1(ℝ^n)]≤ C(T),


    ∂_t P_m_L^1(Q_T)+∇ P_m_L^2(Q_T)≤ C(T),

where R_T=R_0e^C_0T/n.



By the comparison principle,  we directly obtain 

    0≤ c_m≤ c_B,   (x,t)∈ Q_T.

 Inspired by <cit.>, let C_0≥ G_0c_B and

    h(x):=1/2n(1-|x|^2)_+.

Then, for given ϕ:=C_0R^2(t)h(x/R(t)), R(t):=R_0e^C_0t/n, it holds on the support of ϕ that 

    ∂_tϕ-(m-1)ϕ(Δϕ+G_0c_B)_≤ 0-|∇ϕ|^2≥C_0R(t)/n(R'(t)-C_0R(t)/n)=0.

Suppose that P_m,0(x)≤ C_0R_0^2ϕ(x/R_0) with some R_0>0, it follows from the comparison principle that

    P_m(x,t)≤ϕ(x,t),  (x,t)∈ Q_T.

The above conclusion means

    supp(P_m)⊂ B_R_T(0),   P_m_L^∞(Q_T)≤ C(T),

and further concludes for m>2 that 

    supp(ρ_m)⊂ B_R_T(0),   ρ_m_L^∞(Q_T)≤ C(T)
with R_T=R_0e^C_0T/n.
 We multiply (<ref>)_2 by c_m-c_B and obtain
 
    ∫_ℝ^n|∇ c_m|^2dx+∫_ℝ^n[λχρ_m+χ̃(1-ρ_m)_+]|c_m-c_B|^2dx≤∫_ℝ^nχρ_m c_B|c-c_B|dx.

Using the fact λχρ_m+χ̃(1-ρ_m)_+≥min{λχ,χ̃}>0, the estimate (<ref>), and Hölder's inequality yields 
 
    sup_0≤ t≤ Tc_m(t)-c_B_L^2(ℝ^n)+sup_0≤ t≤ T∇ c_m(t)_L^2(ℝ^n)≤ C(T).

 We use the equation (<ref>)_2 and Kato's inequality, it holds  
 
    0≤Δ |∂_t c_m|-|∂_t c_m|(λχρ_m+χ̃(1-ρ_m)_+)+|∂_tρ_m|(G_0χ c_m+χ̃|c_m-c_B|),

 which implies 
 
    ∫_ℝ^n|∂_t c_m|dx≤ C∫_ℝ^n|∂_tρ_m|dx, t>0.

 By means of Kato's inequality for (<ref>)_1, we have 
 
    ∂_t |∂_tρ_m|≤Δ |∂_tρ_m^m|+G_0c_B|∂_tρ_m|+G_0ρ_m|∂_t c_m|

 Taking (<ref>) into consideration,  and integrating the above inequality (<ref>) on [0,t] for any 0≤ t≤ T, we have
 
    sup_0≤ t≤ T∂_tρ_m_L^1(ℝ^n)+sup_0≤ t≤ T∂_t c_m_L^1(ℝ^n)≤ C(T).

 Parallel to the proof of the above estimate, it holds 
 
    sup_0≤ t≤ T∇ρ_m_L^1(ℝ^n)≤ C(T).

 We multiply the inequality equation (<ref>)_1 by φ(x)≥ 0 with -Δφ =χ_B_R_T+1 and integrate on Q_T, then it yields 
 
    ∂_t ρ_m^m_L^1(Q_T)≤   φ_L^∞(ℝ^n)[(∂_tρ_m(0) _L^1(ℝ^n)+∂_tρ_m(T) _L^1(ℝ^n))
       +G_0c_B|∂_tρ_m|+G_0ρ_m|∂_t c_m|_L^1(Q_T)]≤ C(T).

 Hence, we obtain
 
    ∂_t P_m_L^1(Q_T)≤m/2^m-2∂_tρ_m_L^1(Q_T)+2∂_tρ_m^m_L^1(Q_T)≤ C(T).

 In addition, it holds by integrating (<ref>) on Q_T that 

    ∇ P_m_L^2(Q_T)   ≤m-1/m-2G_0c_BP_m_L^1(Q_T)+1/m-1[P_m,0_L^1(ℝ^n)+P_m(T)_L^1(ℝ^n)]
       ≤ C(T).

Combining (<ref>),(<ref>),(<ref>),(<ref>),(<ref>),(<ref>),(<ref>), and (<ref>), we complete the proof.


In the following Theorem, based on the basic energy estimates in Lemma <ref>, we derive some results of convergence in m, and then prove the stiff limit.

Under the same initial assumptions of Lemma <ref>, there exist a pair of functions  (ρ_∞, P_∞, c_∞), satisfying  ρ_∞, P_∞∈ L^1(Q_T)∩ L^∞(Q_T)∩ BV(Q_T), c_∞-c_B∈ W^1,2(Q_T)∩ BV(Q_T), such that, after extracting subsequences, it holds

    ρ_m→ρ_∞,      in L^p(Q_T), 1≤ p<∞,         as m→∞,
       P_m→ P_∞,      in L^p(Q_T), 1≤ p<∞,         as m→∞,
       c_m→ c_∞,      in L^p_loc(Q_T), 1≤ p<∞,         as m→∞, 
       (1-ρ_m)_+→ (1-ρ_∞)_+,      in L^p(Q_T), 1≤ p<∞,         as m→∞.

Moreover, (ρ_∞,P_∞,c_∞) satisfies a Hele-Shaw type system as 

    ∂_tρ_∞=Δ P_∞+G_0ρ_∞ c_∞,   in 𝒟'(Q_T),
    Δ c_∞=Ψ(ρ_∞,c_∞),     in 𝒟'(Q_T),
    
    c_∞→ c_B,  as   |x|→∞.   
    
    0≤ρ_∞≤ 1,   P_∞(1-ρ_∞)=0,   a.e. in Q_T.



Based on the verified estimates (<ref>)-(<ref>), then (<ref>)-(<ref>) hold by the compactness embedding. Furthermore, (<ref>) naturally  supports (<ref>).   Taking (<ref>)-(<ref>) into account, hence (<ref>)_1-3 holds in the sense of distribution.

 Since 0≤ P_m≤ C(T), it holds by passing to limit that 
 
    0≤ρ_∞≤ 1, a.e. in Q_T.

 In addition, we have P_m^m/m-1=(m/m-1)^1/m-1ρ_m P_m, which concludes after taking the limit in m that 
  
    P_∞(1-ρ_∞)=0, a.e. in Q_T,

  and the proof is completed. 



 §.§ Complementarity relationship
 

In this subsection, we give the complementarity relation which describes the limiting pressure P_∞ through a degenerate elliptic equation. To this end, we use the viewpoint of obstacle problem introduced initially in <cit.> to study the incompressible limit. 
Suppose that the same initial assumptions as Lemma (<ref>) hold. For all t>0, let E_t denotes the space 

    E_t:={v∈ H^1(ℝ^n)∩ L^1(ℝ^n) | v≥ 0,⟨ v,1-ρ_∞(t)⟩_H^1,H^-1=0}.

Then, for all t>0, the function x→ P_∞(x,t) is the solution of the minimization problem:

    P_∞∈ E_t,   
    ∫_ℝ^n|∇ P_∞|^2/2-G_0c_∞ P_∞ dx≤∫_ℝ^n|∇ v|^2/2-G_0 c_∞ v dx,   for all v∈ E_t.

Furthermore, P_∞ satisfies the following complementarity condition:

    P_∞(Δ P_∞+G_0c_∞)=0  in 𝒟'(ℝ^n× (0,∞)).




To begin with, we recall the equation for the pressure

    ∂_tP_m=(m-1)P_m(Δ P_m+G_0c_m)+|∇ P_m|^2.

Given t_0>0 and  a function v in E_t_0, we use the equation for the pressure (<ref>) and the density (<ref>)  to write:

    ∫_ℝ^n   ∇ P_m·∇ P_m-ρ_m∇ P_m·  v+G_0c_m(v-P_m)dx
    
    =   -1/m-1[d/dt∫_ℝ^nP_mdx-∫_ℝ^n|∇ P_m|^2dx]+d/dt∫_ℝ^nvρ_mdx
       +∫_ℝ^nG_0c_mv(1-ρ_m)dx.

Integrating on (t_0,t_0+δ) for any δ>0, we have

    ∫_t_0^t_0+δ   -4pt∫_ℝ^n∇ P_m·∇ P_m-ρ_m∇ P_m·  v+G_0c_m(v-P_m)dxdt
    
    =   -1/m-1[∫_ℝ^n(P_m(t_0+δ)-P_m(t_0))dx-∫_t_0^t_0+δ-4pt∫_ℝ^n|∇ P_m|^2dxdt]
       +∫_ℝ^nv(ρ_m(t_0+δ)-ρ_m(t_0))dx
    +∫_t_0^t_0+δ-4pt∫_ℝ^nG_0c_mv(1-ρ_m)dxdt.

Thanks to the  lower semi-continuity of L^2-norm, it holds 

    ∫_t_0^t_0+δ-4pt∫_ℝ^n∇ P_∞·∇ P_∞-∇ P_∞·  v+G_0c_∞(v-P_∞)dxdt
    ≤   lim inf_m→∞∫_t_0^t_0+δ-4pt∫_ℝ^n∇ P_m·∇ P_m-ρ_m∇ P_m·  v+G_0c_m(v-P_m)dxdt
    
    =   lim inf_m→∞∫_ℝ^nv(ρ_m(t_0+δ)-ρ_m(t_0))dx+∫_t_0^t_0+δ-4pt∫_ℝ^nG_0c_∞ v(1-ρ_∞)dxdt
    ≤   lim inf_m→∞∫_ℝ^nv(ρ_m(t_0+δ)-ρ_m(t_0))dx+G_0c_B∫_t_0^t_0+δ-4pt∫_ℝ^nv(1-ρ_∞)dxdt.

We first note that 

    d/dt∫_ℝ^nρ_mvdx=-∫_ℝ^nρ_m∇ P_m·∇ vdx+ ∫_ℝ^nG_0c_mρ_mvdx.

Then, it follows

    |d/dt∫_ℝ^nρ_mvdx|≤ C∫_ℝ^n|∇ P_m| |∇ v|dx+G_0c_B ∫_ℝ^nρ_mvdx.

The first term in the right hand side is bounded in L^2(0,T), and the second term is also bounded in L^∞(0,T). We deduce that the function t→∫_ℝ^nρ_mvdx is bounded in H^1(0,T)⊂ C^1/2(0,T) and thus converges uniformly in [0,T]. Since ∫_ℝ^nvρ_mdx converges to ∫_ℝ^nvρ_∞ dx in 𝒟'(ℝ_+), we have

    ∫_ℝ^nv(·)ρ_m(·,t)dx→∫_ℝ^nv(·)ρ_∞(·,t)dx locally uniformly in ℝ_+.

Consequently,

    lim inf_m→∞   ∫_ℝ^nv(·)(ρ_m(x,t_0+δ)-ρ_m(x,t_0))dx
    
    =   ∫_ℝ^nv(·)(ρ_∞(x,t_0+δ)-ρ_∞(x,t_0))dx
    
    =   ∫_ℝ^nv(·)(ρ_∞(x,t_0+δ)-1)dx≤ 0.

We insert (<ref>) into (<ref>), it holds

    1/δ∫_t_0^t_0+δ-4pt∫_ℝ^n|∇ P_∞|^2/2-G_0c_∞ P_∞ dxdt≤1/δ∫_t_0^t_0+δ-4pt∫_ℝ^n|∇ v|^2/2-G_0c_∞ v dxdt
       +G_0c_B/δ∫_t_0^t_0+δ⟨ v,1-ρ_∞⟩_H^1,H^-1dt.

Due to the fact P_∞∈ BV(Q_T), the trace theorem supports   P_∞^+=P_∞ for the trace P_∞^+ of P_∞. Hence, let δ→ 0^+, we get 

    ∫_ℝ^n|∇ P_∞|^2/2-G_0c_∞ P_∞ dx ≤∫_ℝ^n|∇ v|^2/2-G_0c_∞ v dx,

where G_0c_B/δ∫_t_0^t_0+δ⟨ v,1-ρ_∞⟩_H^1,H^-1dt→⟨ v,1-ρ_∞(t_0)⟩ _H^1,H^-1=0 as t→ 0^+ is used since 1-ρ_∞∈ C([0,T],H^-1).

Finally,  given a test function φ∈ C_0^∞(ℝ^n× (0,∞)), we take v_ϵ=P_∞+ϵ P_∞φ=P_∞(1+ϵφ) with |ϵ|≪ 1 so that 1+ϵφ≥0.  Due to (<ref>), we have

    d/dϵ|_ϵ=0∫_ℝ^n|∇ v_ϵ|/2-G_0c_∞ v_ϵ dx=0,

which yields 

    ∫_ℝ^n∇ P_∞(t)·∇ (P_∞(t)φ)-G_0c_∞ P_∞(t)φ dx=0,  a.e. t∈ℝ_+.

Hence, the complementarity relationship
(<ref>) holds in the sense of distribution.





§ THE ASYMPTOTIC ANALYSIS


Starting from this section, we turn to the study of the boundary instability to (<ref>) on ℝ^2, and thus we further assume x=(x,y)∈Ω(t)⊆ℝ^2. We also employ  (t) to denote the free boundary ∂Ω(t). Recall that the pressure satisfies,



    -p    = G_0 c,       in  Ω(t),

    p    = 0,        on  (t),

    p   <∞,        in  Ω(t),


with the growth rate parameter G_0>0. Then the normal speed σ at a boundary point x∈ is governed by Darcy's law: 

    σ(x)=-∇ p|_B·n̂(x),

where n̂(x) stands for the normal vector at x. While, the nutrient concentration c(x,t)∈ C^1 satisfies:



    -c + c    = 0,        in  Ω(t),

    -c + c    = c_B,        in  ℝ^2∖Ω(t),

    c   <∞,        in  ℝ^2.











The remaining part of this paper aims to discuss the existence of a non-symmetric traveling wave solution in the above model ((<ref>) coupled with (<ref>)).  For the reader's convenience, we begin with reviewing the results in the recent paper <cit.>, in which the authors studied the evolution of small perturbations of this model by using asymptotic analysis. Their calculations are closely related to the bifurcation analysis in Section <ref>. 



 §.§ The symmetric solutions


As the first step in <cit.>, the authors solved the traveling wave solution analytically under symmetry assumption. In which the boundary is given by the vertical line

    _0={(x,y)| x= 0},

with a constant traveling speed to the right, and the corresponding tumor region is the left halfplane,

    Ω_0={(x,y)| x≤ 0}.

We use superscripts (i) and (o) to distinguish the  solutions inside or outside the tumor region. By solving equation (<ref>) and (<ref>) with Ω(t) replaced by Ω_0, one gets the symmetric solutions c_0^(i)(x), c_0^(o)(x), and p_0^(i)(x) as follows:



   c_0^(i)(x)=c_B/√()+1e^√()x,     for  x≤0,

   c_0^(o)(x)=c_B-c_B·√()/√()+1e^-x,     for  x≥0,

   p_0^(i)(x)=G_0·c_B/(√()+1)-G_0·c_B/(√()+1)e^√()x,     for  x≤0.



And the boundary normal speed σ_0 is given by

    σ_0(λ)=-∂ p_0^(i)/∂ x(0)=G_0· c_B/√()(√()+1).

















 §.§ Linearization of single mode perturbation


In the second step, the authors in <cit.> added a perturbation with frequency l∈ℝ^+ and amplitude 0<ϵ≪ 1. Then the perturbed tumor region Ω and boundary  becomes


    Ω={(x,y)|x≤ϵcosly, y∈ℝ};

    ={(x,y)|x=ϵcosly, y∈ℝ}.



Corresponding to the above perturbation, the perturbed solutions have the following asymptotic expansion form 



c^(i)(x,y,t)   =c^(i)_0(x)+ϵ(t) c^(i)_1(x,y)+ϵ(t)^2 c^(i)_2(x,y)+…

c^(o)(x,y,t)   =c_0^(o)(x)+ϵ(t) c_1^(o)(x,y)+ϵ(t)^2 c_2^(o)(x,y)+…

p^(i)(x,y,t)   =p_0^(i)(x)+ϵ(t) p_1^(i)(x,y)+ϵ(t)^2 p_2^(i)(x,y)+…

with the leading order terms given by the symmetric solutions, since the perturbation amplitude is small compared to the tumor size. Moreover, the first-order terms capture the main response to the perturbation. Thus, they are variable-separable in the sense of

    c_1^(i)(x,y)=c_1,l^(i)(x)cosly,   c_1^(o)(x,y)=c_1,l^(o)(x)cosly,  
    p_1^(i)(x,y)=p_1,l^(i)(x)cosly.

Then by plugging the expansion (<ref>) into (<ref>) and (<ref>), one can first obtain the equations for (c_1,l^(i), c_1,l^(o), p_1,l^(i)),


-∂_x^2 c_1,l^(i)+(λ+l^2)c_1,l^(i)   =0,     for  x≤0,

-∂_x^2 c_1,l^(o)+(1+l^2)c_1,l^(o)   =0,     for  x≥0,

-∂_x^2 p_1,l^(i)    = G_0  c_1,l^(i),     for  x≤0,


Secondly, since the perturbation is small one can estimate the values of (c^(i), c^(o), p) on the perturbed boundary via Taylor's expansion, and further derive the boundary conditions for the first-order terms (similar to the calculations in Section <ref>, so we omit it here), 


∂_x c_0^(i)(0)+c_1,l^(i)(0)   =
∂_x c_0^(o)(0)+c_1,l^(o)(0),

∂_x^2 c_0^(i)(0)+∂_xc_1,l^(i)(0)   =
∂_x^2 c_0^(o)(0)+∂_xc_1,l^(o)(0),

∂_x p_0^(i)(0)+p_1,l^(i)(0)   =
0.


Finally, by solving the above boundary value problems for the first-order terms to get their explicit expressions:




c_1,l^(i)(x)   =-√()·c_B/√(+l^2)+√(1+l^2)e^√(+l^2)x,

c_1,l^(o)(x)   =-√()·c_B/√(+l^2)+√(1+l^2)e^-√(1+l^2)x,

p_1,l^(i)(x)   =G_0c_B/√()((1/√(λ)+1-1/√(λ+l^2)+√(1+l^2))e^l x+e^√(+l^2)x/√(+l^2)+√(1+l^2)).






 §.§ The evolution equations for the perturbation amplitude

In the last step, by utilizing the expression of p_0^(i)(x) and p_1,l^(i)(x), one can further derive the evolution equations of the perturbation amplitude from Darcy's law and Taylor's expansion. That is,








    ϵ^-1dϵ/dt   =-(∂^2_xp_0^(i)(0)+∂_xp_1,l^(i)(0))
       =G_0 c_B/√(λ)(√(λ)-l/√(λ)+1+l-√(λ+l^2)/√(λ+l^2)+√(1+l^2)):=E(λ,l).


If we replace the domain ℝ^2 by the following infinite-length tube:

    Ω = {(x,y) | (x,y) ∈ (-∞,+∞) × [-π,π]},

and post periodic boundary conditions for the y variable. The calculations do not change, except one needs to restrict l∈ℕ for the perturbation wave number. 

We plot E as in (<ref>) in Figure 1, from which we claim E(,l) possesses the following properties.




 (1) Given any integer l≥ 2, there exists unique _0^l>l^2 such that 
    
    E(_0^l,l)=G_0 c_B/√(_0)(√(_0)-l/√(_0)+1+l-√(_0+l^2)/√(_0+l^2)+√(1+l^2))=0.


    
   (2) Given any >1, there exists a unique l_0>0, which may not be an integer, such that E(, l_0)=0. Furthermore, E(,l)≠ 0 for any l≠ 0  or  l_0. 



    In the first part, the fact _0^l>l^2 holds, since one can easily check that E(,l)<0 for any ≤ l^2. The existence of _0^l can be shown by checking the limits 
    
    lim_→ 0E(,l)=-∞,  lim_→ +∞E(,l)=0^+.

    And the uniqueness can be checked by directly checking its derivative. 
    In the second part, the existence of l_0 was proved in <cit.> by checking the asymptote. Again, the uniqueness can be verified directly.
    
As we shall show Proposition <ref> allows us to justify the existence of nonsymmetric traveling waves solutions to (<ref>) in the above tube-like domain by applying the Crandall-Rabinowitz Theorem to a functional equation. 



§ BIFURCATION ANALYSIS


In this section, we justify the existence of non-symmetric traveling wave solutions in (<ref>) by relating it to the non-trivial bifurcation branch of a functional equation and conclude by a Crandall–Rabinowitz argument. The framework of utilizing Crandall–Rabinowitz theorem to study the bifurcation behavior of free boundary models was proposed initially by Friedman in <cit.>, then extensively employed to study the bifurcation phenomenon in different tumor growth models, see <cit.>. 
 
 














We begin with introducing the notations and reviewing the model. Recall that Ω denotes the tube-like domain defined in (<ref>), and Ω_0 be the symmetric tumor region in the same manner as (<ref>), and we further restrict l∈ℕ. Then, we consider a perturbed tumor region Ω_ϵ with respect to Ω_0, 

    Ω_ϵ={(x,y)| x≤ϵξ(y),   y∈[-π,π]},

where 0<ϵ≪1, and ξ(y) is a 2π periodic even function that characterize the boundary profile. Thus, the boundary can be represented as 

    ℬ_ϵ={(x,y)| x=ϵξ(y),   y∈[-π,π]}.


We employ (c^(i),c^(o),p^(i)) to denote the solution to equations (<ref>) and (<ref>), but with the domain and boundary replaced by Ω_ϵ and ℬ_ϵ, respectively. More precisely, (c^(i),c^(o),p^(i)) solves the following system:



-c^(i) + c^(i)    = 0,       in  Ω_ϵ,

-c^(o) + c^(o)    = c_B,       in  Ω∖Ω_ϵ,

c^(i)   =c^(o),       on  ℬ_ϵ,

∂/∂nc^(i)   =∂/∂nc^(o),       on  ℬ_ϵ,

c_ϵ=c^(i)+c^(o)   <∞,       in  Ω,

-p^(i)    = G_0 c^(i) ,       in  Ω_ϵ,

p^(i)    = 0,        on  ℬ_ϵ,

p^(i)   <∞,       in  Ω_ϵ.

Then, we further extend the pressure to the whole Ω, denote as p_ϵ, such that 

    p_ϵ=p^(i)+p^(o)

with p^(o) solves the following PDE:


-p^(o)    = G_0 c^(o) ,       in  Ω∖Ω_ϵ,

p^(o)    = 0,        on  ℬ_ϵ,

∂/∂np^(o)   =∂/∂np^(i),       on  ℬ_ϵ.   












Note that we introduce p^(o) for technical requirements  instead of physical, and we do not require p^(o)<∞ in Ω∖Ω_ϵ. Without the perturbation (when ϵ=0), the solution to (<ref>)-(<ref>) is given by

    p_0^(o)(x)=G_0 c_B(-x^2/2+√(λ)/√(λ)+1e^-x+λ-1/√(λ)(√(λ)+1)x-√(λ)/√(λ)+1),     x≥ 0.


To justify the existence of non-symmetric traveling wave solutions in (<ref>), we first introduce the following nonlinear functional map: 

    F(ξ,λ)=-∂ p_ϵ/∂ x|_ℬ_ϵ-σ_0(),

where, as before, the periodic even function ξ(y) stands for the boundary profile, and >0 presents the consumption rate. We emphasize that for the map F(ξ,) we view  as an index parameter and ξ as the independent variable. Regarding the right-hand side of (<ref>), p_ϵ stands for the pressure function associated with the profile ξ(y); and σ_0(), given in (<ref>), represents the traveling speed of the symmetric solution. Then we look for the solution to the functional equation

    F(ξ,λ)=0,

since these solutions correspond to the traveling wave solutions to (<ref>). In particular, the symmetric solutions correspond to the trivial solution (0,). We aim to show that for proper consumption rate , it can induce a non-trivial solution branch (ξ(y),), with ξ(y)≠ 0. The existence of such non-trivial bifurcation branches implies that (<ref>) admits symmetric breaking traveling wave solutions.

To find the non-trivial bifurcation branch to (<ref>), we adopt the framework proposed by Friedman, using the Crandrall-Rabinowitz theorem (see Theorem <ref>). Crandall Rabinowitz theorem is developed to study the bifurcation behavior in nonlinear equations. It provides conditions under which solutions branch off from a trivial solution in nonlinear operator equations. As a useful analysis tool, the Crandall Rabinowitz theorem has been widely applied to study the existence and stability of solutions in nonlinear systems, see, e.g., <cit.>. In particular, Friedman et al. first employed it to study the symmetric breaking solutions to free boundary problems and tumor growth models <cit.>.


To adopt Friedman's framework to our case, we carry out the main steps as follows:

    
  * determine the Fréchet derivative of F(ξ,λ) with respect to ξ∈ X (X is some function space to be specified later) on the line (0,λ), denote it as F_ξ(0,λ):=d/dξF(0,λ); and show the Fréchet derivative F_ξ(0,λ), as a bounded linear operator, can be further characterized in terms of an eigenvalue problem; we then find the complete basis of eigenfunctions with distinct eigenvalues (see equation (<ref>)). 
    
    
    
  * for each 2≤ l∈ℕ, we can determine a bifurcation point (0,λ_0^l) to the functional equation F(ξ,λ)=0 by utilizing the explicit expression of the eigenvalues (see (<ref>));
    
  * conclude {(0,λ_0^l)}_l=2^∞ are indeed bifurcation points to F(ξ,λ)=0 by verifying the Fréchet derivative at these points, F_ξ(0,λ_0^l), satisfies the bifurcation conditions in the Crandall-Rabinowitz theorem. More specifically, we verify our choice of λ_0^l ensures that only the l-th eigenvalue of F_ξ(0,λ_0^l) degenerates (equals zero). Also note that for every fixed λ_0, there is a spectural gap between this mode and other modes. So the limit can pass. 

To accomplish the first two steps above, we need to first look into the linearized system of (<ref>). We investigate it in the following subsection.



 §.§ The linearized system


We devote this section to studying the linearization of system (<ref>), which is closely related to our previous study in <cit.>. 

To begin with, we denote c_ϵ=c^(i)+c^(o), p_ϵ=p^(i)+p^(o), and analogously for

    c_0=c_0^(i)+c_0^(o),  
    c_1=c_1^(i)+c_1^(o),  
    p_0=p_0^(i)+p_0^(o),  
    p_1=p_1^(i)+p_1^(o).

Since the perturbation is small, i.e., 0<ϵ≪ 1, the solutions (c_ϵ, p_ϵ) possess the following asymptotic expansion with respect to ϵ:



    c_ϵ(x,y)   =c_0(x)+ϵc_1(x,y)+O(ϵ^2),

    p_ϵ(x,y)   =p_0(x)+ϵp_1(x,y)+O(ϵ^2),


with the zero-order terms represent the solutions to the unperturbed problem (the symmetric solution). 

Utilizing the expansion (<ref>), we can use Taylor expansion to evaluate c^(i) and ∂/∂ xc^(i) on the perturbed boundary in the following way 


    c^(i)|_ℬ_ϵ
       =c^(i)(ϵξ,y)
   =c_0^(i)(ϵξ)+ϵc_1^(i)(ϵξ,y)+O(ϵ^2)
   =c_0^(i)(0)+ϵξ∂/∂xc_0^(i)(0)+ϵc_1^(i)(0,y)+O(ϵ^2),

    ∂/∂xc^(i)|_ℬ_ϵ
       =∂/∂xc^(i)(ϵξ,y)
   =∂/∂xc_0^(i)(ϵξ)+ϵ∂/∂xc_1^(i)(ϵξ,y)+O(ϵ^2)
   =∂/∂xc_0^(i)(0)+ϵξ∂^2/∂x^2c_0^(i)(0)+ϵ∂/∂xc_1^(i)(0,y)+O(ϵ^2).


Similarly, for c^(o), p^(i), and p^(o) we have 


c^(o)|_ℬ_ϵ
   =c_0^(o)(0)+ϵξ∂/∂xc_0^(o)(0)+ϵc_1^(o)(0,y)+O(ϵ^2),

∂/∂xc^(o)|_ℬ_ϵ
   =∂/∂xc_0^(o)(0)+ϵξ∂^2/∂x^2c_0^(o)(0)+ϵ∂/∂xc_1^(o)(0,y)+O(ϵ^2),

p^(i)|_ℬ_ϵ
   = p^(i)_0(0)+ϵξ∂/∂xp^(i)_0(0)+ϵp^(i)_1(0,y)+O(ϵ^2),

∂/∂xp^(i)|_ℬ_ϵ   =∂/∂xp^(i)_0(0)+ϵξ∂^2/∂x^2p^(i)_0(0)+ϵ∂/∂xp^(i)_1(0,y)+O(ϵ^2),

p^(o)|_ℬ_ϵ
   = p^(o)_0(0)+ϵξ∂/∂xp^(o)_0(0)+ϵp^(o)_1(0,y)+O(ϵ^2),

∂/∂xp^(o)|_ℬ_ϵ   =∂/∂xp^(o)_0(0)+ϵξ∂^2/∂x^2p^(o)_0(0)+ϵ∂/∂xp^(o)_1(0,y)+O(ϵ^2).


Plugging the expansion (<ref>) into (<ref>), the zero-order terms are canceled out, and we collect the terms of order O(ϵ). Regarding the nutrient, the first order terms solve the following boundary value problem



-c_1^(i) + c_1^(i)    = 0,

-c_1^(o) + c_1^(o)    = 0,

c_1^(i)(0,y)   =c_1^(o)(0,y),

ξ·∂^2/∂x^2c_0^(i)(0)+∂/∂xc_1^(i)(0,y)
   =ξ·∂^2/∂x^2c_0^(o)(0)+∂/∂xc_1^(o)(0,y),

c_1^(i)(-∞,y)   <∞,

c_1^(o)(+∞,y)   <∞.


While, for pressure, the first order terms solve



-p^(i)_1    = G_0 c_1^(i),

ξ·∂/∂xp^(i)_0(0)+p^(i)_1(0,y)    = 0,

p^(i)_1(-∞,y)   <∞,

-p^(o)_1    = G_0 c_1^(o),

ξ·∂/∂xp^(o)_0(0)+p^(o)_1(0,y)    = 0,

ξ·∂^2/∂x^2p_0^(i)(0)+∂/∂xp_1^(i)(0,y)
   =ξ·∂^2/∂x^2p_0^(o)(0)+∂/∂xp_1^(o)(0,y).


On the other hand, the first-order terms capture the main reaction to the perturbation and, therefore, variable-separable,



    c_1^(i)(x,y)   =c_1^(i)(x)ξ(y),    c_1^(o)(x,y)   =c_1^(o)(x)ξ(y),

    p^(i)_1(x,y)   =p^(i)_1(x)ξ(y),    p^(o)_1(x,y)   =p^(o)_1(x)ξ(y).


In particular, when ξ(y)=cosly the first order terms reduce to the form of (<ref>), and the solutions to (<ref>) and (<ref>) are given by the single mode perturbation problem solved in Section <ref>. Note that through we did not provide the expression of p^(o)_1 in (<ref>), it is solvable via (<ref>) once c_1^(o), p_0^(i), p_0^(o), and p_1^(i) are determined. 

In the following sections, we will show that the above 
linearization study helps us to characterize the Fréchet derivative, F_ξ(0,), for any >0. 





 §.§ Derivation and characterization of F_ξ(0,)

In this section, we further determine and characterize the Fréchet derivative F_ξ(0,) based on the calculations in the previous subsection. To begin with, we introduce the following Banach spaces



    X^m+α   ={f(y)∈C^m+α: f(y) is 2π periodic and even};

    X_1^m+α   =closure of the linear space spanned by  {cosjy, j=0,1,2,…}  in  X^m+α.


Note that all modes are included. Thus, any ξ(y)∈ X_1^m+α can be represented as Fourier series. Then, to determine the Fréchet derivative F_ξ(0,), one needs to justify the expansion (<ref>) rigorously, which is equivalent to proof of the following two lemmas.


If ξ belongs to C^3+α(ℝ), c_ϵ and p_ϵ be the solution to the system (<ref>), then

    ‖ c_ϵ-c_0‖_C^1+α(Ω)   ≤ C|ϵ|‖ξ‖_C^3+α(ℝ),
    ‖ p_ϵ-p_0‖_C^3+α(Ω)   ≤ C|ϵ|‖ξ‖_C^3+α(ℝ),

where C is a constant independent of ϵ, and (c_0, p_0) stands for the unperturbed solutions given in (<ref>) and (<ref>).



If ξ belongs to C^3+α(ℝ), c_ϵ and p_ϵ be the solution to (<ref>), then

    ‖ c_ϵ-c_0-ϵĉ_1‖_C^1+α(Ω)   ≤ C|ϵ|^2‖ξ‖_C^3+α(ℝ),
    ‖ p_ϵ-p_0-ϵp̂_1‖_C^3+α(Ω)   ≤ C|ϵ|^2‖ξ‖_C^3+α(ℝ),

where C is a constant independent of ϵ, and (c_0, p_0) stands for the unperturbed solutions, (ĉ_1, p̂_1) corresponds to the Hanzawa transformation (see (<ref>)) of the first order terms.

The proof of the above lemmas is standard but cumbersome. To avoid the reader's distraction, we provide a sketch proof in Section <ref>. For now, we directly use them to determine the Fréchet derivative of F.

The calculations in Section <ref> yields

    F(ξ,λ)
           =-∂ p_ϵ/∂ x|_ℬ_ϵ-σ_0()
       =-(∂_x^2 p_0(0)+∂_x p_1(0))ϵξ+O(ϵ^2).

Furthermore, Lemma <ref> and Lemma <ref> implie F maps from X^m+α×ℝ to X^m+α for any m≥ 3. And, according to (<ref>), for ξ(y)=cosly the above identity reduces to

    F(cosly,λ)
       =-(∂_x^2 p_0^(i)(0)+∂_x p_1,l^(i)(0))ϵcosly+O(ϵ^2)
       =E(λ,l)ϵcosly+O(ϵ^2),

equivalently,

    F(cosly,λ)-F(0,λ)-E(λ,l)ϵcosly_X^m+α≤ Cϵ^2.

Thus, given any >0 the Fréchet derivative F_ξ(0,), as a bounded linear operator, is fully characterized by the following eigenvalue problem

    [F_ξ(0,)]cosly=E(,l) cosly,

where the eigenvalue E(,l)∈ℝ. In the following subsection, we show that based on the above eigenvalue problem and the properties of E(,l) in Proposition <ref>, we can determine bifurcation points and further conclude the existence of non-trivial bifurcation branches to equation (<ref>).
















 §.§ Existence of non-trivial bifurcation branches
 In the seminal work <cit.>, Friedman et al. employed the Crandall-Rabinowitz theorem to show the existence of non-radial symmetric solutions to a tumor growth model developed from <cit.>. Although, as discussed in the introduction section, (<ref>) is derived from the incompressible limit of PME and, therefore, essentially different from the models developed from <cit.>, the bifurcation analysis framework established by Friedman remains applicable. We carry out the bifurcation analysis in this subsection. To begin with, for the reader's convenience, we present the Crandall-Rabinowitz theorem below.



    Let W,Z be real Banach spaces and (w,μ) a C^p map, p≥ 3, of a neighborhood (0,μ_0) in W×ℝ into Z. For any ω∈ W and ν∈ℝ, _w(ω,ν):=d/dwF(ω,ν)∈ℒ(W, Z), where ν is viewed as a parameter. Suppose
    
        
  * (0,μ)=0 for all μ in a neighborhood of μ_0,
        
        
  * The kernel space of the partial derivative _w:W→ Z at (0,μ_0) is of one dimensional spanned by w_0∈ W, i.e.,  Ker{_w(0,μ_0)}=span{w_0}.
        
  * The range of _w(0,μ_0) has codimension 1, i.e., dim{Z/Z_1}=1 with Z_1=Img{_w(0,μ_0)}⊆ Z.
        
  * The derivative _wμ(w,r), w∈ W, r∈ℝ, satisfies [_wμ(0,μ_0)]w_0∉ Z_1.
    
    Then, (0,μ_0) is a bifurcation point of the equation (w,μ)=0 in the following sense: In a neighborhood of (0,μ_0), the set of solutions of (w,μ)=0 consists of two C^p-2 smooth curves 𝒞_1 and 𝒞_2 which intersect only at the point (0,μ_0). Moreover, 𝒞_1 is the curve (0,μ) and 𝒞_2 can be parameterized as follows:
    
    𝒞_2: (w(ϵ),μ(ϵ)), ϵ small, (w(0),μ(0))=(0,μ_0), w'(0)=w_0.


To apply Theorem <ref> to the nonlinear map (<ref>), the main ingredient is to find a bifurcation point (0,_0) such that the partial derivative F_ξ(0,_0) satisfies the assumptions in Theorem <ref>. According to Friedman's framework, it is crucial to check this _0 leads the eigenvalue to (<ref>) vanish in only one direction, i.e., E(_0,l)=0 holds for only one specific l∈ℕ. We show that this can be done by utilizing Proposition <ref>. And therefore, (<ref>) posses non-symmetric traveling wave solutions. We summarize this main result in the following theorem.


Consider the nonlinear map (<ref>), which maps X_1^3+α×ℝ to X_1^3+α. Assume 0<ϵ≪ 1. Then for each integer l≥ 2, there exists a _0^l>l^2 such that (0,_0^l) is a bifurcation point to F(ξ,)=0 in the sense of: In a neighborhood of (0,_0^l), the set of solutions of F(ξ,)=0 consists of two smooth curves 𝒞_1 and 𝒞_2 which intersect only at the point (0,_0^l). Moreover, 𝒞_1 is the curve (0,) and 𝒞_2 can be parameterized as follows:
    
    𝒞_2: (ξ(ϵ),(ϵ)), ϵ small, (ξ(0),(0))=(0,_0^l), ξ'(0)=cosly.



According to Proposition <ref>, given any integer l≥ 2, we can find an unique _0^l>l^2 such that E(_0^l,l)=0. Then we show (0,_0^l) is indeed a bifurcation point to (<ref>) by verifying the map F(ξ,) indeed satisfies the conditions for applying Theorem <ref> with the setting W=X_1^3+α, Z=X_1^2+α, w=ξ, w_0=cosly, and μ=λ.

For the differentiability of F, it is equivalent to establishing the regularity of the corresponding PDEs. Firstly, note that the structure of the PDEs guarantees that F maps even 2π-periodic functions to even 2π-periodic functions. Then, Lemma <ref> and Lemma <ref> imply that F maps X_1^3+α×ℝ into X_1^3+α. Secondly, by using classical elliptic estimates and Sobolev imbedding theory, one can justify F(ξ,μ) is differentiable to any order by repeating the process in the same manner as Lemma <ref> and Lemma <ref>. Therefore, F(ξ,) is C^p with p≥ 3.

Next, we verify the assumptions (1) to (4) hold for F(ξ,) at the point (0,_0^l). Firstly, (1) obviously holds since these trivial solutions correspond to the symmetry solutions. Regarding assumptions (2) and (3), recall that F_ξ(0,) as a bounded linear operator is characterized by the eigenvalue problem (<ref>). Thus, to check (2) and (3), it is sufficient for us to check that our choice of _0^l ensures: 

    E(_0^l,j)≠ 0  for any   j≠ l;  and   E(_0^l,l)=0.

Based on Proposition <ref> and the way of chosen _0^l, condition (<ref>) indeed holds. Also note that for every fixed _0>0, there is a spectral gap between the l-th mode and other modes. So the limit can pass. Finally, for assumption (4), it is sufficient for us to show ∂_E(_0^l,l)≠ 0. Indeed,

    ∂_E(_0^l,l)
       =1/2√(_0^l)(l+1/√(_0^l)(√(_0^l)+1)^2-l+√(1+l^2)/√(_0^l+l^2)(√(_0^l+l^2)+√(1+l^2))^2)
       >l+1/2√(_0^l)(1/√(_0^l)(√(_0^l)+1)^2-1/√(_0^l+l^2)(√(_0^l+l^2)+√(1+l^2))^2)>0,

where we used condition (<ref>) to derive the first identity. Thus,
 
    [F_ξ(0,λ_0^l)]cosly=[∂_E(_0^l,l)]cosly∉Img[F_ξ(0,λ_0^l)].

By now, we have finished verifying all the assumptions in the Crandall-Rabinowitz theorem. Therefore, (0,_0^l) is a bifurcation point to (<ref>) and generates a non-trivial solution branch. As we interpreted before, this non-trivial solution branch corresponds to the non-symmetric traveling wave solutions to (<ref>).





 §.§ Justification of the expansion



We devote this section to the proof of Lemma <ref> and Lemma <ref>. 









To begin with, recall that Ω is the tube-like domain defined in (<ref>). Ω_0 and Ω_ϵ corresponds to the unperturbed and perturbed tumor region respectively. For concision, we denote the complementary sets as

    Ω_0^c=Ω∖Ω_0,    Ω_ϵ^c=Ω∖Ω_ϵ.










Now, we provide the proof of Lemma <ref> as follows. 











Note that if we denote c^δ=c_ϵ-c_0, then it satisfies


   -Δc^δ+λc^δ   =0,       in  Ω_ϵ∩Ω_0:=Ω_1;
 
   -Δc^δ+λc^δ   =(1-λ)c_0-c_B,       in  Ω_ϵ∩Ω_0^c:=Ω_2;
 
    -Δc^δ+ c^δ   =(λ-1)c_0+c_B,       in  Ω_ϵ^c∩Ω_0:=Ω_3;
 
    -Δc^δ+ c^δ   =0,       in  Ω_ϵ^c∩Ω_0^c:=Ω_4.


Write them in a single equation, one gets

    -Δ c^δ+(λ·χ_Ω_ϵ+χ_Ω_ϵ^c)c^δ=((1-λ)c_0-c_B)·(χ_Ω_2-χ_Ω_3),    in  Ω.

Observe the facts that λχ_Ω_ϵ+χ_Ω_ϵ^c can be treated as a function in L^∞(Ω), c_0 has already been solved explicitly on Ω. Furthermore, the areas |Ω_2| and |Ω_3| are both bounded by ϵ‖ξ‖_C^3+α(ℝ). Then, the classical W^2,p estimate of elliptic equations and Sobolev embedding theory together yield the first inequality in Lemma <ref>. More precisely, for any p>2 and α=1-2/p one has

    ‖ c^δ‖_C^1+α(Ω)≤‖ c^δ‖_W^2,p(Ω)≤‖((1-λ)c_0-c_B)·(χ_Ω_2-χ_Ω_3)‖_L^p(Ω)≤ C|ϵ|‖ξ‖_C^3+α(ℝ).

Finally, send p→∞ to complete the proof.

For the second inequality in Lemma <ref>, one can easily write down the equation and boundary condition for p^δ=p_ϵ-p_0 on Ω̅_ϵ, that is


-Δp^δ   =G_0 c^δ,       in  Ω_ϵ,

p^δ   =-p_0,       on  ℬ_ϵ.


Note that p_0=p_0^(i)+p_0^(o) has already been solved on Ω, in particular for p_0|_ℬ_ϵ. Thus, by using Schauder estimate one has

    ‖ p^δ‖_C^3+α(Ω̅_ϵ)≤ G_0‖ c^δ‖_C^1+α(Ω_ϵ)+‖ p_0‖_C^3+α(ℬ_ϵ)≤ C|ϵ|‖ξ‖_C^3+α(ℝ).


Next, observe the fact that (c^(i), c^(o), p^(i), p^(o)) are defined
on Ω_ϵ or Ω_ϵ^c respectively. However, the first-order terms (c_1^(i), c_1^(o), p_1^(i), p_1^(o)) are only defined on Ω_0 or Ω_0^c. Therefore, we need to transform them to Ω_ϵ or Ω_ϵ^c by Hanzawa transformation _ϵ, which is defined as follows:

    (x,y)=_ϵ(x',y')=(x'+ (x')ϵξ,y'),

where  is defined by:

    ∈ C^∞,  (ζ)={[                   0,  if  |ζ|≥3/4δ,;                   1,  if  |ζ|<3/4δ, ].,   with  |d^k/dζ^k|<C/δ^k,

where δ is a small positive scalar. Thus, _ϵ maps Ω_0 onto Ω_ϵ, and maps Ω_0^c onto Ω_ϵ^c. We further denote



    ĉ_1^(i)(x,y)   =c_1^(i)(_ϵ^-1(x,y)),    ĉ_1^(o)(x,y)   =c_1^(o)(_ϵ^-1(x,y)),

    p̂_1^(i)(x,y)   =p_1^(i)(_ϵ^-1(x,y)),    p̂_1^(o)(x,y)   =p_1^(o)(_ϵ^-1(x,y)).


Now, we turn to the proof of Lemma <ref>. The detail of the proof is cumbersome, but the idea is quite simple and in the same manner as the proof of Lemma 3.1. Therefore, we only provide a sketch of it.










The proof is similar to that of Lemma <ref>. Denote c^δ:= c_ϵ-c_0-ϵĉ_1 and similarly for p^δ:=p_ϵ-p_0-ϵp̂_1. Then, one is able to write out the equation for c^δ on the whole Ω. Then, employ W^2,p estimate of the elliptic equations and the embedding theory to obtain the estimate for the nutrient first, as we did in Lemma <ref>. However, to do this, one needs to compute the first and second derivatives of ĉ_1 with respect to (x,y), which further requires us to consider the change of variables induced by the Hanzawa transformation. This process is cumbersome but standard. Therefore, we refer the reader to Theorem 4.5 in <cit.> for a similar proof. Once the estimate of the nutrient is obtained, one can further obtain the pressure estimate by Schauder estimate in the same manner as Lemma <ref>.




















































































































































































































































































































































































































































§ ACKNOWLEDGMENTS

The work of Y.F. is supported by the National Key R&D Program of China, Project Number 2021YFA1001200. J.-G.L is supported by NSF under award DMS-2106988. The work of Z.Z. is supported by the National Key R&D Program of China, Project Number 2021YFA1001200, and NSFC grant number 12031013, 12171013. We thank Jiajun Tong (BICMR) for the helpful discussions. 

















































































 







plain
