Numerical continuation for a fast-reaction system and its cross-diffusion limit

In this paper we investigate the bifurcation structure of the triangular SKT model in the weak competition regime and of the corresponding fast-reaction system in 1D and 2D domains via numerical continuation methods. We show that the software pde2path can be adapted to treat cross-diffusion systems, reproducing the already computed bifurcation diagrams on 1D domains. We show the convergence of the bifurcation structure obtained selecting the growth rate as bifurcation parameter. Then, we compute the bifurcation diagram on a 2D rectangular domain providing the shape of the solutions along the branches and linking the results with the Turing instability analysis. In 1D and 2D, we pay particular attention to the fast-reaction limit by always computing sequences of bifurcation diagrams as the time-scale separation parameter tends to zero. We show that the bifurcation diagram undergoes major deformations once the fast-reaction systems limits onto the cross-diffusion singular limit. Furthermore, we find evidence for time-periodic solutions by detecting Hopf bifurcations, we characterize several regions of multi-stability, and improve our understanding of the shape of patterns in 2D for the SKT model.


Introduction
Systems with multiple time-scales appear in a wide variety of mathematical areas but also in many applications [38]. The basic class of multiple time-scale systems has two time-scales, which are so-called fast-slow systems given by du dt where u, v are the unknowns and ε > 0 is a small parameter so that u is fast while v is slow. There is quite substantial theory for the case of fast-slow ordinary differential equations (ODEs). Although there are many links of the ODE theory to multiscale partial differential equations (PDEs) [38,39], a lot less is known about fast-slow PDEs. One important sub-class are fast-reaction PDEs given by posed on a domain R + × Ω, where ∆ is the usual spatial Laplacian with respect to x ∈ Ω, t ∈ [0, T ) for some T > 0, u = u(t, x), and v = v(t, x). In this work, one of the two key equations we study is a particular version of (1.2) arising in mathematical ecology [30] ∂ t u 1 = d 1 ∆u 1 + ( where the quantities u 1,2 (t, x), v(t, x) ≥ 0 represent the population densities of two species at time t and position x, confined and competing for resources on a bounded domain Ω ⊂ R N . The positive coefficients d i , r i , a i , b i (i = 1, 2) describe the diffusion, the intrinsic growth, the intra-specific competition and the inter-specific competition rates respectively. The population u is split into two states, distinguishing between less active and active individuals, respectively denoted by u 1 and u 2 (thus u = u 1 + u 2 ). We assume throughout that 0 ≤ u 2 (t, x) ≤ M on [0, T ] × Ω for a fixed constant M > 0. The coefficient d 12 > 0 stands for competition pressure between the two sub-classes. Each individual in the sub-classes converts its state into the other one depending on the spatial distribution of the competitor v.
Since ε > 0 is assumed to be small, it is natural to consider the fast-reaction limit ε → 0. This limit aims to model the dynamics of the fast variables as instantaneous and reduce the analysis to the slow time-scale dynamics. Yet, one easily sees that this limit is very singular. To motivate this observation, suppose we discard the diffusion terms in (1.3) and linearize the fast-reaction terms, then we get a matrix D u f (u 1 , u 2 , v, 0) ∈ R 2×2 (1.4) which always has a zero eigenvalue. This means that the classical normal hyperbolicity condition [22] for fast-slow ODEs fails. Yet, the intuition is that the diffusion terms should help to still obtain a well-defined fast-reaction limit as ε → 0.
In the literature, several results for such fast-reaction limits for various PDEs exist. One of the first works is the paper [29] presenting the fast-reaction limit in a system of one parabolic and one ordinary differential equation. A reaction-diffusion system which models a fast reversible reaction between two mobile reactants was then treated in [4] and the limiting problem yields a nonlinear diffusion term. Also fast irreversible reactions (in which two chemical components form a product) were considered, where the limiting system is a Stefan-type problem with a moving interface at which the chemical reaction front is localized [5]. Furthermore, in [27] a dynamical boundary condition has been interpreted as a fast-reaction limit of a volume-surface reaction-diffusion system. It turns out that when the fast-reaction system has three components, the limiting system are often of two types: free boundary problems [42] and cross-diffusion systems [6], which also arise in population dynamics [11,17,30]. In this framework, individuals of one or more species exist in different states and the small parameter represents the average switching time. In our context, it is known that the limiting system of (1.3) is a cross-diffusion PDE [28]. The limiting triangular cross-diffusion system [30] is given by on Ω, (1.5) where the quantities u(t, x), v(t, x) ≥ 0 again represent the population densities of the same two species at time t and position x, confined and competing for resources on a bounded domain Ω ⊂ R N .
The coefficients as above and are all supposed to be positive constants. The model (1.5) is known as Shigesada-Kawasaki-Teramoto (SKT) model as it was initially proposed in [45] in 1979 to account for stable non-homogeneous steady states in certain ecological systems. These states describe, for suitable parameters sets, spatial segregation that is a situation of coexistence of two competitive species on a bounded domain. Historically, the SKT model was first postulated without any reference to a larger system such as (1.3). In 2006, the fast-reaction system (1.3) was introduced in [30] to approximate bounded solutions to the cross-diffusion system (1.5) In this paper, we aim to contribute the understanding of both variants of the SKT model: the three-component fast-reaction PDE as well as its cross-diffusion singular limit. A primary focus from the ecological and mathematical viewpoints is evidently played by the steady states (or equilibria) of both of the models. So we briefly recall the well-known structure of homogeneous steady states starting from the cross-diffusion SKT model (1.5).
-weak competition or strong intra-specific competition, namely a 1 Under the additional condition on the growth rates b 1 /a 2 < r 1 /r 2 < a 1 /b 2 , without diffusion the coexistence homogeneous steady state is stable, while the non-coexistence ones are unstable. With standard diffusion in a convex domain and with zero-flux boundary conditions, any non-negative solution generically converges to homogeneous one, and this implies that the two species coexist but their densities are homogeneous in the whole domain [36]. With cross-diffusion, the model exhibits stable non-homogeneous steady states if d 1 , d 2 are small enough or d 12 is large enough [30]. If r 1 /r 2 < b 1 /a 2 or r 1 /r 2 > a 1 /b 2 , the coexistence homogeneous steady state is not positive and the non-coexistence ones are stable for the reaction part (in absence of diffusion).
-strong competition or strong inter-specific competition, namely a 1 a 2 − b 1 b 2 < 0. When a 1 /b 2 < r 1 /r 2 < b 1 /a 2 , the coexistence homogeneous steady state is unstable, while the noncoexistence ones are stable. With only standard diffusion, in a convex domain and with zero-flux boundary conditions, it has been shown that if positive and non-constant steady states exist, they must be unstable [36], and numerical simulations suggest that any non-negative solution generically converges to either (ū, 0) or (0,v), that is, the competitive exclusion principle occurs between the two species. However, adding the cross-diffusion term, it is numerically shown that even if the domain is convex, there exists stable non-homogeneous solutions exhibiting spatially segregating coexistence when d 12 is suitable large [32]. Analogously to the previous regime, if r 1 /r 2 < a 1 /b 2 or r 1 /r 2 > b 1 /a 2 , the coexistence homogeneous steady state is not positive and the non-coexistence ones are stable for the reaction part (in absence of diffusion).
We remark that the stationary problem of the classical Lotka-Volterra competition model with linear diffusion endowed with Dirichlet boundary conditions has been extensively studied (see [3,12,13] and references therein), and it shows different features to the SKT cross-diffusion model. For the fast-reaction PDE (1.3), the homogeneous steady state is where u * , v * have the same expression as in the cross-diffusion case. In particular, the homogeneous steady state turns out to be independent of ε, which is a nice starting point for a comparative study as it remains to investigate also the non-homogeneous steady states. Before proceeding to summarize the main results of our work, we briefly put it into context with the existing literature. Both systems, the cross-diffusion model and the fast-reaction model, are interesting from different mathematical viewpoints. From the modeling point of view the justification of cross-diffusion terms by means of semilinear reaction-diffusion systems including simple reactions and linear diffusion is fundamental to the understanding of the hidden processes that they can capture. Vice versa, their approximation with simpler models (not in the number of equations but due to the possibility to remove the nonlinearity in the diffusion) is useful both for the analysis and the numerics.
On the one hand, theoretical results require sophisticated techniques. Regarding the cross-diffusion system, existence, smoothness and uniqueness of solutions have been widely investigated (see [1,2,16,23,33] and the references therein). The convergence of the solution of the stationary problem have been treated in [32], while the convergence of the solutions of the fast-reaction system (1.3) towards the solutions of the cross-diffusion system (1.5) was shown in [10] in dimension one, and then generalized to a wider set of admissible reaction terms and in any dimension [19]. Similar results were obtained for a class of non-triangular cross-diffusion systems [41], assuming the same time-scale for all the species. Finally, it has been proven that the limiting system inherits the entropy structure with an entropy that is not the sum of the entropies of the components [14].
On the other hand, the capability of the cross-diffusion system (1.5) to model the spatial segregation of competing species is related to the appearance of non-homogeneous solutions. Although the system (1.5) does not have an activator-inhibitor structure (the most important mechanism in the Turing instability theory for pattern formation), the cross-diffusion term turns out to be the key ingredient to destabilize the homogeneous equilibrium [24,30,49]. In this framework bifurcations diagrams are useful to present and explore the bifurcating branches and the non-homogeneous steady state solutions. Numerical continuation techniques are going to allow us to obtain a more global picture far from the homogeneous branch. The bifurcation diagram of the triangular cross-diffusion system on a 1D domain was presented in [30]; the existence of these non-homogeneous steady states significantly far from being perturbations of the homogeneous solutions, was proved in [7] applying a computer assisted method. A mathematically rigorous construction of the bifurcation structure of the three-component system was obtained in [9], while the convergence of the bifurcation structure on a 1D domain was also investigated in [32] with respect to two different bifurcation parameters appearing in the diffusion part, both in the weak and strong competition case. Only recently, the influence of the additional cross-diffusion term in the system has been investigated, combining a detailed linearized analysis and numerical continuation. To the best of the authors' knowledge, not much has been done on 2D domains: the possible pattern admitted by the non-triangular SKT model was explored in [25], but the bifurcation structure on 2D domain is not known.
To study the bifurcation structure, we use and extend the continuation software for PDEs pde2path [20,46,48], based on a FEM discretization of the stationary problem. Since the general class of problems which can be numerically analyzed by the software does not include the cross-diffusion term appearing in (1.5), it requires an additional setup to be able to compare the fast-reaction PDE results with the singular limit cross-diffusion systems. Yet, it has been shown in the past for other classes of PDEs related to classical elliptic problems that pde2path can be extended beyond its standard setting [37]. This is the key reason, why we have chosen to carry out numerical continuation within this framework.
Here is a summary of the main contributions of this paper.
-We explain how to set up the continuation software pde2path in order to treat cross-diffusion terms.
Then we cross-validate and extend previous computations for various regimes of the time-scale separation parameter for 1D domains.
-We compute a new bifurcation diagram with respect to a parameter appearing in the reaction part. We show how the bifurcation structure of the fast-reaction system modifies when ε → 0. In particular, a novel "broken-heart" structure of non-homogeneous steady state bifurcation branches is observed in the singular limit.
-Then we compute and interpret the various bifurcation diagrams and the non-homogeneous solutions of the triangular SKT model on a 2D rectangular domain. Also in this case the convergence towards the singular limit is analyzed carefully.
-We show a link between the computed bifurcation diagrams in 1D and 2D domains, and the Turing instability analysis as a tool to fully understand and validate the numerical continuation calculations.
The paper is organized as follows. In Section 2 we present the the numerical continuation results obtained with pde2path in the 1D case: we provide a detailed picture of the different types of stable unstable non-homogeneous solutions to the cross-diffusion system arising at each bifurcation point, and we quantify the convergence of the fast-slow system to the cross-diffusion one. Section 3 is devoted to the 2D case: we show the bifurcation diagram and different patterns, as well as the bifurcation diagrams of the fast-reaction system. In Section 4 some concluding remarks can be found. Appendix A contains the pde2path setup for cross-diffusion systems, while in Appendix B we report the Turing instability analysis of the cross-diffusion system for reference and validation.

Numerical continuation on a 1D domain
The numerical analysis of systems (1.5) and (1.3) is performed using the continuation software pde2path, originally developed to treat standard reaction-diffusion systems and here adapted to investigate crossdiffusion systems. The software setup required for system (1.5) can be found in Appendix A. For the numerical results we use the set of parameter values widely used in literature for the weak competition regime [7,9,30,32], here reported in Table 1. We set d 1 = d 2 =: d and use d as one main bifurcation parameter. It follows that which is one possible starting point for the continuation. In this section we consider a 1D domain (interval) Ω = (0, 1), as in [7,8,9,30,32]. We provide in Section 2.1 a detailed characterization of different steady state types of the cross-diffusion system (1.5). In Section 2.2 we study the convergence of the bifurcation structure of the fast-reaction system (1.3) to the cross-diffusion system (1.5), taking the standard diffusion coefficient d as bifurcation parameter. In this framework, reproducing such diagrams is crucial to test the numerical continuation software; we add to the existing literature a quantification of the convergence of the bifurcation structure when ε → 0, a clear picture of the behavior of the non-homogeneous solutions along the bifurcation branches and their stability properties, and the corresponding bifurcation diagram in the L 1 -norm. Finally, in Section 2.3 we consider the growth rate r 1 as bifurcation parameter: also in this case we present the  Table 1: Parameter values used in numerical continuation. The set r i , a i , b i , (i = 1, 2) corresponds to the weak competition case (or strong intra-specific competition).
bifurcation diagram of the cross-diffusion system and the stable steady states appearing beyond the usual range of parameters, as well as the convergence of the bifurcation structure of the fast-reaction system to the cross-diffusion one.

Bifurcation diagram of the cross-diffusion system
We numerically compute the bifurcation diagram of the cross-diffusion system (1.5) using the continuation software pde2path and we provide a clear picture of the different solution types.
In Figure 1 the bifurcation diagram of the cross-diffusion system is plotted with respect to different quantities on the y-axis, namely v(0), corresponding to the density of the second species at the left boundary value of the domain, and the L 1 -norm of species u. From now on, the homogeneous solution is denoted by the black line, while the other branches correspond to non-homogeneous solutions originating by successive bifurcations. The bifurcations corresponding to branch points are marked with circles, while fold/limit points are marked with crosses. Thick and thin lines denote stable and unstable solutions, respectively. Note that in the (d, v(0))-plane at each bifurcation point, two separate branches appear, corresponding to two different solutions. In the (d, ||u|| L 1 )-plane, the two branches are overlayed, since the solutions on the branches are symmetric on the domain. The shape of the non-homogeneous steady states originating along the branches are reported in Figure 2 Starting from d = 0.04 and decreasing its value, we can see that the homogeneous solution is stable and no other solutions are present. At B 1 the homogeneous solution undergoes to a primary bifurcation losing its stability, and two stable non-homogeneous solutions appear (blue lines). Along those branches the density of species v is greater in a part of the domain, while species u occupies the other one. Note that the stable solutions on those branches are symmetric (solid and dashed lines in Figure 2a). At B 2 , further non-homogeneous solutions appear (red lines), initially unstable; the solutions are again symmetric but now the species density is concentrated either in the central part of the domain or close to the boundary. Those solutions become stable for smaller values of d at a further bifurcation point. We observe the same behavior at the successive bifurcation points from the homogeneous branch: new branches appear (green, yellow and cyan), one each branch the new solutions add half a bump to their shape (Figures 2c-2e). Along the bifurcation branches, the differences between peaks and valleys increases, as the bifurcation parameter d becomes smaller. Finally, for small values of the bifurcation parameter d there can be many different locally stable non-homogeneous solutions. Moreover, there are bifurcation curves connecting three different branches of the homogeneous solutions (magenta and orange lines): along the branches the solution changes shape in order to match the solution profile on the primary branches (Figures 2g, 2h). Black circles on the homogeneous branch for small values of d indicate the presence of further bifurcation points that we have not continued.

Convergence of the bifurcation diagram
The idea to study bifurcation diagrams and the associated singular limit bifurcation diagram as ε → 0 has been successfully carried out for fast-slow ODEs in several examples [26,31]. Yet, more systematic studies for fast-reaction limits for PDEs are missing. Here we numerically compute the convergence of the bifurcation structure of the fast-reaction system (1.3) to the one of the cross-diffusion system (1.5), as already reported in [32], quantifying the convergence of the bifurcation points.
In Figure 3   (or they exist for small values of the parameter d and it is difficult to numerically detect them), and the homogeneous steady state is stable for all the values of the bifurcation parameter d. For ε = 0.05, the homogeneous steady state becomes unstable for small values of d but we can see that the bifurcation structure corresponding to the fast-reaction system is already qualitatively similar to the one of the cross-diffusion system, but it is squeezed into a small region near d = 0. The stability properties also match with the cross-diffusion bifurcation structure. For smaller values of ε the bifurcation structure is expanding to the right, towards the bifurcation structure of the cross-diffusion system. Note that Figure 3c, obtained with ε = 0.001, is almost indistinguishable from Figure 3d corresponding to the cross-diffusion system (ε = 0). Hence, from the viewpoint of the bifurcation structure, the threecomponent fast-reaction system is indeed a good approximation for the cross-diffusion one (at least concerning stationary steady states).
In Figure 4 we quantitatively show the convergence of the first three bifurcation points on the homogeneous branch, namely B 1 , B 2 , B 3 , using the difference between the bifurcation values d 0 Bi and d ε Bi . The order of convergence is approximately one.

Changing the bifurcation parameter
After we have set up the system in pde2path, we can also change the bifurcation parameter. One possible choice is r 1 , the growth rate of population u. This parameter appears in the reaction part of the system and the homogeneous coexistence state (u * , v * ) depends on its value. Note that, without diffusion, in the weak competition case (namely a 1 a 2 − b 1 b 2 > 0) the homogeneous coexistence state is positive (and stable) when b 1 /a 2 < r 1 /r 2 < a 1 /b 2 , otherwise it is not meaningful and the noncoexistence states (ū, 0) and (0,v) are stable. With cross-diffusion, it has been shown in [8] that the homogeneous equilibrium can be destabilized and stable non-homogeneous solutions arise, which survive in a region of the parameter space in which the homogeneous solution is no longer admissible, namely r 1 /r 2 > a 1 /b 2 .
In Figure 5 we report the bifurcation diagram of the cross-diffusion system with respect to the bifurcation parameter r 1 and some solutions, obtained with the set of parameters of Table 1 and d = 0.02. Note that it is possible to compare it with the one reported in Figure 1a, since they are different cross-sections of a two-parameter bifurcation surface. The bifurcation diagram is composed  Figure 2: (a)-(h) Different solution types along the branches in the bifurcation diagram (i); the concentrations of u, v are shown in black and blue respectively. Solid lines correspond to points (marked with yellow triangles) located above the homogeneous branch in the bifurcation diagram, while dotted lines to points located below (marked with yellow squares). In the bifurcation diagram, thick lines correspond to stable solutions and thin lines to unstable ones.   by three rings and the solutions profile is shown in Figures 5a-5d. In particular, the two outer (blue) rings contain qualitatively similar solutions (Figures 5a and 5b). Furthermore, Figure 5b shows a stable non-homogeneous steady state corresponding to a parameter value outside of the usual weak competition regime. Furthermore, the (red) branches originating from the second bifurcation point are non-symmetric regarding the stability properties even if they correspond to symmetric solutions (with respect to the homogeneous one) on the domain (Figures 5c and 5d).
In Figure 6, we show how the bifurcation diagram with respect to r 1 behaves when the fast timescale parameter ε decreases. On the vertical axis the value v(0) is reported. Also in this case we observe convergence of the bifurcation structure of the fast-reaction system to the one of the crossdiffusion system of Figure 5e (well approximated by taking ε = 10 −4 , corresponding to Figure 6f). However, we observe that the bifurcation structure is not only expanding as in the previous section, but it is qualitatively changing as ε decreases. In particular, in Figure 6a and 6b there are only two non-homogeneous branches connecting two primary branch points on the homogeneous branch and forming a bifurcation ring. For smaller values of ε other bifurcation points (and consequently other rings) appear inside the main bifurcation ring; the inner rings expand while the outer ring folds in the middle part forming a heart shape. Then the outer ring interacts with the inner ring and separates into two rings, giving rise to the bifurcation structure of the cross-diffusion system, i.e., the heart-shaped structure breaks.
Also in this case we quantitatively show the convergence of the first bifurcation points on the homogeneous branch in Figure 7. Although the bifurcation structure is qualitatively changing, the convergence rate is comparable to Figure 4. This clearly shows that locally one can expect a convergence rate near the homogeneous branch. In fact, this leads one to the conjecture that even the global diagram could be captured asymptotically by a convergence rate towards the singular limit but since we have only captured part of the full diagram, this is hard to validate completely numerically. Furthermore, many different convergence metrics are conceivable if one moves beyond single points so we leave this conjecture as a topic for future work.   We clearly observe a highly non-trivial deformation of the bifurcation diagram as ε is decreased, starting from a ring and then heart-shape, we eventually have a "broken-heart" structure leading to two rings in the cross-diffusion limit.

Numerical continuation on a 2D rectangular domain
In this section we consider a 2D rectangular domain with edges of length L x = 1 and L y = 4. This choice reduces the presence of multiple branch points. As in the previous section, we set the standard diffusion coefficient d as our main bifurcation parameter. The bifurcation analysis close to the homogeneous branch can be performed (see Appendix B). It allows to compute the values of the parameter d at the bifurcation points on the homogeneous branch, which can be compared to the values numerically obtained. It also predicts the shape of the steady state for each bifurcation point close to the homogeneous branch, by looking at the eigenvalues of the Laplacian at the bifurcation point and the associated eigenfunctions [35,39].
In Figure 8 we show part of the bifurcation diagram close to the homogeneous branch with respect to the L 2 -norm of the species u. As in the 1D case, for decreasing values of d the homogeneous solution destabilizes at the first bifurcation point, and stable non-homogeneous steady states appear (blue branch). For smaller values of d other bifurcation points occur. We show in the figure just some of the successive branches and we report the corresponding solution at the gray points (Figures 8a-8f), which turn out to be unstable. Note that certain cross-sections of the 2D-solutions have a similar shape as the steady states of 1D steady states.
Different from the 1D case, in Figure 9 the enlargement of the initial part of the first (blue) branch is reported. It shows that the first bifurcation point is supercritical, but then two successive fold bifurcations (indicated with a cross) appear leading to multi-stability of solutions: in a small range of the bifurcation parameter the system admits four (two symmetric) stable non-homogeneous steady states. Their shape is reported in Figures 9a-9d. Close to the homogeneous branch (Figures 9a and 9b) the steady states have half a bump on both edges (in agreement with the eigenvalues corresponding to the first bifurcation point, see Appendix B), while going further along the branch the shape modifies.
Finally, Figure 8 shows three different branches, computed far away from the homogeneous one. The first (blue) branch undergoes a further bifurcation and a secondary stable branch arises (marked in green). Then, on this branch a Hopf bifurcation point has been detected. Note that this type of bifurcation is not present in the 1D case in the weak competition regime and yields the existence of time-periodic solutions. Furthermore, the bifurcation diagram in the 2D case seems to be even more intricate compared to the 1D case, many branches are curly and they swirl back and forth. The solution along the branches deforms (see Figures 10a-10i), and various patterns are evident. Stripes occur in 10b, 10f and spots in 10g, 10i, although they are unstable.
x y (f)

Convergence of the bifurcation structure
As in the 1D case, we investigate how the bifurcation structure deforms on a 2D domain with respect to the time-scale separation parameter ε. In Figure 11 the bifurcation diagrams for different values of ε are reported showing the convergence of the bifurcation structure of the fast-reaction system (1.3) to the cross-diffusion one (1.5), shown in Figure 12i. For the sake of clarity of the visualization we only show stable branches. For ε = 0.05 the instability region of the homogeneous branch is reduced, and the two primary branches roughly sketch the cross-diffusion one: the first bifurcation point is supercritical, the first branch is stable and the other (red) represents a stable part. However, there is no a secondary bifurcation branch giving rise to a Hopf point. For decreasing values of ε, the cross-diffusion bifurcation structure is well recognizable.
Finally, in Figure 12 we show three different stable stationary solutions of the fast-reaction system with ε = 0.0001: they corresponds to points on the blue, green and red branches of Figure 11d, showing different (stable) patterns. It is interesting to observe the distribution of the two different states of species u on the domain. It turns out that the "excited individuals" and the competing species v are more abundant in the same region of the domain. Furthermore, the steady state solutions satisfy the algebraic equation derived by a quasi-steady state ansatz from the three-species system (1.3) In addition, we have found other stable non-homogeneous solutions of the cross-diffusion system, and the corresponding steady states are stripe patterns.
x y  (a) u1

Conclusion and Outlook
In this paper we have investigated the bifurcation structure of the triangular SKT model and of the corresponding fast-reaction system in 1D and 2D domains in the weak competition case via numerical continuation. Despite the fact that part of the bifurcation structure of this system in 1D has been already computed [7,9,30,32], the key points of our work can be here summarized.
-We have adapted the continuation software pde2path to treat cross-diffusion terms. Providing the code in the appendix, this work can be used as a guide to implement such class of problems.
-We have reproduced the already computed bifurcation diagrams for the triangular SKT model on a 1D domain with respect to d. Even though this is not per se a new result, we have quantitatively checked, how accurate the computation of bifurcation points is. It is worthwhile to note that we have also provided information about the stability of non-homogeneous steady states, which are in agreement with previous results [7,30] and confirms the reliability of our new software setup for cross-diffusion systems.
-Once the software setup has been established, we have changed the bifurcation parameter to obtain novel structures. We have selected as new bifurcation parameter the growth rate of the first species, which appears in the reaction part. We have shown that the bifurcation structure qualitatively changes as the small parameter tends to zero leading to a "broken-heart" bifurcation structure.
-With respect to both of the considered bifurcation parameters, we have also provided a novel precise quantification of the convergence of the bifurcation points on the homogeneous branch of the fastreaction system to the cross-diffusion ones.
-A major new contribution is that we have computed the bifurcation diagram and the non-homogeneous steady states of the triangular SKT model on a 2D domain (rectangular). This case is intricate; the resulting bifurcation diagram is not as clear as in 1D. We have highlighted the main characteristics of the diagrams, and we have presented different types of non-homogeneous steady states, with a focus on stability and pattern formation. We have seen that solutions can exhibit spots or stripes, depending on the parameters, but such solutions turn out to be unstable (at least as far as we have computed the bifurcation branches).
-We have provided a link between the computed bifurcation diagrams in 1D and 2D domains, and the Turing instability analysis as a tool to fully understand and validate the software output. It is worthwhile to note that the results obtained with the Turing instability analysis only provide insights close to the homogeneous branch. However, the global bifurcation structure has to be numerically computed to achieve a full comprehension of the possible (stationary) outcomes of the system.
-Our numerical calculations can now provide also guidelines and conjectures for analytical approaches to cross-diffusion systems, e.g., where in parameter space one can expect entropy structures to behave differently, or where multi-stability and deformation of global branches has to be taken into account.
Several research directions arise at this point. On the one hand, it is natural to further combine analytical methods and numerical continuation techniques to the study of the full cross-diffusion systems in order to to understand the role of cross-diffusion terms on pattern formation and their influence on the bifurcation diagram [8]. On the other hand, continuing our work to numerically investigate of the convergence of the bifurcation structure of a four-equation fast-reaction system leading to the full cross-diffusion SKT system, not only the triangular one, is a straightforward continuation of this work. In this context, analytical proofs of convergence of the solutions of the four species system with two different time-scales to the non-triangular cross-diffusion one is an open problem [17].
Furthermore, the 2D case we have computed has shown that Hopf bifurcations can appear, which can give rise to time-periodic solutions. Their presence could then also be investigated further using the continuation software pde2path, in order to provide a clearer picture on possible long-time asymptotic dynamics of the model. Another direction is to explore other cross-diffusion systems which have been derived by time-scale arguments [11,17,18], or simply proposed to describe different processes [34,40]. To this end, we have shown that the continuation software is easily applicable to different models and it can be adapted to treat systems that do not exactly belong to the class of problems for which it has been developed. For instance, another direction could be its extension to other non-standard diffusion processes such as fractional reaction-diffusion equations [21].
In [44] the software setup for the quasilinear Allen-Cahn equation is explained, and only recently this approach was extended to treat a chemotaxis reaction-diffusion system involving a quasi-linear cross-diffusion term of the form ∇·(u∇v) [47]. However, the previous approach is not directly applicable to the cross-diffusion system (1.5), since the cross-diffusion term is not written in divergence form. Then we need to rewrite it as The steady state problem reads and on the FEM level it becomes where K is the standard one-component Neumann-Laplacian (stiffness matrix), K 21 (v)u and K 12 (u)v implement ∇ · (c(v)∇u) and ∇ · (c(u)∇v) respectively, while F 1 , F 2 belong to the reaction part. In detail, depend on v and on u. Hence, they have to be computed at each step. In the OOPDE setting, we can employ the routine assema to compute those matrices, but this needs c(v) andc(u) on each element center, which is obtained interpolating u and v from the nodes to the element centers, as it can be seen in Listings 1 and 2, which show the main files implementing the triangular cross-diffusion system (1.5) in pde2path.
If we rewrite (B.1) as a second-order polynomial in d it is possible to locate the bifurcation points solving det M * k = 0, we obtain In Table 2 we report the bifurcation values d B (λ k ) corresponding to the eigenvalue λ k of the Laplacian for the 1D domain (0, 1) and obtained by formula (B.3), for different values of n, compared with the numerical values estimated by the software pde2path. The obtained values are in good agreement. The scenario is more complicated in a rectangular 2D domain, reported in Table 2 for different values of n and m. Also in this case we compare the values obtained by formula (B.3) with the numerical values. The software is not able to detect all the bifurcation values, but it fails when they are too close to each other or with multiplicity more than one. Note that the values are ordered with respect to the bifurcation values d B , which does not translate into an order of the indices n, m: for instance the couple (0, 3) corresponds to a smaller bifurcation value than (0, 7). The reason is evident in Figure 3, due to a non-monotonicity of the function d B (λ k ) w.r.t. λ k .