{'fraction_non_alphanumeric': 0.058664359606818095, 'fraction_numerical': 0.028421013879173875, 'mean_word_length': 4.002486187845304, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Cahn--Hilliard--Darcy model for tumour growth with chemotaxis and active transport

Using basic thermodynamic principles we derive a Cahn--Hilliard--Darcy model for tumour growth including nutrient diffusion, chemotaxis, active transport, adhesion, apoptosis and proliferation. The model generalises earlier models and in particular includes active transport mechanisms which ensure thermodynamic consistency. We perform a formally matched asymptotic expansion and develop several sharp interface models. Some of them are classical and some are new which for example include a jump in the nutrient density at the interface. A linear stability analysis for a growing nucleus is performed and in particular the role of the new active transport term is analysed. Numerical computations are performed to study the influence of the active transport term for specific growth scenarios.


Introduction
In the last decades the understanding of tumour related illnesses has undergone a swift development. Nowadays tumour therapy can be adapted to the genetic fingerprint of the tumour, resulting in a "targeted therapy" that has dramatically improved the prognosis of many illnesses. While some important mutations in tumour genomes have been identified and exploited by modern tumour drugs, basic growth behaviours of tumours are still far from being understood, e.g. angiogenesis and the formation of metastases. The complexity of oncology has also attracted increasing interest of mathematicians, who are trying to find the appropriate equations to provide additional insights in certain aspects of tumour growth, see for example [6] and [15]. In this paper we want to introduce a new diffuse interface model for tumour growth, and compare the resulting system of partial differential equations to some other recent contributions [13,14,22,28,29,30,32,37,38,43,44].
In order to obtain a tractable system of partial differential equations, we will in this paper neglect some effects which could be addressed in further research and which then would lead to more complete theories. From a medical point of view we will hence make the following assumptions as foundations for our further considerations: 1. Tumour cells only die by apoptosis. Hence we neglect the possibility of tumour necrosis, where we would have to take account of the negative effects of chemical species from the former intracellular space on the surrounding tumour cells.
2. The tissue around the tumour does not react to the tumour cells in any active way.
In particular, we neglect any response of the immune system to the tumour tissue.
3. Larger tumour entities are actually enforcing blood vessel growth towards themselves by secreting vessel growth factors. This is a phenomenon that could be addressed in future in a generalised model. 4. We postulate the existence of an unspecified chemical species acting as a nutrient for the tumour cells. This nutrient is not consumed by the healthy tissue. We will also introduce terms which will reflect chemotaxis, which is the active movement of the tumour colony towards nutrient sources. Additionally, the introduction of chemotaxis will also lead to the opposite process, meaning that the nutrient is moving towards the nearby tumour cells. As we will point out later, this could be seen as a correlate of a nutrient uptake mechanism.
Here we state a slightly simplified version of the general system, which will be derived in Section 2 from thermodynamic principles. We will derive and analyse a two-component mixture model of tumour and healthy cells, whose behaviour is governed by the system div v = αΓ, Here, v denotes the volume-averaged velocity of the mixture, p denotes the pressure, σ denotes the concentration of an unspecified chemical species that serves as a nutrient for the tumour, ϕ ∈ [−1, 1] denotes the difference in volume fractions, with {ϕ = 1} representing unmixed tumour tissue, and {ϕ = −1} representing the surrounding healthy tissue, and µ denotes the chemical potential for ϕ. The particular simple form of (1.1a) is different to earlier modelling attempts and is based on the fact that we use volumeaveraged velocities. The positive constants K, β, P, A, and C denote the permeability, surface tension, proliferation rate, apoptosis rate, and consumption rate, respectively. The constants ρ S and α are related to the densities of the two components (see (2.32) below), in particular, for the case of matched densities we have α = 0. Meanwhile m(ϕ) and n(ϕ) are nonnegative mobilities for ϕ and σ, respectively, and Ψ(⋅) is a potential with two equal minima at ±1. In addition, we choose h as an interpolation function with h(−1) = 0 and h(1) = 1. The simplest choice is given as h(ϕ) = 1 2 (1 + ϕ). We denote χ σ ≥ 0 as the diffusivity of the nutrient, and χ ϕ ≥ 0 can be seen as a parameter for transport mechanisms such as chemotaxis and active uptake (see below for more details). Finally, the parameter ε is related to the thickness of the interfacial layers present in phase field systems. The system (1.1) is a Cahn-Hilliard-Darcy system coupled to a convection-diffusion-reaction equation for the nutrient.
(1.1a) and (1.1b) model the mass balance using a Darcy-type system, and in the situation of unmatched densities (α ≠ 0), the gain and loss of volume resulting from the mass transition Γ leads to sources and sinks in the mass balance. In (1.1c) and (1.1d), ϕ is governed by a Cahn-Hilliard type equation with additional source terms. The mass transition from the the healthy cells to the tumour component and vice versa is described in (1.1f), where tumour growth/proliferation is represented by the term Pσh(ϕ), and the process of apoptosis is modelled by the term Ah(ϕ). In (1.1e), the nutrient is subjected to an equation of convection-reaction-diffusion type, and the term Cσh(ϕ) represents consumption of the nutrient only in the presence of the tumour cells. As in [10], we could also consider the situation where the tumour possesses its own vasculature and the nutrient may be supplied to the tumour via a capillary network at a rate B(σ−σ B ), where σ B is the constant nutrient concentration in the vasculature and B is the bloodtissue transfer rate which might depend on ϕ and x. This leads to the following nutrient balance equation instead of (1.1e) ∂ t σ + div (σv) = div (n(ϕ)(χ σ ∇σ − χ ϕ ∇ϕ)) − Cσh(ϕ) + B(σ B − σ).
Under appropriate boundary conditions the system (1.1) allows for an energy inequality (see (2.27) below) and we believe that this inequality will allow the well-posedness of the above system to be rigorously shown.
We now motivate the particular choices for the modelling of proliferation, apoptosis, chemotaxis, and mass transition in (1.1).
• In (1.1f), we obtain that Γ = Pσ −A holds in the tumour region {ϕ = 1}. The implicit assumption that the tumour growth is proportional to the nutrient supply can be justified by the fact that malign tumours have the common genetic feature that certain growth inhibiting proteins have been switched off by mutations. Hence, we can assume that while in healthy cells the mitotic cycle is rather strictly inhibited, tumour cells often show unregulated growth behaviour which is only limited by the supply of nutrients.
Moreover, implicit in the choice of zero mass transition Γ = 0 in the healthy region {ϕ = −1} is the assumption that the tumour proliferation rate is more significant than that of the healthy tissue.
• In (1.1c) and (1.1e), the fluxes for ϕ and σ are given by respectively. It has been pointed out by Roussos, Condeelis and Patsialou in [40] that the undersupply of nutrient induces chemotaxis in certain tumour entities. This is reflected in the term m(ϕ)∇(χ ϕ σ) of q ϕ , which drives the cells towards regions of high nutrient.
On the other hand, we note that the term n(ϕ)∇(χ ϕ ϕ) in q σ drives the nutrient to regions of high ϕ, i.e., to the tumour cells, which indicates that the nutrient is actively moving towards the tumour cells. This may seem to be counter-intuitive at first glance. However, this term will only contribute to the equation significantly in the vicinity of the interface between the tumour and healthy cells. This allows the interpretation that the term n(ϕ)∇(χ ϕ ϕ) reflects active transport mechanisms which move the nutrient into the tumour colony. Here we use the term "active transport" in the biological sense in order to indicate that some kind of mechanism is needed to maintain the transport (in contrast to passive transporters which are driven only by the concentration gradient of the substance). The additional mechanism allows cells to establish persisting concentration differences between different compartments. In particular we can expect that tumours, which have these active transporters on their cell membrane, are not dependent on diffusion but can establish high concentration of the vital nutrient even against the nutrient concentration gradient.
Here we briefly give an example, where mechanisms like this have already been observed: Malign tumour cells often have a significantly increased need for glucose, a fact that is sometimes referred to as the Warburg effect. As a consequence of several mutations in the tumour genome, these cells can adapt to their high rate of glucose consumption in several ways. Apart from angiogenesis, which leads to a well perfused tumour environment providing large amount of glucose, the tumour cells can also express (i.e., build) more glucose transporters, which provide an improved glucose transport through the cell membrane. Recently, both passive glucose transporters, so-called GLUT proteins, and active glucose transporters called SGLTs, have been observed on the cell membrane of several tumour entities. For a more detailed description regarding the GLUT transporters we refer to [11], whereas SGLT expression of tumours has been described by [31] and [41]. In the system we will derive in this paper, the existence of passive nutrient transporters like the GLUTs is implicitly assumed by including nutrient diffusion. Apart from that, it will become more obvious in the corresponding sharp interface system that the term n(ϕ)∇(χ ϕ ϕ) represents an active nutrient transport towards the tumour.
We note that in (1.1), the mechanism of chemotaxis and active transport are connected via the parameter χ ϕ . In principle, it is possible to decouple the two mechanisms. In order to do so, we introduce the following choice for the mobility n(ϕ) and diffusion coefficient χ σ (see also Section 3.3.3 below): For λ > 0 and a non-negative mobility D(ϕ), we set Then, the corresponding fluxes for ϕ and σ are now given as For this choice, we can switch off the effects of active transport by sending λ → 0, while preserving the effects of chemotaxis. We now compare the new model (1.1) and some of the previous diffuse interface models in the literature: • An important difference to other models is that we take a volume-averaged velocity which leads to the simple form div v = αΓ, as a consequence of mass balance. In other models, this equation has to be replaced by a more complicated transport equation.
• Another significant difference is the presence of the term − div (n(ϕ)χ ϕ ∇ϕ) in (1.1e), which only appears in cases where chemotaxis and active transport are taken into account. As we have pointed out before, it represents active nutrient transport towards the tumour. The corresponding nutrient equations in [13,14,22,28,30,38,44] do not include an equivalent term. However, we point out that this active transport mechanism is present in the nutrient equation of [29], who however used different source terms and no Darcy-flow contributions.
• Our choice of the mass transition term Γ in (1.1f) can also be found in [14,28,38,44]. Alternatively, one may consider equations of the form where the chemical potential µ enters as a source term for the equations of ϕ and σ. Here χ ≥ 0 is a constant, and P (⋅) denotes a non-negative proliferation function. This type of mass transition term appears in [29] and in [13,22,30] with χ = 0.
• In [32,37,44], the nutrient does not enter the Darcy law for v like in (1.1b).
In the diffuse interface model (1.1), the parameter ε is related to the thickness of the interfacial layer, which separates the tumour cell regions {ϕ = 1} and the healthy cell regions {ϕ = −1}. Hence, it is natural to ask if a sharp interface description of the problem will emerge in the limit ε → 0. This means in the limit the interface between the tumour cells and the healthy cells is represented by a hypersurface of zero thickness.
For convenience, suppose we take the mobilities m(ϕ) = m 0 , n(ϕ) = n 0 to be constant. A formally matched asymptotic analysis will yield the following sharp interface limit from (1.1) (see Section 3 for more details): Let Ω T and Ω H denote the tumour cell region and the healthy cell region, respectively, which are separated by an interface Σ. Then it holds that Here, γ is a constant related to the potential Ψ (see (3.17) below), V denotes the normal velocity of Σ, κ is the mean curvature of Σ, [f ] T H denotes the jump of f from Ω T to Ω H across Σ (see (3.16)), and ν is the outward unit normal of Σ, pointing towards Ω T .
In comparison to the formal sharp interface limits of [14,30,44], the most significant difference is the jump condition (1.4e) 2 . Let us remark on its physical meaning. Let σ T and σ H denote the limiting values of the nutrient on the interface Σ from tumour cell regions and from the healthy cell regions, respectively. Then, (1.4e) 2 implies that Thus, if χ ϕ is positive, then (1.4e) 2 tells us that the tumour cells will experience a higher level of nutrient concentration than the healthy cells on the interface, which reflects the effect of the active transport mechanism in (1.1e), attracting nutrients from the healthy cell regions into the tumour. If we consider the fluxes (1.3) in (1.1), then one obtains the sharp interface model (1.4) with the following modification (see Section 3.3.3 for more details): Instead of (1.4e) 2 , we now have In particular, the parameter λ only enters explicitly in the jump condition for σ, which relates to the above discussion regarding the physical interpretation of (1.4e) 2 . Note that λ ≠ 0 is a consequence of the active transport term we have discussed in the phase field model. In the matched asymptotics expansion, this term directly leads to a jump of the nutrient concentration at the tumour interface. Therefore, we obtain exactly the situation one would expect from active transport mechanisms: Close to the tumour surface, we observe a higher nutrient concentration inside the tumour than on the outside of the tumour, a situation that is only possible due to the transporter molecules. Hence it should be considered if λ could be referred to as a density parameter for the active transport proteins. The plan of this paper is as follows: In Section 2 we derive the new phase field model from thermodynamic principles and compare with previous phase field models of tumour growth in the literature. In Section 3 we perform a formal asymptotic analysis to derive certain sharp interface models of tumour growth. In Section 4 we investigate the stability of radial solutions to a particular sharp interface model via a linear stability analysis, and highlight the effect of the active transport parameter on the stability. In Section 5 we present quantitative simulations for radially symmetric solutions and qualitative simulations for more general scenarios.

Model Derivation
Let us consider a two component mixture consisting of tumour and healthy cells in a bounded domain Ω ⊂ R d , d = 1, 2, 3.
We denote the first component as the component of healthy tissues, and the second component as the tumour tissues. Let ρ i , i = 1, 2, denote the actual mass of the component matter per volume in the mixture, and letρ i , i = 1, 2, be the mass density of a pure component i. Then, ρ ∶= ρ 1 + ρ 2 denotes the mixture density (which is not necessarily constant), and we define the volume fraction of component i as We expect that physically, ρ i ∈ [0,ρ i ] and thus u i ∈ [0, 1]. In addition to the considerations stated in Section 1, we make the following modelling assumptions: • There is no external volume compartment besides the two components, i.e., • We allow for mass exchange between the two components. Growth of the tumour is represented by mass transfer from component 1 (healthy tissues) to component 2 (tumour tissues), while tumour cells are converted back into the surrounding healthy tissues when they die.
• We choose the mixture velocity to be the volume-averaged velocity: where v i is the individual velocity of component i.
• We model a general chemical species which is treated as a nutrient for the tumour tissues. Its concentration is denoted by σ and it is transported by the volumeaveraged mixture velocity and a flux J σ .

Balance laws
The balance law for mass of each component reads as Observe that by (2.1), we can write (2.4) in the following way: For i = 1, 2, We see that by (2.2), (2.3), and (2.5), We introduce the fluxes: Then, we see that and so, upon adding the equations in (2.4) we obtain the equation for the mixture density: We now want to derive an equation for the phase field variable ϕ. Recalling ρ i =ρ i u i , we obtain from (2.5) that We define the order parameter ϕ as the difference in volume fractions: then, subtracting the equation for u 1 from the equation for u 2 , and using (2.7), we obtain the equation for ϕ: We point out that from the constraint (2.2), we obtain Thus, the region of the tumour tissues is represented by {x ∈ Ω ∶ ϕ = 1} and the region of healthy tissues is represented by {x ∈ Ω ∶ ϕ = −1}. In particular, the mixture density ρ can be expressed as a linear function of ϕ: For the nutrient, we postulate the following balance law: where S denotes a source/sink term for the nutrient. In addition, σv models the transport by the volume-averaged velocity and J σ models other transport mechanisms like diffusion and chemotaxis.

Energy inequality
We postulate a general energy density of the form: e(ϕ, ∇ϕ, σ) = f (ϕ, ∇ϕ) + N (ϕ, σ). (2.14) Here, we neglected inertia effects, and so the kinetic energy does not appear in e. Instead we refer the reader to [43] for the derivation of a model that includes inertia effects, leading to a Navier-Stokes-Cahn-Hilliard version of (1.1). The first term f in (2.14) accounts for interfacial energy and unmixing tendencies, while the second term N describes the chemical energy of the nutrient and energy contributions resulting from the interactions between the tumour tissues and the nutrient. The latter will, for example, lead to chemotatic effects which are of particular interest as they result in the tumour tissue growing towards regions with high nutrient concentration.
In the following, we will consider f to be of Ginzburg-Landau type: For constants A, B > 0, we choose where Ψ(s) is a potential with equal minima at s = ±1.
We will now derive the diffuse interface model based on a dissipation inequality for balance laws with source terms which has been used similarly by Gurtin [25,26] and Podio-Guidugli [39] to derive phase field and Cahn-Hilliard type equations. These authors used the second law of thermodynamics which in an isothermal situation is formulated as a free energy inequality. We also refer to Chapter 62 of Gurtin, Fried, and Anand [27] for a detailed discussion of situations with source terms. The second law of thermodynamics in the isothermal situation requires that for all volumes V (t) ⊂ Ω which are transported with the fluid velocity the following inequality has to hold (see [25,26,27,39]): where ν is the outer unit normal to ∂V (t) and J e is an energy flux yet to be specified. Following [27], we postulate that the source terms Γ v , Γ ϕ and the nutrient supply (−S) carry with them a supply of energy described by for some c v , c ϕ and c S yet to be determined. Using the transport theorem and the divergence theorem, we obtain the following local form We now use the Lagrange multiplier method of Liu and Müller, see for example Section 2.2 of [2] and Chapter 7 of [36]. Let λ v , λ σ and λ ϕ denote Lagrange multipliers for the divergence equation (2.6), the nutrient equation (2.13) and the order parameter equation (2.11). We require that the following inequality holds for arbitrary where we used the notation as the material derivative of ϕ with respect to v. Using the identity We use the following notation: Applying the product rule to the divergence term in (2.19), we then obtain (2.20) Employing the following identities

Constitutive assumptions and the general model
We are now seeking for a model fulfilling the second law of thermodynamics in the version of a dissipation inequality stated in Section 2.2. We do not aim for the most general model but will state certain constitutive assumptions which take the most relevant effects into account. We hence make the following constitutive assumptions: where n(ϕ) and m(ϕ) are non-negative mobilities. We introduce a pressure-like function p and choose and for a positive constant K, Eq. (2.22a) makes a constitutive assumption for the energy flux J e which guarantees that the divergence term in (2.21) vanishes. It contains classical terms like µJ and N ,σ J σ which describe energy flux due to mass diffusion and the non-classical term B∂ t ϕ∇ϕ which is due to moving phase boundaries, see also [3,4] where this term is discussed. The last term in (2.22a) will result in energy changes due to work by macroscopic stress, compare [2]. Meanwhile, (2.22b), (2.22c), (2.23) and (2.24) are considered in order for the right hand side of (2.21) to be non-positive for arbitrary values of (ϕ, σ, v, . We mention that (2.24) is a Darcy law with force (µ − N ,ϕ )∇ϕ. Thus, the model equations for tumour growth are Supplemented with the boundary conditions then the above model satisfies the following energy equality: This follows from integrating (2.21) over Ω and using the definition of −D from (2.18), the constitutive assumptions (2.22), (2.23), and (2.24), and applying the divergence theorem.
Here, we have not prescribed boundary conditions for N ,σ and v. We will look at suitable boundary conditions for them later. We point out that using (2.8), (2.12), (2.25a), (2.25c), and the definition of Γ ϕ and Γ v , we obtain Thus, we identify J = − ρ 2 −ρ 1 2 m(ϕ)∇µ, and the equation for ρ becomes Remark 2.1 (Reformulations of the pressure and Darcy's law). In the above derivation, we may consider the following pressure-type functions: We point out that (2.29) can also be obtained from the momentum balance of the Navier-Stokes-Cahn-Hilliard equations by neglecting the inertia terms and replacing the viscous term with a multiple of the velocity. This is consistent with the classical derivation of Darcy's law. Meanwhile, in (2.30) we have the gradient of the primary variables (ϕ, σ) multiplied by their corresponding chemical potentials (µ, N ,σ ), and vice versa in (2.31) (compare with the interfacial term K in Section 3 of [2] and Eq. (2.34) of [33]). It is common to reformulate the pressure as above to obtain equations of momentum balance in the Navier-Stokes-Cahn-Hilliard equations or the Cahn-Hilliard-Darcy equations that are more amenable to further analysis. See for instance [1,20].

Zero velocity, zero excess of total mass and equal densities
Suppose the volume-averaged mixture velocity v is zero, the excess of total mass Γ 1 + Γ 2 is zero and the densities are equal. Then, substituting v = 0 in (2.34) and neglecting the Darcy system (2.34a,b), we obtain

Boundary conditions for velocity and nutrient
For the nutrient, we may prescribe a Robin type boundary condition: where c ≥ 0 is a constant, and σ ∞ denotes a given supply at the boundary. When c = 0, we obtain the zero flux boundary condition: If we formally send c → ∞, then we obtain the Dirichlet boundary condition: We may consider a boundary condition for the normal component of the velocity (which corresponds to a Neumann boundary condition for the pressure):  However, if the mass source terms Γ i depend on ϕ, σ or µ, then considering (2.40) as a boundary condition would imply that ϕ, µ and σ have to satisfy Alternatively, we can prescribe a boundary condition for the pressure. Recall the reformulated pressurep and the Darcy's law (2.30). We can prescribe a Dirichlet boundary condition:p for some given function g 1 . We may also consider the mixed boundary condition as in Section 2.3.3 of [12] (which corresponds to a Robin boundary condition for the pressure): for constants a, b ≥ 0 and a given function g 3 .

Absence of nutrients
Scaling mass, permeability, and mobility appropriately, by setting The existence of strong solutions in 2D and 3D have been studied in [37] for the case Γ = 0. For the case where Γ ≠ 0 is prescribed, existence of global weak solutions and unique local strong solutions in both 2D and 3D can be found in [32]. We also refer the reader to [9] for the study of weak solutions to a related system, denoted as the Cahn-Hilliard-Brinkman system, where an additional viscosity term η div D(v) is added to the left hand side of the velocity equation (2.43b) and the mass exchange Γ is set to zero.
We have elaborated on the physical motivations for the particular forms of Γ and S in Section 1. For the choice of N (ϕ, σ), if both χ ϕ and χ σ are positive constants, then for physically relevant parameter values, i.e., σ ≥ 0, and ϕ ∈ [−1, 1], Thus, this choice of the flux ∇N ,σ provides two transport mechanisms for the nutrient σ.
The first term χ σ ∇σ results in a diffusion process along negative gradients of σ, while the second term −χ ϕ ∇ϕ is a chemotactic term that drives the nutrient towards the tumour cell regions. In particular, in the tumour cell regions {ϕ = +1}, the nutrient will only experience diffusion, while in the healthy cell regions {ϕ = −1}, the nutrient will experience diffusion and active transport to the tumour. We point out that for this particular form of N ,σ , together with the zero Neumann boundary condition for ϕ, we have With these choices, (2.36) becomes We remark that (2.46) is similar to Eq. (68)-(73) of [14], the two-phase diffuse interface tumour model (Eq. 5.27) of [38], and model M 2 of [28]. The only difference between these three models and (2.46) is that the flux for the nutrient equation (2.46c) consists of an advection term and a Fickian diffusion term for [14], while in [38,28], the nutrient is in a quasi-steady state and the flux for the nutrient equation is a Fickian diffusive flux. We point out that in [14,28,38], h(ϕ) is replaced by ϕ in the definition of Γ and S. Since, in their notation, ϕ ∈ [0, 1] denotes the tumour volume fraction instead of the difference of volume fractions. Next, choosing N (ϕ, σ) as in (2.44b) above, and where P (ϕ) is a non-negative function, then (2.36) becomes This is similar to the model derived in [29], where the chemical potentials N ,σ and µ enter as source terms in (2.47a) and (2.47c). The specific form for Γ is motivated by linear phenomenological constitutive laws for chemical reactions. The non-negative function P (ϕ) takes on the form for positive constants δ and P 0 . Subsequently, if we choose This is the model studied in [22], for a more general function P (ϕ) than (2.48), while a viscosity regularised version of (2.49) (where there is an extra α∂ t µ term on the left hand side of (2.49a) and an extra α∂ t ϕ term on the right hand side of (2.49b) for a positive constant α) is studied in [13]. A formal asymptotic limit for the viscosity regularised version of (2.49) is derived in [30].

Sharp Interface Asymptotics
We consider Model (2.25) with the following choices and assumptions: Assumption 3.1.
• The mass exchange terms Γ i , i = 1, 2, and the nutrient consumption term S depend only on σ, µ, and ϕ, and not on any derivatives.
• The mobilities m(ϕ) and n(ϕ) are strictly positive and continuously differentiable.
• The potential Ψ is chosen to be either the smooth double-well potential Ψ(ϕ) = 1 4 (1 − ϕ 2 ) 2 or the double-obstacle potential With these choices, Model (2.25) becomes We point out that in the case of the double-obstacle potential, the "derivative" Ψ ′ is to be understood in the sense of subdifferentials, i.e., and (3.2d) will have to be formulated in terms of the following variational inequality: for all ψ ∈ K ∶= {η ∈ H 1 (Ω) ∶ η ≤ 1}. We perform a formal asymptotic analysis on Model (3.2) in the limit ε → 0. Details of the method can also be found in [2,7,8,23]. The following assumptions are considered: • We assume that for small ε, the domain Ω can be divided into two open subdomains Ω ± (ε), separated by an interface Σ(ε) that does not intersect with ∂Ω.
• We assume that there is a family (ϕ ε , v ε , p ε , µ ε , σ ε ) ε>0 of solutions to (3.2), which are sufficiently smooth and have an asymptotic expansion in ε in the bulk regions away from Σ(ε) (the outer expansion), and another expansion in the interfacial region close to Σ(ε) (the inner expansion).
• We assume that the zero level sets of ϕ ε converge to a limiting hypersurface Σ 0 moving with normal velocity V.
The idea of the method is to plug the outer and inner expansions in the model equations and solve them order by order, in addition we have to define a suitable region where these expansions should match up.
We will use the following notation: (3.2e) α O and (3.2e) α I denote the terms resulting from the order α outer and inner expansions of (3.2e), respectively.

Outer expansion
We assume that for f ε ∈ {ϕ ε , v ε , p ε , µ ε , σ ε }, the following outer expansions hold: The solutions to (3.5) corresponding to minima of Ψ are ϕ 0 = ±1, and thus, we can define the tumour tissues and the healthy tissues region by Then, thanks to ∇ϕ 0 = 0, we obtain from the equations to zeroth order: For the double-obstacle potential, we obtain from For this to hold for all ψ 0 ≤ 1, we require that ϕ 0 = ±1, and thus we can define Ω T and Ω H as before, and also recover (3.7), (3.8), (3.9), and (3.10).

Inner expansions and matching conditions
By assumption, Σ 0 is the limiting hypersurface of the zero level sets of ϕ ε . In order to study the limiting behaviour in these parts of Ω we introduce a new coordinate system. We introduce the signed distance function d(x) to Σ 0 , and set z = d ε as the rescaled distance variable, and use the convention that d(x) < 0 in Ω H , and d(x) > 0 in Ω T . Thus, the gradient ∇d points from Ω H to Ω T , and we may use ∇d on Σ 0 to denote the unit normal of Σ 0 , pointing from Ω H to Ω T .
Let g(t, s) denote a parametrization of Σ 0 by arc-length s, and let ν denote the unit normal of Σ 0 , pointing into the tumour region. Then, in a tubular neighbourhood of Σ 0 , for sufficiently smooth function f (x), we have In this new (t, s, z)-coordinate system, the following change of variables apply, compare [24]: , where V is the normal velocity of Σ 0 , ∇ Σ 0 g denotes the surface gradient of g on Σ 0 and h.o.t. denotes higher order terms with respect to ε. In particular, we have where κ = − div Σ 0 ν is the mean curvature of Σ 0 . Moreover, if v is a vector-valued function with V (t, s, z) = v(x) for x in a tubular neighbourhood of Σ 0 , then we obtain We denote the variables ϕ ε , v ε , p ε , µ ε , σ ε in the new coordinate system by Φ ε , V ε , P ε , Ξ ε , C ε , respectively. We further assume that they have the following inner expansions: The assumption that the zero level sets of ϕ ε converge to Σ 0 implies that Furthermore, we make the following assumption: For the double-obstacle potential, we assume that the inner variable Φ ε is monotone increasing with z and the interfacial layer has finite thickness of 2l, where the value of l will be specified later. For the double-well potential, we take l = ∞. Furthermore, we assume that Φ ε (t, s, z = +l) = +1, Φ ε (t, s, z = −l) = −1. (3.12) In order to match the inner expansions valid in the interfacial region to the outer expansions of Section 3.1 we employ the matching conditions, see [24]: where f ± 0 (t, x) ∶= lim δ↘0 f 0 (t, x ± δν) for x ∈ Σ 0 . Moreover, we use the following notation: Let δ > 0 and for x ∈ Σ 0 with x − δν ∈ Ω H and x + δν ∈ Ω T , we denote the jump of a quantity f across the interface by For convenience, we define the constant γ > 0 to be for the double-well potential, for the double-obstacle potential. (3.17)

Expansions to leading order
To leading order (3.2d) −1 I gives Using (3.11) we obtain that Φ 0 can be chosen to be independent of s and t, i.e., Φ 0 is only a function of z, and solves For the double-well potential, we have the unique solution (3.20) Furthermore, multiplying (3.19) by Φ ′ 0 (z), integrating and applying the matching conditions (3.13) and (3.14) to Φ 0 gives the so-called equipartition of energy: Similarly, for the double-obstacle potential, we obtain from (3.4) −1 I , For (3.22) to be satisfied, it suffices to consider Φ 0 as a function only in z which solves (3.24) so that l = π 2 for the double-obstacle potential, and we deduce from (3.12) that for the double-obstacle potential, Moreover, we obtain the equipartition of energy (3.21) via a similar argument to the double-well potential. Thanks to the equipartition of energy (3.21), and the definition of γ (3.17), we point out that For the rest of this section, we do not differentiate between the two cases of potentials, and use the notation that l = π 2 represents the case of the double-obstacle potential and l = ∞ represents the case of the double-well potential. Next, Integrating from −l to l with respect to z, and applying the matching condition (3.13) to V 0 yields We have from (3.2c) −2 I , Upon integrating and using the matching condition (3.14) applied to Ξ 0 , we obtain m(Φ 0 )∂ z Ξ 0 (t, s, z) = 0 for all z < l.
As n(Φ 0 ) > 0, we see that Integrating once more with respect to z from −l to l, and applying the matching condition (3.13) to C 0 and Φ 0 then gives Integrating and applying the matching condition (3.13) to P 0 and Ξ 0 leads to Thanks to (3.33), we see that Then, (3.36) becomes

Expansions to first order
For the double-well potential, to first order, we obtain from (3.2d) 0 I , We multiply (3.39) with Φ ′ 0 and integrate with respect to z from −∞ to ∞, which gives Applying integration by parts and the matching conditions (3.13) and (3.14) applied to Φ 0 , and using that Ψ ′ (±1) = 0, we see that and so the first two terms on the right hand side of (3.40) are zero. Then, using (3.26), (3.30), and (3.37), we obtain from (3.40), Moreover, together with (3.38), we obtain Meanwhile, for the double-obstacle potential, to first order, we obtain from (3.4) 0 I , Since Φ 0 (z) < 1 for z < π 2 , we can test with ψ = Φ 0 + λ with either non-positive or non-negative λ ∈ K, leading to the equality Multiplying with Φ ′ 0 and integrating with respect to z from − π 2 to π 2 , and applying matching conditions leads to Upon integrating by parts and using (3.23), the matching conditions (3.14) for Φ 0 , (3.15) for Φ 1 , and (3.25), we see that Then, using (3.26), and (3.37), we obtain from (3.44) the following solvability condition for Φ 1 : Lastly, thanks to (3.30), we obtain from (3.2c) −1 I and (3.2e) −1 I , respectively, and Thanks to (3.27), upon integrating from −l to l with respect to z, and applying the matching condition (3.15) to Ξ 1 , we obtain from (3.45) Similarly, thanks to ∇ϕ 0 = 0, upon integrating from −l to l with respect to z, and applying the matching condition (3.15) to C 1 and Φ 1 , we obtain from (3.46)

Sharp interface limit for Model (3.2)
In summary, we obtain the following sharp interface limit from Model (3.2): together with the free boundary conditions where γ is defined in (3.17). Note that we can write where σ 0 ∶= 1 2 (σ T 0 + σ H 0 ) denotes the average of the nutrient concentrations from both sides of Σ 0 . Thus, using (3.34), we can rewrite (3.41) and (3.50b) 2 as

Specific sharp interface models
In this section, we take Ψ as the double-well potential.

Sharp interface limit of the new active transport model
Choosing as before N (ϕ, σ) = χσ 2 σ 2 + χ ϕ σ(1 − ϕ) and for some positive constants P, A and C in Model (2.33), we obtain the Cahn-Hilliard-Darcy model (1.1) in Section 1. Then, the sharp interface limit of Model (2.33) with (3.53) is given by The active transport term n(ϕ)∇(χ ϕ ϕ) in the flux for the nutrient results in the jump term 2 χϕ χσ in (3.55d) 2 which is a new feature of the proposed model.

Linear constitutive laws for chemical reactions
Let us consider Model (2.47) with m(ϕ) = n(ϕ) = 1, and P (ϕ) defined as in (2.48). Then, we obtain that This was the setting introduced in [29]. Hence, we obtain from (2.47) the following sharp interface limit: in Ω H , (3.57b) We point out that the diffuse interface model studied in [30] takes a different mass transition term Γ and a different consumption term S. More precisely, the following choices are considered: (3.60) * In [30], the choice of P (ϕ) is actually This presents some difficulties in the analysis of the outer expansions, as P ′ (±1) ≠ 0 for the choice Ψ(ϕ) = (1 − ϕ 2 ) 2 . Thus, we do not recover (3.61a) and (3.61b). However, the formal analysis in [30] is different compared to what we present here, and it turns out that the analysis in [30] has to be modified and will only work for (3.60).
The term 1 ε Ψ(ϕ) acts as a regularisation on the Hausdorff measure restricted to the limiting hypersurface Σ 0 , and hence the reaction term σ − δµ will appear in the interfacial relations of µ and σ rather than in the bulk equations. More precisely, we obtain as a sharp interface limit of Model (2.36) with ρ = 1, m(ϕ) = n(ϕ) = 1, Γ and S as in (3.58) with P (ϕ) chosen as in (3.60) and N (ϕ, σ) = 1 2 σ 2 . This is similar to the sharp interface limit (Eq. (1.9)) of [30] with α = 0.
We briefly state the derivation of the sharp interface limit for Model (3.63). From From the leading order inner expansion (3.63b) −1 I , we obtain (3.18), and subsequently the profile (3.20) and the equipartition of energy (3.21). From (3.63c) −2 I we have Integrating and using the matching conditions (3.14), we obtain and upon matching, we obtain Integrating and using the matching condition (3.14) applied to Ξ 0 we deduce that Since Φ 0 (z) < 1 for z < ∞, we obtain that ∂ z Ξ 0 (t, s, z) = 0 for z < ∞.
To first order, we obtain from (3.63b) 0 I , Multiplying by Φ ′ 0 , using that Ξ 0 is independent of z, and applying integration by parts and matching conditions (3.13) and (3.14) to Φ 0 , we obtain, in the same spirit as (3.40), where we have used (3.64). Applying (3.51), and (3.65), we see that where we recall that σ 0 ∶= 1 2 (σ T 0 + σ H 0 ) is the average of the nutrient concentration at the interface. Meanwhile, thanks to (3.64) we obtain from (3.63c) −1 I , Integrating with respect to z from −∞ to ∞ and applying the matching condition (3.15) to C 1 and Φ 1 leads to Lastly, thanks to ∂ z Ξ 0 = 0, we obtain from (3.63a) −1 I , Integrating from −∞ to ∞ with respect to z, and applying the matching condition (3.15) to Ξ 1 gives Thus, the sharp interface limit of Model (3.63) is In addition, we can use (3.66c) 1 to rewrite (3.66c) 2 as Next, sending λ → 0 in (3.66) leads to and we define the bulk velocity and pressure via the relations: and from (3.66d) and (3.68), Thus, we obtain in Ω H , (3.73c) which coincides with the sharp interface model (Eq. (79)-(81), (83)-(86)) of [14]. We point out that the same sharp interface limit (3.73) can be recovered if we set λ = ε in (3.63). We introduce the parameter λ in (3.63) in order to study the effect of active transport on the linear stability of radial solutions to (3.66), see Section 4 below. Let us also remark that the mobility m(ϕ) = 1 2 (1 + ϕ) 2 is degenerate in the region {ϕ = −1}, and thus the bulk equation for µ H 0 in Ω H and the interfacial condition for ∇µ H 0 ⋅ν on Σ 0 remain undetermined in (3.66). Furthermore, if m(ϕ) is chosen to be degenerate in the bulk regions {ϕ = ±1}, then we obtain from the outer expansion (3.63a) 0 O the following equations 0 = Pσ 0 (ϕ 0 + 1) − A(ϕ 0 + 1) in Ω T ∪ Ω H .
In particular, we see that σ T 0 = A P is a constant in Ω T , which is inconsistent with (3.63c) 0 O . Hence, it is necessary that the mobility m(ϕ) is not degenerate in the tumour region {ϕ = 1}.

Linear Stability Analysis
Let us consider the sharp interface model (3.66). By sending the active transport parameter λ to zero, we recover the sharp interface model (Eq. (79)-(81), (83)-(86)) of [14]. In this section, we extend the linear stability analysis of [14,34] to include the effects of active transport. For the linear stability analysis of a one-phase model, we refer to [16,35].

Radial solutions
We now drop the index 0 in (3.66), and let Ω = B R (0) denote the d-dimensional ball, d = 2, 3, of radius R centered at the origin. We assume that the interface Σ is a (d − 1)sphere of radius q(t), partitioning the domain Ω into Ω T and Ω H as follows: The outer unit normal ν(p) at a point p ∈ Σ is given as while the normal velocity V is given as The mean curvature κ for a (d − 1)-sphere radius r 0 is given by where d denotes the dimension. Then, for radially symmetric solutions ϕ( x ) = ϕ(r), µ( x ) = µ(r), σ( x ) = σ(r), (3.66) becomes We complete (4.4) with the following boundary conditions: where σ ∞ denotes the concentration of a nutrient supply from the boundary ∂Ω. Upon solving the differential equations and applying the interface and boundary conditions, we arrive at the following radial solutions: where, for α ∈ R, I α (x) denote the modified Bessel functions of the first kind : (4.7) Here, Γ(⋅) denotes the Gamma function. Together with the modified Bessel functions of the second kind, , the pairs {I α , K α } are the two linearly independent solutions to the modified Bessel's equation: Moreover, for the case α = 0, the modified Bessel functions I 0 (x), K 0 (x) satisfy the following properties Furthermore, the coefficients in (4.6) are given as a 2 (t) = q(t)ΛI 1 (Λq(t))(σ ∞ + 2λ) DI 0 (Λq(t)) − Λq(t) log(q(t) R)I 1 (Λq(t)) , , (4.10g) and the differential equation satisfied by q(t) is (4.11) We point out that, thanks to the boundary condition σ T (r = 0, t) < ∞, the solution σ T does not contain any terms involving K 0 (Λr) (in d = 2) and cosh(Λr) r (in d = 3).

Perturbation of radial solutions
We now consider a perturbation of a radially symmetric tumour, whose radius w is given by where q(t) is the radius of the unperturbed interface, δ(t) is a dimensionless perturbation size, Y l,m is a spherical harmonic with l and θ denoting the polar wavenumber and angle, and m and φ denoting the azimuthal wavenumber and angle, respectively. We will denote the radial solutions in (4.6) by σ * H , σ * T , and µ * T , and consider where we assume that (σ H , σ T , µ T ) solve (3.66). Therefore, we get Here, we used the more convenient form (3.67) of (3.66c) 2 . Next, we linearise (4.14) about the original interface r = q to derive the equations satisfied by U , V , W and δ. We introduce the Laplace-Beltrami operator on the (d − 1)-sphere, for d = 2, 3: (4.15) so that the Laplace operator can be decomposed into Moreover, the function Z(θ, φ) defined in (4.12) satisfies (4.17) From the bulk equation (4.14a) we obtain ∆µ * T + δ∆(W Z) = A − Pσ * T − δPV Z, and so, using that ∆µ * T = A − Pσ * T , we deduce that For the interface conditions, we employ Taylor's expansion and neglect terms of order O(δ 2 ). For instance, from (4.14d), we see that Then, by (3.66c) 1 , (3.66d) 1 , and (4.12), we obtain We use the following expansion for the mean curvature (compare with Eq. (4.12) of [19], page 647 of [17] and page 12 of [21], where instead of (4.3), the mean curvature of a (d − 1)-sphere radius r 0 is defined to be 1 r 0 ): so that the linearisation of (4.14e) about r = q is Finally, by the relation ∇f ( x ) ⋅ ν = f ′ (r) for x ∈ Σ, we have that and so we obtain the following system for the perturbations U, V, W and δ from linearising (4.14) about the unperturbed interface r = q: We complete (4.18) with the following boundary conditions:

Solutions to the perturbed system
Recalling the definition of ζ l,d in (4.17), we see that the general solution for (4.18c) is (4.20) We observe that the ODE (4.18b) in d = 2 is a scaled modified Bessel's equation (see (4.8)), while (4.18b) in d = 3 is a scaled modified spherical Bessel's equation. Due to the boundary condition (4.19), we see that the general solution to (4.18b) is given by where I l (x) is the modified Bessel function of the first kind, defined in (4.7), while i l (x) is the modified spherical Bessel function of first kind that satisfies Again, due to the boundary condition (4.19) 2 , V (r, t) does not contain any terms involving the modified spherical Bessel function of the second kind. For (4.18a), we see that W is a sum of the solution to the homogeneous equation (4.18c) and the particular solution P C V . Due to the boundary condition (4.19) for W , we find that the general solution to (4.18a) is then from (4.18d), (4.18e), and (4.18f) we obtain the following relations: where we have used that (4.24) and by (4.4d) 1 , Also, from (4.4d) 1 , we observe that the following relations hold Together with the relation and the relations (4.23), we can simplify (4.18g) in order to obtain the following differential equation for the perturbation size δ: and (4.28) Consequently, using (4.11), and (4.26), we obtain the following differential equations for the shape perturbation δ q : (4.29) and the relations b 2 I 0 (q) = DI 0 (q) DI 0 (q) − q log(q R)I 1 (q) = I 0 (q) I 1 (q) we find that (4.32) Substituting F 0 = 0 and λ = 0 in (4.32) and (4.29), we obtain the differential equation for the shape perturbation as derived in Eq. (89) of [14] with the notationG −1 ∶= 1 2 βγ. Next, we find, for given P, D, χ ϕ , and β, a critical value A c such that d dt δ q = 0, i.e., the shape perturbation ( δ q ) is a constant. This critical value A c is given by Numerically, we find that C 2 , C 3 , X and Y are positive for q ∈ (0, 13]. Moreover, We note that A is the apoptosis parameter and A c divides the phase portrait into regions of stable growth for low apoptosis (the region A < A c ) and regions of unstable growth for high apoptosis (the region A > A c ) for a given mode l. Thus, from (4.33) and (4.34), we observed the following: 1. In the absence of chemotaxis, χ ϕ = 0, increasing λ will increase the value of A c . From  Figures 1(a) and 1(b), the curves are pushed upwards, and so the region of stable growth for low apoptosis is enlarged. In particular, active transport has a stabilising effect on the perturbations in the absence of chemotaxis.
2. In dimension d = 2, while χ ϕ < 0.2, active transport has a stabilising effect on the perturbations. When χ ϕ > 0.2, the perturbations are now amplified by the presence of active transport. In Figure 1(c), we see that, as λ increases, the curves are pushed up for χ ϕ = 0.1, while the curves are pulled down for χ ϕ = 0.3. Similarly, in dimension d = 3, we find that 0.25 q − 0.15Y (q) 1 q − Y (q) ∈ (0.400, 1.459) for q ∈ [0.01, 13], and from Figure 1(d), we see that, as λ increases, the curves are pushed up for χ ϕ = 0.3, while the curves are pulled down for χ ϕ = 1.7.

Numerical Computations
In this section we first derive a finite element approximation of (3.63) and then we display some numerical results obtained using this approximation. We concentrate on (3.63), however approximations of other variations of the model follow in a natural way. In the approximation we take Ψ(ϕ) to be the double obstacle potential given in (3.1). This choice of potential leads to (3.63b) taking the form of a variational inequality (3.4).

Finite element approximation
Let T be a regular triangulation of Ω into disjoint open simplices, associated with T is the piecewise linear finite element space where we denote by P 1 (T ) the set of all affine linear functions on T . We now introduce a finite element approximation of (3.63) in which we have taken homogeneous Neumann boundary conditions for ϕ and µ, and the Dirichlet boundary condition σ = σ B ∈ R on ∂Ω: Find where m(ϕ) = 1 2 (1 + ϕ) 2 , τ denotes the time step, (η 1 , η 2 ) denotes the L 2 inner product and (η 1 , η 2 ) h ∶= ∫ Ω π h (η 1 (x)η 2 (x)) dx where on each triangle π h is taken to be an affine interpolation of the values of η 1 η 2 at the nodes of the triangle.
We note that since the interfacial thickness is proportional to ε in order to resolve the interfacial layer we need to choose h ≪ ε, see [18] for details. Away from the interface h can be chosen larger and hence adaptivity in space can heavily speed up computations. In fact we use the finite element toolbox Alberta 2.0, see [42], for adaptivity and we implemented the same mesh refinement strategy as in [5], i.e., a fine mesh is constructed where ϕ n−1 h < 1 with a coarser mesh present in the bulk regions ϕ n−1 h = 1. We begin our numerical results by following the authors in [30] in comparing solutions obtained from a simplified form of the diffuse interface model with exact solutions to a sharp interface limit model.

Comparison with a sharp interface limit solution
In Figures 2 and 3 we display results obtained from the growing circle tumour test case introduced in Section 4.2 of [30]. To this end we consider the simplified model on a circular domain Ω with radius R: Here ϕ and µ satisfy homogeneous Neumann boundary conditions, and σ satisfies the Dirichlet boundary condition σ = σ R ∈ R on ∂Ω. We take the radially symmetric case of an initial circular tumour with initial radius 0.25. From [30] we have that the solution to the sharp interface limit of (5.2) is given by , with µ being constant and ρ(t), which is the radius of the tumour, is determined by numerically solving the ODE ρ ′ (t) = √ 2σ(ρ(t), t) with initial condition ρ(0) = 0.25.
We set R = 10 and σ R = 2, however for the diffuse interface computations we did not solve the problem in the whole of Ω instead we solved it on a circular domain with radius 2 with the time dependent Dirichlet boundary condition σ(x, t) = σ D ( x , t) computed from (5.3) with r = 1. We set τ = 1.0e −4 , the minimal diameter of an element h min = 7.8125⋅10 −3 and the maximal diameter h max = 3.125 ⋅ 10 −2 .
In Figure 2 we display the diffuse interface solutions ϕ and σ at t = 0, 0.2, 0.4 obtained with ε = 0.05. In the plots of ϕ we include the sharp interface limit solution of the tumour position. In Figure 3 we examine the convergence of the diffuse interface solution to the sharp interface limit solution as ε tends to zero. In Figure 3(a) we plot the radius of the growing tumour for the diffuse interface model with ε = 0.1, 0.075, 0.05 together with the sharp interface limit solution ρ(t). In Figure 3(b) we plot the solution σ of the diffuse interface model with ε = 0.1, 0.075, 0.05 together with the sharp interface limit solution σ at t = 0.1. From this figure we see that as ε decreases the diffuse interface solution converges to the sharp interface limit solution.

Influence of the proliferation rate P
In Figures 4 and 5 we investigate the influence of P. We set χ ϕ = 5 and λ = 0.03. In Figure  4 we set P = 0.5 while in Figure 5 we set P = 0.1, and in both sets of figures we display ϕ (top row) and σ (bottom row) at times t = 5, 10, 13. From this figure we see taking P = 0.5 gives rise to fingers that are thicker than the ones resulting from P = 0.1.  Influence of the chemotaxis parameter χ ϕ In Figures 6 and 7 we investigate the influence of χ ϕ . We set P = 0.1 and λ = 0. In Figure  6 we set χ ϕ = 5 while in Figure 7 we set χ ϕ = 10, and in both sets of figures we display ϕ (top row) and σ (bottom row). The results for χ ϕ = 5 are displayed at times t = 5, 10, 20, while the results for χ ϕ = 10 are displayed at times t = 2.5, 5, 10. From these figures we see that, akin to the results in [14], for both values of χ ϕ after some time fingers develop, and thereby increasing the surface area of the tumour to allow for better access to the nutrient. For the larger value of χ ϕ the formation and evolution of the fingers is quicker and the fingers are slimmer. Influence of the active transport parameter λ In Figures 8 -10 we investigate the influence of λ. We set P = 0.1 and χ ϕ = 5. In Figure  8 we show ϕ (top row) and σ (bottom row) at t = 4, with λ = 0 (left), λ = 0.07 (centre) and λ = 0.09 (right). From this figure we see that when λ = 0 the variation of σ across the interfacial region is smooth while taking λ > 0 leads to a drastic change in σ. This change in σ can be seen better in Figure 9 where we show plots of ϕ and σ along a line that spans the interfacial region. The scales for σ and ϕ are shown on the left and right axes respectively.
Here we see that the change in σ across the interfacial region is more pronounced for larger values of λ. In Figure 10 we display the influence of ε on the change in σ across the interfacial region, we set λ = 0.07 and plot σ along a line that spans the interfacial region for ε = 0.04, 0.02, 0.01. From this figure we see the convergence of σ as ε decreases. In Figure 10 the jump in σ across the interfacial region for ε = 0.01 is 0.1327 ≈ 2λ which is consistent with the formal asymptotic analysis, recall (3.65). where we recall that α ∶= 1

Numerical computations with Darcy flow
, Γ = 1 2 (Pσ − A)(ϕ + 1), and D is defined in (3.62). As additional boundary condition we prescribe p = 0 on ∂Ω, while we take homogeneous Neumann boundary conditions for ϕ and µ, and the Dirichlet boundary condition σ = σ B ∈ R on ∂Ω. Recalling the finite element spaces K h , S h , S B h and S 0 h defined at the start of Section 5, for the double-obstacle potential (3.1), we propose the following scheme for the above system: Find ϕ n h ∈ K h , µ n h ∈ S h , σ n h ∈ S B h , p n h ∈ S 0 h such that for all η h ∈ S h , ζ h ∈ K h and χ h ∈ S 0 h , As initial condition for p and µ, we always choose p 0 h = 0 and µ 0 h = 0. We perform three different numerical simulations in which we vary the tumour and healthy cell densities. The three cases are given as follows: