Turing instability analysis of a singular cross-diffusion problem

The population model of Busenberg and Travis is a paradigmatic model in ecology and tumour modelling due to its ability to capture interesting phenomena like the segregation of populations. Its singular mathematical structure enforces the consideration of regularized problems to deduce properties as fundamental as the existence of solutions. In this article we perform a weakly nonlinear stability analisys of a general class of regularized problems to study the convergence of the instability modes in the limit of the regularization parameter. We demonstrate with some specific examples that the pattern formation observed in the regularized problems, with unbounded wave numbers, is not present in the limit problem due to the amplitude decay of the oscillations. We also check the results of the stability analysis with direct finite element simulations of the problem.


Introduction
Busenberg and Travis [4] introduced a class of singular cross-diffusion problems under the assumption that the spatial relocation of each species is caused by a diffusion flow which depends on the densities of all the involved species. In the case of two species, if u 1 and u 2 denote their densities, the flow, in its simplest form, may be assumed to be determined by the total population u 1 + u 2 , and thus the conservation laws for both species lead to the system ∂ t u 1 − div u 1 (∇u 1 + ∇u 2 ) = f 1 (u 1 , u 2 ), (1.1) ∂ t u 2 − div u 2 (∇u 1 + ∇u 2 ) = f 2 (u 1 , u 2 ). (1.2) Here, the functions f 1 and f 2 capture some ecological features of the populations, such as growth, competition, etc., and the equations are complemented with nonnegative initial data and non-flow boundary conditions. System (1.1)-(1.2) is called a cross-diffusion system because the flow of each species depend upon the densities of the other species. We call it singular because the resulting diffusion matrix is singular. Indeed, when rewriting (1.1)-(1.2) in matrix form, for u = (u 1 , u 2 ), we obtain the equation where the divergence is applied by rows. The full and singular structure of A introduces serious difficulties in the mathematical analysis of the problem, as we shall comment later.
In his seminal paper [14], Turing introduced a mechanism explaining how spatially uniform equilibriums may evolve, small perturbations mediating, into stable equilibriums with non-trivial spatial structure. He considered a system of the type with σ > 0, and proved that when σ is small or large enough then the stable equilibrium of the dynamical system is not stable for the diffusion system (1.3)-(1.4), being replaced by a non-uniform equilibrium with spatial structure or patterns. This mechanism is known as Turing instability or Turing bifurcation.
In this article we study Turing instability for the cross-diffusion singular system (1.1)-(1.2). We already know that some cross-diffusion systems, such as the paradigmatic SKT model introduced by Shigesada, Kawasaki and Teramoto [13], exhibit Turing instability when cross-diffusion coefficients are large in comparison with self-diffusion coefficients, see e.g. [10,11]. However, the singularity of the diffusion matrix of the system (1.1)-(1.2), not present in the SKT model, introduces important mathematical difficulties to the analysis of this system.
Regarding the existence of solutions of (1.1)-(1.2), it has been proved only in some special situations: for a bounded spatial domain Ω ⊂ R (Bertsch et al. [2]) and for Ω = R n (Bertsch et al. [3]). In their proofs, the following observation is crucial: adding the two equations of (1.1)-(1.2) we deduce that if a solution of this system does exist then the total population, u = u 1 + u 2 , satisfies the porous medium type equation for which the theory of existence and uniqueness of solutions is well established. In particular, if the initial data of the total population is bounded away from zero and if f is regular enough with f (0) ≥ 0, it is known that the solution of (1.7) remains positive and smooth for all time. This allows to introduce the change of unknowns w 1 = u 1 /u into the original problem (1.1)-(1.2) to deduce the equivalent formulation for certain well-behaved functions F 1 and F 2 . This change reveals the parabolichyperbolic structure of the problem, which is handled by Bertsch et al. [2,3] by the vanishing viscosity method, adding the regularization term −δ∆w 1 to the left hand side of (1.9), and justifying the passing to the limit δ → 0 in the corresponding sequence of solutions.
In [9] we followed a different approach to prove the existence of solutions of the original system (1.1)-(1.2) for a bounded domain Ω ⊂ R. We introduced a parabolic regularization of the original system by adding a cross-diffusion perturbation term while keeping the porous medium type equation satisfied by u. More concretely, we considered the system 11) and then used previous results for cross-diffusion systems [7,6,12] to establish the existence of solutions of the approximated problems. Then, estimates of the bounded variation (BV) similar to those in [2] allowed us to prove the convergence of the sequence (u 2 ) to a solution of the original problem. Let us mention that the system (1.1)-(1.2) is a limit case of a general type of problems with diffusion matrix given by for which, if a ii > 0, for i = 1, 2, and a 11 a 22 > a 12 a 21 , then the existence of solutions is ensured for any spatial dimension of Ω, see [8]. In addition, it has been recently shown that this kind of systems, when set in the whole space Ω = R n , may be rigorously deduced as mean field limits [5]. Concerning Turing instability, since the diffusion matrix corresponding to the system (1.1)-(1.2), A(u), is singular, the linearization of this system at an equilibrium of the dynamical system (1.5)-(1.6) does not provide any information on the behavior of the equilibrium in the spatial dependent case. Thus, our approach to Turing instability for the system (1.1)-(1.2) relies on the study of this property for regularized problems like (1.10)-(1.11) and their limit behavior.
For the regularized problems, we prove that linear instability is always present for δ small enough. Interestingly, the linear analysis also establishes that the main instability wave number is unbounded as δ → 0.
For understanding this convergence of a increasingly oscillating sequence of functions to a BV function (the solution of (1.1)-(1.2) ensured in [2,9]), we perform a weakly nonlinear analysis (WNA) which allows to gain insight into the behavior of the amplitude of the main instability mode as δ → 0. As expected, we find that the amplitude of the main instability mode vanishes in the limit δ → 0 and is, therefore, coherent with the BV convergence. In addition, this result also suggests that the uniform equilibrium is stable for the original problem. We, finally, illustrate our analytical results with some numerical examples.

Main results
For simplicity, we study Turing instability for the one-dimensional spatial setting which has also the advantage of a well established existence theory for the case of a bounded domain [2,9]. Redefining the functions f 1 and f 2 , we can fix without loss of generality Ω = (0, π) and then rewrite problem (1.1)-(1.2) together with the auxiliary conditions as u 1 (0, ·) = u 10 , u 2 (0, ·) = u 20 in Ω, (2.4) where Q T = (0, T ) × Ω and the initial data u 10 , u 20 are non-negative functions. We assume a competitive Lotka-Volterra form for the reaction term, this is, , for i = 1, 2, and for some non-negative parameters α i , β ij , for i, j = 1, 2.
To deal with several types of regularized problems we introduce, for positive δ and b, the uniformly parabolic cross-diffusion system u 1 (0, ·) = u 10 , u 2 (0, ·) = u 20 in Ω, (2.8) where the diffusion matrix D δ (u) = (d δ ij (u)) and the Lotka-Volterra function b (u) = (f b 1 (u), f b 2 (u)) satisfy the following assumptions. Assumption (H1): (1) D δ (u) is linear in u and affine in δ, so that it allows the decomposition for some matrices D δi for i = 1, 2, being the coefficients of D δ (u) given by ij u 2 δ, for some non-negative constants d mn ij , for i, j, m = 1, 2 and n = 0, 1. (2) We assume that d δ ii (u) > 0 for i = 1, 2, and that det(D δ (u)) is an increasing function with respect to δ satisfying det(D δ (u)) > 0 if δ > 0 and u ∈ R 2 + . (2.10) Observe that (2.10) ensures the existence of a stable coexistence equilibrium for the dynamical system (1.5)-(1.6), given by There are two examples of D δ (u) in which we are specially interested. The first is chosen because of its simplicity for the calculations, with det(D δ (u)) = δ(2 + δ)u 1 u 2 . According to [8], the second hypothesis of (H1) guarantees the well-posedness of the problem (2.5)-(2.8) corresponding to this diffusion matrix. The second example corresponds to the approximation used in [9] for proving the existence of BV solutions of the original problem (2.1)-(2.4): Our first result provides conditions under which linear instability arises. The following notation is used, (2.14) then there exists δ c > 0 such that if δ < δ c then u * is a linearly unstable equilibrium for problem (2.5)-(2.8). In such situation, the wave number of the main instability mode tends to infinity as δ → 0.
and introduces a further restriction on the competence coefficients. Roughly speaking, for B b to fulfill (2.10) and (2.16), its elements must be such that intrapopulation joint competence is larger than inter-population joint competence (condition (2.10)) and one of the inter-population competence coefficients is large in comparison with the others (condition (2.16)). A numeric example we shall work with along the article is Assuming the forms of D δ (u * ) given in Experiments 1 and 2, see (2.12) and (2.13), we have that the conditions (2.10) and (2.16) on B b are satisfied if δ < b/4 (Experiment 1) or δ < bu * 1 u * 2 /(u * 1 + u * 2 ) 2 (Experiment 2), which are always true for the case of interest: δ close to zero.
Our second result allows to estimate not only the instability wave numbers provided by the linear analysis but also the amplitude corresponding to these modes. The approximation of the steady state solution is obtained by using a weakly nonlinear analysis (WNA) based on the method of multiple scales.
where k c ∈ Z is the critical wave number corresponding to δ c , A ∞ is a positive constant and ρ, v 20 and v 22 are constant vectors.
Our third result focuses on the limit behavior of the critical parameters and the amplitude when δ → 0, this is, when the solutions of the regularized problems converge to the solution of the original singular problem. For the sake of simplicity, we limit our study to the following example: 20) whose solutions we approximate by the two-parameter family of solutions of On one hand, Theorem 2.1 ensures the existence of δ c > 0 such that, for any b ≥ 0, the equilibrium u * = 1 1−b (1 − 2b, 2) of (2.21)-(2.22) becomes unstable for δ < δ c , with an associated critical wave number such that k c → ∞ as δ → 0.
On the other hand, for δ < b/4 and b → 0, the sequence of solutions of (2.21)-(2.22) converges to a solution of (2.19)-(2.20) in the space BV (0, T, L ∞ (Ω)) ∪ L ∞ (0, T ; BV (Ω)). Therefore, for the approximation (2.18) provided by the weakly nonlinear analysis to remain valid for all δ > 0, the corresponding amplitude A ∞ must vanish in the limit δ → 0, making in this way compatible the increase of oscillations with its BV regularity.
and the amplitude provided by the weakly nonlinear analysis satisfies In particular, the weakly nonlinear approximation v given by (2.18) satisfies v → u * uniformly in Ω as b → 0.

Numerical experiments
To analyze the quality of the approximation provided by the WNA, and the properties stated in Theorems 2.1-2.3, we compare it to a numerical approximation of the evolution problem computed with the finite element method (FEM).
For the FEM approximation, we used the open source software deal.II [1] to implement a time semi-implicit scheme with a spatial linear-wise finite element discretization. For the time discretization, we take in the experiments a uniform time partition of time step τ = 0.01. For the spatial discretization, we take a uniform partition of the interval Ω = (0, π) with spatial step depending on the predicted wave number of the pattern, see Table 1.
Since (3.1)-(3.2) is a nonlinear algebraic problem, we use a fixed point argument to approximate its solution, (u n 1 , u n 2 ), at each time slice t = t n , from the previous approximation (u n−1 . Then, for k ≥ 1 the linear problem to solve is: Find (u n,k 1 , u n,k 2 ) such that for for all χ ∈ S h 1 τ u n,k 1 − u n−1 We use the fixed point stopping criteria for values of tol F P chosen empirically, and set (u n 1 , u n 2 ) = (u n,k 1 , u n,k 2 ). Finally, we integrate in time until a numerical stationary solution, (u S 1 , u S 2 ), is achieved. This is determined by where tol S is chosen empirically too. In the following experiments we always fix tol F P = 1.e − 07 and tol S = 1.e − 12. Experiment 1. We investigate the behavior of the instabilities arising in the solutions of the approximated problems (2.5)-(2.8) when δ → 0. Our main objective is checking if the predictions of the weakly nonlinear analysis stated in Theorem 2.3 are captured by the FEM approximation too. Thus, we use the diffusion matrix D δ (u) and the competence parameters B b given by (2.12) and (2.17), respectively.
We run three simulations according to the choice of b, see Table 1, and fix δ = 0.95δ c (b) in all of them, so that u * is unstable and pattern formation follows. In Figure 1 we show the typical onset and transmission of disturbances found in all the experiments. In this figure and in the following we plot only the first component of the solution, being the behavior of the second component similar. After a fast decay of the initial data towards the unstable equilibrium, a perturbation with the wave number predicted by the linear analysis grows from one side of the boundary to the rest of the domain until reaching the steady state, see Figure 2. In the latter figure, we may check the good accordance between the FEM and the WNA approximations which, in numeric figures, have a relative difference of the order 10 −5 .
In Figure 3 we show three interesting behaviors of solutions when δ → 0. In the left panel, the shrinking amplitude of the stationary patterns while the wave number increases. The equilibrium has been subtracted from the solution to center the pattern in y = 0. The center panel shows the time evolution of the amplitude (log scale) as given by the exact solution of the Stuart-Landau equation (4.17). We readily see that the stabilization time is a decreasing function of δ. This fact together with the increment of the wave number when δ → 0 results in very high execution times, see Table 1. Finally, the third panel shows how the variation of the numerical stationary solution is an increasing function of δ and tends to zero as δ → 0, in agreement with the regularity of solutions stated by the theoretical results. Experiment 2. We repeated Experiment 1 replacing the diffusion matrix D δ (u) by that defined in (2.13). In Table 2 we show the relative differences in L p , given by in terms of the critical bifurcation parameter, δ c , the stationary solution of the FEM approximation, u(T, ·), the WNA approximation, v, and the pattern amplitude, A ∞ , corresponding to both approximations of the original diffusion matrix. We see that although the critical bifurcation parameter is clearly affected by the approximation scheme, the FEM and WNA approximations provided by both schemes are in a very good agreement, as well as the amplitudes of the instability patterns, suggesting that in the limit δ → 0 both sequences of approximations converge to the same limit.

Proofs
We use the decomposition of the nonlinear problem (2.5)-(2.8) in terms of its linear and nonlinear parts. Let v = u − u * , where u is a solution of (2.5)-(2.8). Then, v satisfies where we split the reaction-diffusion terms into their linear parts with K given by (2.14), and their nonlinear parts Proof of Theorem 2.1. We study the linearization of (4.1), this is, the equation satisfying Neumann homogeneous boundary conditions and with initial data w 0 = u 0 − u * . This linear problem is well-posed due to the second assumption of (H1). The type of boundary conditions lead to seek for solutions of the form w = e λt cos(kx)w, with k = 1, 2, . . ., where w is a constant vector. Replacing w in (4.3) we obtain the matrix eigenvalue problem Since, by hypotheses, tr(A k ) = tr(K) − k 2 tr(D δ (u * )) < 0 for all k = 0, 1, . . . ; an eigenvalue with positive real part (instability) may exist only if det(A k ) is negative for some wave number k. We introduce the notation h(k 2 ) = det(A k ): . The minimum of the convex parabola h is attained at requiring q δ (u * ) < 0, which is true in view of (2.15). A necessary condition for .
In this expression, det(K) is a positive constant and q δ (u * ) 2 > 0 for all δ ≥ 0. Thus, since q δ (u * ) and D δ (u * ) are monotone with respect to δ and det(D δ (u * )) → 0 as δ → 0, we deduce the existence of an uniqueδ c > 0 such that h(k 2 m (δ c )) = 0. Therefore, for δ <δ c we have h(k 2 (δ)) < 0 if k 2 (δ) ∈ (k 2 − (δ), k 2 + (δ)), where Because of the boundary conditions, the onset of instabilities only occurs when one of the extreme values of the interval (k − (δ), k + (δ)) is an integer number. Since k + (δ) → ∞ as δ → 0, this will certainly hold for δ small enough. We define the critical bifurcation parameter, δ c , as such number, and the critical wave number, k c ∈ Z, as the corresponding root of h(k 2 ). Finally, the last assertion of the theorem is a consequence of the infinite limit of k + (δ) as δ → 0.
Proof of Theorem 2.2. We retake the whole nonlinear equation (4.1) for v = u−u * . The idea of the weakly nonlinear analysis is to look for an approximation of v for a value of δ near the critical bifurcation parameter δ c . This approximation is defined as an expansion in terms of the small parameter 2 = (δ c − δ)/δ c , for δ < δ c . We consider the expansions and then introduce these expressions in equation (4.1) and collect the resulting equations in terms of powers of . Since this procedure is standard, we give the results and omit intermediate calculations for the sake of brevity. We obtain Order : Order 2 : Order 3 : (4.6) Here, D 1 (u * ) is given by (2.9), , where we introduced the notation d δm ij are the elements of the matrices D δ1 , D δ2 introduced in the first assumption of (H1). Observe also that (4.2) may be written as We now compute the solutions corresponding to each order in the expansion.
Order : The solution of (4.4) is given by , where A is the amplitude of the pattern, unknown at the moment. Observe that ker(A δc kc ) is a one-dimensional subspace, implying that the vector ρ is defined up to a multiplicative constant. We shall fix this constant later.
Order 2 : We start by expressing F in terms of A and ρ. We have Using standard trigonometric identities, we obtain ). Gathering the above expressions, we obtain . By Fredholm's alternative, (4.5) admits a solution if and only if F, ψ L 2 = 0, where ·, · L 2 denotes the scalar product in L 2 (0, π), and ψ ∈ ker((L δc ) * ) is of the form Observe that η, for similar reasons than ρ, is defined up to a multiplicative constant.
We fix η at the end of this proof, and also show that ρ, η = 0. The compatibility condition implies Since the solution to this equation is an exponential function, we do not obtain from it any useful indication on the asymptotic behavior of the pattern amplitude. Therefore, to suppress the secular terms appearing in F, we impose t 1 ≡ 0 and δ 1 ≡ 0. (4.8) In particular, this implies A ≡ A(t 2 ). Assuming these restrictions, the Fredholm's alternative is satisfied, and motivated by the functional form of F, we seek for a solution of (4.5) of the form where v 2j are constant vectors. The linear operator L δc may be decomposed as Then, L δc v 2 = F if the vectors v 2j are the solutions of the linear systems Order 3 : We have to solve L δc v 3 = G, where, taking into account (4.8), Replacing the solutions obtained for the orders and 2 , i.e. v 1 = A(t 2 )ρ cos(k c x) and v 2 = A(t 2 ) 2 (v 20 + v 22 cos(2k c x)) in G yields The solvability condition for problem (4.6) is G, ψ L 2 = 0, with ψ = η cos(k c x) given by (4.7). This condition leads to the differential equation Thus, we deduce the cubic Stuart-Landau equation for the amplitude We, finally, fix the vectors ρ ∈ ker(K − k 2 c D δc (u * )) , and η ∈ ker((K − k 2 c D δc (u * )) * ). Since all the elements of both matrices are negative, we may set ρ = (1, M ) t and η = (1, M * ) t for some M, M * < 0, implying ρ, η > 0. Thus, the asymptotic behavior of the solution to (4.10) is fully determined by the signs of the numerators in (4.11). When σ and are positive, the amplitude stabilizes to a positive value, this is, A(t 2 ) → A ∞ := σ/ as t 2 → ∞. Therefore, in this case, the corresponding An example of this situation is studied in Theorem 2.3.
Proof of Theorem 2.3. Our objective is to compute the coefficients of the Stuart-Landau equation (4.10). Specifically, we are interested in the ratio Determination of G 1 , η . For the given data, we obtain q δ (u * ) = −u * 1 u * 2 b 2 − 2δ , which is negative if b > 4δ. The corresponding roots of h(k 2 m ) are positive and, therefore, we take δ c = δ − , so that for any δ < δ c we have h(k 2 m ) < 0. The corresponding critical wave number is the minimum of h(k 2 ), given by .