Traveling waves of an FKPP-type model for self-organized growth

We consider a reaction–diffusion system of densities of two types of particles, introduced by Hannezo et al. (Cell 171(1):242–255.e27, 2017). It is a simple model for a growth process: active, branching particles form the growing boundary layer of an otherwise static tissue, represented by inactive particles. The active particles diffuse, branch and become irreversibly inactive upon collision with a particle of arbitrary type. In absence of active particles, this system is in a steady state, without any a priori restriction on the amount of remaining inactive particles. Thus, while related to the well-studied FKPP-equation, this system features a game-changing continuum of steady state solutions, where each corresponds to a possible outcome of the growth process. However, simulations indicate that this system self-organizes: traveling fronts with fixed shape arise under a wide range of initial data. In the present work, we describe all positive and bounded traveling wave solutions, and obtain necessary and sufficient conditions for their existence. We find a surprisingly simple symmetry in the pairs of steady states which are joined via heteroclinic wave orbits. Our approach is constructive: we first prove the existence of almost constant solutions and then extend our results via a continuity argument along the continuum of limiting points.


Motivation and result
The mechanics of tissue-growth have drawn the attention of the scientific community. A central question is, how the cells are organized, how they react to and communicate with their environment on the microscopic level, and how their behavior during the growth phase gives rise to distinct macroscopic structures. Mathematical models can help to understand these processes and works regarding organoids, wound healing or tumor growth are abundant (Montes-Olivas et al. 2019;Mammoto and Ingber 2010;Falco et al. 2021;d'Alessandro et al. 2021;Jain et al. 2021). However, for most of these models, our numerical skills far predominate the possibility to analyze them rigorously. Hence, for understanding the basic mechanics of the underlying biological processes, the need for simplified models arises.
Especially when studying spatiotemporal effects and macroscopic pattern formation, reaction-diffusion systems and their traveling waves have proven insightful. One of the oldest and most studied models is the FKPP-equation (Fisher 1937;Kolmogorov et al. 1937), describing the advance of an advantageous population. The arise of more complex spatial patterns due to the instability of a homogeneous state was first described in Turings groundbreaking paper The Chemical basis of Morphogenesis (Turing 1952). More recently, systems of Keller-Segel type have been studied extensively, where growth, movement and self-organization of a population are driven by chemotactic guidance (Keller and Segel 1971;Perthame 2004;Painter 2019). Broad introductions to mathematical modeling of pattern formation in developmental biology have been written by Painter (2019) and Othmer et al. (2009), among others.
The group of Hannezo et al. proposed A Unifying Theory of Branching Morphogenesis in epithelial tissues (Hannezo et al. 2017). They introduced a stochastic model, related to branching and annihilating random walks (Cardy and Täuber 1996). In this model, a branched structure is represented by a network. This network undergoes stochastic growth dynamics, where each branch of the network grows independently from the others and follows a set of simple, local rules. At its tip, each branch elongates or splits up at certain rates and these tips are called active. When an active tip comes too close to a different branch, it irreversibly ceases any activity and becomes inactive. The numerical results of Hannezo et al. reveal that this stochastic growth process self-organizes: the active tips concentrate at the boundary of the network and form a rather sharp layer of growth. The center of the network is static and -rather surprisingly -exhibits a homogenous geometry, in particular a constant density of branches. Remarkably, as mentioned by the authors, this model self-organizes without any signaling gradients. Even a directional bias of the branches can be achieved, as the result of an appropriate spatial boundary. Moreover, the authors observed that their simulations were in good agreement with biological data from mammary glands, kidneys and the human prostate (Hannezo et al. 2017).
To study their model analytically, Hannezo et al. proposed the following system, which corresponds to the diffusive limit of the above stochastic dynamics. We restrict ourselves to the one-dimensional case. Due to a simple linear rescaling ("Appendix C"), we only need to consider the normalized reaction-diffusion system (1.1) Here, A, I : R × R + → R + are the densities of active particles and inactive particles. The diffusion term describes the movement of the active particles, all other terms encode a growth process where the active particles eventually become inactive: the active particles branch with rate 1, produce inactive particles with rate r ≥ 0, and become inactive upon collision with either an active or an inactive particle. The active particles grow logistically, which implies that the inactive particles grow at most exponentially. The resulting simple time-dependent bounds on A and I yield uniqueness and existence of smooth solutions, a suitable fixed-point theorem is presented in chapter 14 of Smoller (1994). More details about the underlying stochastic processes together with a non-rigorous derivation of this PDE can be found in Hannezo et al. (2017). Note that without the inactive particles, i.e. when I = 0, the remaining equation for A reduces to the well-known FKPP-equation (Fisher 1937;Kolmogorov et al. 1937).
The System (1.1) can be interpreted as a twofold degenerate Keller-Segel system (Keller and Segel 1971;Arumugam and Tyagi 2020): the active particles are not guided by a chemotactic gradient, but explore the space solely diffusively, and the inactive particles do not diffuse at all. Still, simulations of System (1.1) show that general solutions of (1.1) self-organize, a phenomenon which is typical for many different Keller-Segel systems (Painter 2019). The invading front of the system converges to a fixed shape: a pulse of active particles, that represents a layer of growth, is accompanied by a monotone wave of inactive particles, the resulting static tissue, as demonstrated in Fig. 1. This self-organization of the Reaction-Diffusion System (1.1) resembles that of the stochastic dynamics.
For a wave speed c > 0, a right-traveling wave solves Eq. (1.1) via the Ansatz A(x, t) = a(x − ct), I (x, t) = i(x − ct). We substitute z = x − ct, such that any traveling wave must be a solution of Fig. 1 Simulation of the Reaction-Diffusion System (1.1) for r = 0. Given a small initial heap of active particles A(x, 0) = 1/2 exp(−x 2 ) and I (x, 0) = 0, two identical traveling fronts arise, the right one is shown. After the separation of the two fronts away from the origin, the density of the remaining inactive particles is given by I = 2 and the front moves asymptotically with speed c = 2 0 = a zz + ca z + a − a(a + i), 0 = ci z + a(a + i) + ra. (1. 2) The occurrence of these seemingly stable traveling waves is quite surprising, since the System (1.1) features a continuum of steady state solutions: which is due to the fact that the inactive particles do not degrade. This continuum of steady states represents the difficulty when studying the system: we first need to find out which limiting states are chosen by the growth process. Hannezo et al. presented a rich discussion of the Wave-Equation (1.2) along with numerics and several heuristics that show a deep connection with the original FKPPequation, and predicted some of the following results. The goal of this paper is to give necessary and sufficient conditions for the existence of such traveling wave solutions and to analyze the shape of the wave form. Our main result characterizes a family of pulled traveling waves: Theorem 1.1 Let r ≥ 0, c > 0 and consider the System (1.1) and its traveling wave solutions given by (1.2). Set i c : there exists a unique bounded and positive traveling wave a, i ∈ C ∞ (R, R 2 ) with speed c such that (1.5) The function i(z) is decreasing, whereas a(z) has a unique local and global maximum. If c 2 4 + i +∞ − 1 = 0, then convergence as z → +∞ is sub-exponentially fast and of order z · e − c 2 z . If c 2 4 + i +∞ − 1 > 0, then convergence as z → +∞ is exponentially fast. Convergence as z → −∞ is exponentially fast in all cases. The corresponding rates are Moreover, these are all bounded, non-negative, non-constant and twice differentiable solutions of Eq. (1.2).
Notice that this result is independent of the reproduction rate r , which affects the shape of the wave, but neither its limits nor the minimal speed of a positive solution. Hence, all non-negative and bounded traveling waves resemble the ones depicted in Fig. 2, consisting of a pulse of active particles and a monotone wave of inactive particles. These traveling wave solutions share many similarities with classical FKPPwaves of a single type of particles. Among other mathematical aspects, this will be discussed at the end of the paper, in Sect. 9. Notably, Theorem 1.1 analytically connects two continua of fixed points via a continuum of traveling waves. Our constructive approach is a novelty: we first prove the existence of almost constant solutions and then continuously deform these solution along the continuum of possible limits. Figure 1 shows a simulation of the System (1.1), starting with a small initial amount of active particles. After a short transition phase, we observe a front with fixed shape. Asymptotically, it equals the unique traveling wave with limits i −∞ = 2, i +∞ = 0 and speed c = 2, which is the minimal possible wave speed for this pair of limits. We observed this behavior for all compact initial data that we chose. Moreover, this wave seems to be stable against perturbations, as briefly discussed in the concluding Sect. 9. Even though it is only a first step into this direction, this paper sheds light at the ability of the Growth-Process (1.1) to self-organize and at the robustness of this mechanism, e.g. against errors of individual particles. Our theoretical analysis fortifies the numerical and biological findings of Hannezo et al., where a simple set of local rules organizes the growth of a complex epithelial structure. The underlying assumption of a logistic growth is quite natural, so similar rules might drive and regulate other growth processes as well, without the need for guiding gradients.

Outline of the paper
A sketch of the central ideas and techniques is presented in Sect. 2. The identity i −∞ + i +∞ = 2 is proved in Sect. 3. The asymptotic behavior around the stable and unstable set of the traveling waves is analyzed in Sect. 4. A non-negative trapping region of a lower-dimensional sub-system is analyzed in Sect. 5. We use our knowledge about the sub-system to construct a suitable attractor of the full system in Sect. 6. Then, we connect the unstable manifold of the unstable set with this attractor, see Sect. 7. We complete the proof of Theorem 1.1 in Sect. 8. In Sect. 9, we highlight the similarity of the traveling waves with those of the original FKPP-equation and give a short outlook at their stability.

Identifying the correct limits
We reformulate the System (1.2) for a traveling wave as an equivalent system of firstorder ODEs. Denoting differention with respect to z by a prime, we introduce the auxiliary variable a = b, so that (1.2) becomes , then it is also in C ∞ by a simple induction. This equation has a continuum of non-negative fixed points, similar to that of the PDE, cf. (1.3): Thus, in the first place, we need to find out which of these fixed points can be considered as limits of right-traveling waves. Any bounded and non-negative solution of System (2.1) can not be periodic and must converge since ci = −a(a + i + r ) ≤ 0. It is now evident that the limits at z = ±∞ must be fixed points of Eq. (2.1), thus we denote them as (a, b, i) = (0, 0, i ±∞ ). Under mild assumptions regarding integrability, we can interrelate two different points on a given traveling wave, see Sect. 3. Most importantly, this leads to the correspondence of the limits In view of this, monotonicity of i implies that i −∞ ∈ (1, 2] and i +∞ ∈ [0, 1). The fixed points of the ODE System (2.1) are not isolated, hence its Jacobian D is degenerate there. It is easily verified that D is given by Hence, we can not apply the classical Theorem of Grobmann-Hartmann to linearize the asymptotic behavior around the fixed points. We apply center manifold theory to work out the higher moments of the approximation, see Sect. 4. The center manifold coincides with the continuum of fixed points. This implies that asymptotically, there is no flow along the direction of the eigenvector (a, b, i) = (0, 0, 1) with zero eigenvalue.
Hence, the asymptotic flow around any fixed point is two-dimensional and the stability of the fixed point (0, 0, K ) is dictated by the two eigenvalues λ ± . When K > 1, the fixed point is unstable, while for K < 1, it is stable. At the same time, the analysis of the asymptotic behavior also yields a necessary condition on the speed c of a non-negative wave. A traveling wave can only be nonnegative if a(z) does not spiral while converging to 0. Therefore, the two eigenvalues λ ± at the limiting fixed point must be real-valued. In view of (2.5), for a fixed point (0, 0, K ), this is given if Thus, if the stable fixed point (0, 0, i +∞ ) is the limit of a non-negative traveling wave, where i +∞ ∈ [0, 1), it must by (2.6) further hold that as in Theorem 1.1. In other words, i c is the minimal limiting density of inactive particles that is necessary for the existence of a non-negative traveling wave with speed c.

Construction of a traveling wave
We will explicitly construct a non-negative traveling wave such that the two necessary conditions i +∞ ≥ i c and i +∞ + i −∞ = 2 are fulfilled. Two key features of the model make it tractable: first, the monotonicity of i(z) allows us to investigate the convergence of the sub-system that arises for a fixed value of i, and then lift our result to almost constant solutions of the full system. Second, for extending this result to non-small solutions, we lean on an integral equation that allows us to interrelate two points on a given trajectory. However, the central Proposition 3.2 depends essentially on the logistic growth of the active particles. Apart from this, our general approach seems to be applicable to a broader class of systems. Regarding the ODE System (2.1), our analysis of the flow around the fixed points in Sect. 4 reveals a suitable unstable set and a suitable stable set Each point (0, 0, i −∞ ) ∈ S −∞ has an unstable manifold of dimension one. Its restriction to a ≥ 0 is the only possible candidate for the tail of a non-negative traveling wave as z → −∞. Each point (0, 0, i +∞ ) ∈ S +∞ is Lyapunov stable, which can also be seen in Fig. 3. To begin with, we let (0, , and follow its unstable manifold in positive direction of a. We prove that there exists a finite phase-time z 0 such that b(z 0 ) = 0: the trajectory reaches a local maximum of active particles, again see Fig. 3. We denote it as (a z 0 , 0, i z 0 ). This is carried out in Sect. 7.
Thus, for finding a suitable attractor of S +∞ , we analyze solutions that start in points of type (a, b, i) = (a 0 , 0, i 0 ). We first analyze the lower-dimensional subsystem in coordinates (a, b), imposing a fixed value of i, which is done Sect. 5. We construct a trapping region, wherein a converges and stays non-negative. The monotonicity of i(z) allows us to lift this result to the full system, see Sect. 6. Here, the continuum of fixed points comes at help: we first prove the existence of almost constant solutions, where a 1 and i ∼ i 0 , that stay non-negative and converge. Then, we continuously deform these solutions: the Lyapunov-stability of the limits in S +∞ implies continuity of the entire trajectory up to z = +∞ in initial data. We use this to derive sharp conditions regarding (a 0 , 0, i 0 ) such that the trajectory stays non-negative and converges.
We show that the first local maximum (a z 0 , 0, i z 0 ) along the instable manifold of (0, 0, i −∞ ) does fulfill these conditions, see again Sect. 7. The technique is the same as for proving the identity i +∞ + i −∞ = 2, which is the starting point of our analysis and presented in the next section. The proof of Theorem 1.1 is completed in Sect. 8, where we bring together all the different pieces. The resulting continuous family of solutions is sketched in Fig. 3. a(z) ≡ 0. If not a(z) ≡ 0, then there is at least one local maximum of active particles, which we denote as (a 0 , 0, i 0 ).
At this point, a = b = a 0 (a 0 + i 0 − 1) ≤ 0, so either a 0 = 0 and the wave is constant, or a 0 + i 0 ≤ 1. In the second case case, assume that there is also a local minimum of a(z), denoted as (a m , 0, i m ). Since a(z) vanishes as z → ±∞, we may assume that this be the first local minimum after passing through (a 0 , 0, i 0 ). As before, (a m , 0, i m ) is already a fixed point or a m + i m ≥ 1. Since i(z) is decreasing, a(z) must have been increasing, a contradiction to the assumption that this is the first local minimum after the maximum (a 0 , 0, i 0 ). Thus, there is only one local maximum of active particles, which is also the global one. Further, this implies a = b ∈ L 1 (R). By ci = −a(a + i + r ) ≤ 0, we know that a(a + i + r ) is also in L 1 (R). We integrate b + cb + a = a(a + i) over the finite interval [−M, M], then send the boundaries to ±∞: (3.1) We know that the right-hand is integrable since i ∈ L 1 (R), and that both a(±M) and b(±M) vanish as M → +∞. This implies Hence also a ∈ L 1 (R), since a ≥ 0. Finally, as a sum of integrable terms, also b ∈ L 1 (R).
The following Proposition 3.2 will be used several times to interrelate two points (a 1 , 0, i 1 ), (a 2 , 0, i 2 ) on a traveling wave, where b i = 0. By the previous Lemma, the necessary conditions regarding integrability are always verified for non-negative and bounded solutions.

Proposition 3.2 Let a(z), b(z), i(z) be a smooth and bounded solution of the ODE
The following three identities hold: Proof Any solution of the ODE System (2.1) also fulfills the original Wave Equations (1.2). We integrate these over [z 1 , z 2 ], substitute a = b and use that b(z i ) = 0. This directly proves (3.3) and (3.4). Regarding Eq. (3.5), note that by integration by parts: (3.6) Remark Equations (3.3) and (3.4) encode a mass transfer from the active to the inactive particles and are not specific for the chosen reactions. It is the quadratic Eq. (3.5) that relies on a logistic saturation mechanism, we do not see a (direct) way to generalize this result.
Given Proposition, the identity i −∞ + i +∞ = 2 is a mere

Corollary 3.3 (Limits of traveling waves) Let a(z), b(z), i(z) be a non-negative and bounded traveling wave that solves the ODE System (2.1), and denote its limits as
or the following identity holds: Proof We apply Proposition (3.2). Since a(±∞) = 0, the Eqs. (3.3), (3.4) and (3.5) can be simplified to Either A (+∞) = 0 implies i +∞ = i −∞ , or we divide the first equation by A (+∞) and solve the resulting linear system, which proves the claim.

Asymptotics around the fixed points
Let us recall that the Jacobian D (0,0,K ) of the ODE System (2.1) at a fixed point (0, 0, K ) has eigenvalues The existence of λ 0 implies the existence of a center manifold. In the present case, it locally coincides with the set of fixed points a = b = 0. This implies that there is no flow along the center manifold, so the asymptotics are fully described by the remaining two linear terms. The calculations are standard and presented in "Appendix A", along with a short review of the underlying theory. We only state the results here. First, regarding the unstable set S −∞ , as defined in (2.8): Locally, there exists a smooth unstable manifold of dimension one. Its restriction to {a ≥ 0} is the unique trajectory that emerges from the fixed point such that a(z), i(z) > 0 as z → −∞. It has the following properties: (4.2) Proof By choice of i −∞ > 1, the eigenvalue λ + is positive, whereas λ − is negative. Denote by u, v, w the coordinates in the system of eigenvectors e 0 , e + , e − of the Jacobian at the fixed, where the fixed point is shifted to the origin. This transformation is done explicitly in Lemma A.11. By Theorem A.12, the dynamics in an open neighborhood around the fixed point are equivalent to Hence, there is a stable and an unstable manifold, each of dimension one. The eigenvector e + describes the asymptotic direction of the unstable manifold, in coordinates a, b, i it is given by (4.4) Since λ − < 0 and λ + > 0, asymptotically along the branch of the unstable manifold in direction e + , where a > 0 and b > 0, it also holds that The choices of c change the type of convergence towards the origin: spiraling for c = 1, one stable manifold with eigenvalue −c/2, which has algebraic multiplicity 2 and geometric multiplicity 1 for c = 2, two stable manifolds for c = 3 Next, we prove Lyapunov-stability of the points in S +∞ , defined in (2.9). Figure 4 shows how the phase lines converge to (0, 0) in the (a, b)-plane. For technical reasons, we require that λ + = λ − . Later, we deal with this degenerate case via a continuity argument.

Theorem 4.2 (Stable set) For all c > 0 and i
Proof By choices of c and i +∞ , both non-zero eigenvalues (4.1) of the Jacobian are real-valued and negative and it holds that λ + = λ − . As before, denote by u, v, w the coordinates in the system of eigenvectors e 0 , e + , e − of the Jacobian at the fixed point, see Lemma A.11. By Theorem A.12, the dynamics of the system in a neighborhood around the fixed point are equivalent to (4.5) Take some small enough initial data ( u , v , w ): in view of Eq. (4.5), u does not vanish, but also does not propagate, whereas v and w converge to zero exponentially fast. Since a and b are represented in terms of v and w, see (A13), they vanish exponentially fast.

Proposition 4.3 Let c > 0 and i +∞ < i c . There is no non-negative and non-constant traveling wave that converges to
Proof As in the previous Theorem, the asymptotic behavior of Eq. (2.1) around the limiting fixed point (a, b, i) = (0, 0, i +∞ ) is described by the linear System (4.5). But now, since i +∞ < i c = max{0, 1 − c 2 /4}, either i +∞ < 0 or both eigenvalues λ ± have a non-vanishing imaginary part and thus, v and w spiral. Since a and b are represented in terms of v and w, see (A13), any trajectory that converges to (0, 0, i +∞ ) can not stay non-negative in its a-component.

Construction and result
We begin our search for a non-negative attractor of S +∞ in an easier setting: we fix i(z) = i = const. and investigate the two-dimensional sub-system in the remaining coordinates. To separate it from the full system, we write it asā(z),b(z). For this system, we prove the existence of a suitable attractor. This set will be denoted as T c (i), to emphasize that it depends on the chosen value of i, which will be constant only in this section. The flow of the sub-system and the region T c (i) are drawn in Fig. 5.
There are only two fixed points of (5.1), . We denote the eigenvalues and eigenvectors of the Jacobian at (0, 0) as Note that λ ± are identical to the non-zero eigenvalues of the full system around the fixed point (0, 0, i), see (4.1). The eigenvectors l ± are the projections of the corresponding three-dimensional eigenvectors into the (a, b)-plane. The Jacobian at (1 − i, 0) has eigenvalues and eigenvectors and it holds that β − (i) < 0 < β + (i). These have no direct correspondence to the three-dimensional system. We now define the region T c (i). It is a triangle, spanned by the two fixed points (0, 0) and (1 − i, 0) and two adjacent eigenvectors: Visually, it can easily be seen in Fig. 5 that the set T c (i) is invariant under Dynamics (5.1): the flow at the boundary of T c (i) points inwards.

Proposition 5.3 (Invariant region of the reduced system) The set T c (i) is an invariant region of Dynamics
converges to 0 monotonically as z → +∞.
For any non-negative solution of the full Wave System (2.1), it holds that i ≤ 0. Thus, we are interested in how T c (i) changes when i decreases: This proposition holds due to an easy geometric argument, again take a look at Fig. 5: when i decreases, the point (1 − i, 0) moves to the right and the two internal angles γ l (i) and γ r (i) increase. The computations for both propositions are performed in the following Sect. 5.2. The reader might skip those and proceed with Sect. 6, where we investigate the full system.

Invariance and monotonicity of T c (i)
We analyze the (ā,b)-system and the set T c (i) in detail. We prove that the Flow (5.1) at the boundary of T c (i) points inwards, and that the sets T c (i) are increasing in −i. We formalize what is sketched Fig. 5, and begin by examining the eigenvector l + (i) at the fixed point (ā,b) = (0, 0):

The quotient of the absolute values of theb-component andā-component of l
is the angle between the two vectors (1, 0) and −l + (i), the previous Lemma directly implies For the invariance of T c (i), we need The claim of the Lemma is now equivalent to First compute the Flow (5.1): .
(5.9) Its part in direction L inw is given by , we get similar results concerning its unstable eigenvector r + (i):

The quotient of the absolute values of theb-component andā-component of r
Now, since γ r (i) is the angle between the two vectors (−1, 0) and −r + (i), the previous Lemma implies Proof Let R inw := i − 1, −β − (i) be orthogonal to r + (i) and point inwards T c (i). The claim of the Lemma is now equivalent to ā b , R inw > 0.
Considering the invariance of T c (i), we conclude the

Attractor of the full system
We now analyze solutions of the full Wave System (2.1) under initial condition (a, b, i) = (a 0 , 0, i 0 ), such that (a 0 , 0) ∈ T c (i 0 ) as defined in the previous Section. In Sect. 6.1, we apply the results about the two-dimensional subsystem to the full system. Theorem 6.1 states that as long as i(z) ≥ i c , the (a, b)-components of the full system stay within the triangle T c (i c ). Thus, it suffices to control i(z) ≥ i c , which we do in two steps. In Sect. 6.2, we prove via some rough bounds that i(z) ≥ i c for sufficiently small initial values 0 ≤ a 0 1. This result is refined in Sect. 6.3: the Lyapunov-stability of the limiting point at z = +∞ implies that the entire trajectory including its limit is continuous in initial data. Carefully increasing a 0 , we increase the known attractor of the stable set S +∞ , resulting in Theorem 6.10. This procedure is sketched in Fig. 6.

Invariant region of the full system
Theorem 6.1 (Invariant region of the full system) Assume that i(z) ≥ i c for all z ∈ [0, ∞). We then can control the two remaining coordinates a(z), b(z) of the wave. It holds that Within T c (i c ), a ≥ 0 and b ≤ 0. Notice that while a, i ≥ 0, it holds that ci = −a(a + i + r ) ≤ 0. This directly implies the following Corollary 6.2 Under the assumption that i(z) ≥ i c for all z ∈ [0, ∞), the trajectory stays non-negative and converges as z → +∞: Proof of Theorem 6.1 In the full System (2.1) with coordinates (a, b, i), neither b nor b depend on i , but only on a and i. Thus, we can easily compare the full system to the two-dimensional System (5.1) in coordinatesā,b. At a phase-time z, the two vector fields (a , b ) and (ā ,b ) for fixed value i = i(z) are equal, compare (2.1) and (5.1). By Proposition 5.4, this implies that (a , b ) points strictly inwards T c i(z) . There are two irrelevant exceptions: for (a, b) = (0, 0), the system has already reached its limiting state. The point (a, b) = (1 − i 0 , 0) is a fixed point of the reduced, but not of the full system. In this case, since b = 0, b = a(a + i − 1) = 0, b = ai < 0 and ci < 0, a Taylor-expansion reveals that for small times > 0: i( ) changes much faster than a( ) and a( ) ∈ T c i( ) .
Importantly, the set T c i(z) is not decreasing as a function of z. Hence, at each phase-time z ≥ 0, the two components a(z), b(z) can not escape T c i(z) . In fact, since i ≤ 0 and T c (i) is increasing in −i, the set T c i(z) is increasing in z, at most up to T c (i c ). Thus, a(z), b(z) remain within T c (i c ) for all z ≥ 0.
With a similar argument, we can determine the rate of convergence: Proposition 6.3 Assume that i(z) ≥ i c for all z ∈ [0, ∞), such that (a, b, i) → (0, 0, i +∞ ) as z → +∞ for some i +∞ ∈ [i c , 1). If i +∞ > c 2 /4 − 1, then convergence is exponentially fast with rate Further, if i +∞ = c 2 /4 − 1, which can only happen if i +∞ = i c , then the system converges sub-exponentially fast. As z → +∞, the distance to the limit is of order z · e − c 2 z .

Proof
In the case i +∞ ∈ (i c , 1), all eigenvalues of the limit are simple, we refer to Sect. 4. The system converges exponentially fast, as shown in Theorem 4.2. It remains to determine the rate of convergence. The two candidates are λ ± = − c 2 ± c 2 4 + i +∞ − 1. Corresponding to λ ± , the projections of the eigenvectors into the (a, b)-plane are given by We know that a(z), b(z) ∈ T c (i +∞ ) for all z ≥ 0. At (0, 0), the triangle T c (i +∞ ) is bounded by the line −l + , see Definition 5.2 and Fig. 5. Since 0 > λ + > λ − , the direction of l − is steeper than that of l + , such that the line {q · l − | q ∈ R} lies outside T c (i +∞ ) for all q = 0. Thus, the two components a(z), b(z) cannot converge towards (0, 0) along −l − . But since they converge exponentially fast, the only possible remaining rate of convergence is λ + . We do not have a complete description of the asymptotics around the fixed point for the degenerate case λ + = λ − = −c/2. However, under the assumption that the system stays non-negative and converges, the relation (a, b, i) = − c 2 ·(a, b, i)+o(|(a, b, i)| 2 ) holds asymptotically. This will be proven in Theorem 6.10 (which does not rely on the type of convergence). The eigenvalue −c/2 has algebraic multiplicity 2, but geometric multiplicity 1. It is well-known that this results in sub-exponential convergence, cf. chapter 9 in Boyce et al. (2017).

A small attractor
In view of the previous paragraph, convergence and non-negativeness follow if we can show that i(z) ≥ i c for all z ≥ 0. If we choose a 0 small enough, some rough bounds do the trick. We control the total mass of active particles via Proof With the help of Theorem 6.1, we can use that a (6.6) By monotonicity: 1 − i 0 ≤ 1 − i(s) and 0 ≤ a(s) ≤ a 0 . It follows that where we need a 0 + i 0 < 1 to avoid a blow-up, which is true by our choice of a 0 . It remains to bound −b(z). Take a look at the flow in the (a, b)-plane in Fig. 5

. It holds that a(z), b(z) stay within T c (i c ). Within the triangle T c (i c ), it holds for the left inner angle
Since we can bound the total mass of active particles, we can also bound the change of i(z): Proposition 6.5 (Small attractor of S +∞ ) Fix c > 0 and i 0 ∈ (i c , 1). There exists a constant M(c, i 0 , r ), 0 < M 1, such that for all 0 ≤ a 0 ≤ M: Hence, also a(z), b(z) ∈ T c (i c ) for all z ≥ 0. The trajectory is non-negative and converges to S +∞ as z → +∞.
Proof As long as i(z) ≥ i c , it must be that a(z), b(z) ∈ T c (i c ) by Theorem 6.1. Assume there exists finite phase-time τ := inf z≥0 {i(z) < i c }: The upper bound a * is given in Definition 6.8. Trajectories with such initial data converge and stay non-negative, since where we used a(s) + i(s) ≤ 1. For z ≤ τ and a 0 sufficiently small, we can apply Lemma 6.4. This implies that there is a finite constant L, which does not depend on a 0 , such that The right-hand side is strictly larger than i c for sufficiently small a 0 , say a 0 ≤ M, and the bound is independent of the phase-time. Thus, there is no such τ for a 0 ≤ M.

Extending the attractor
The previous section ended with a condition of type a 0 1, under which the system stays non-negative and converges. However, given a 0 and i 0 and under the assumption that the system converges, we can explicitly calculate its limit i +∞ . Then, for fixed i 0 , we continuously deform the trajectory while increasing a 0 up to some upper bound a * (i 0 ), as sketched in Fig. 6. This results in Theorem 6.10.
We apply Proposition 3.2 to interrelate the limit (0, 0, i +∞ ) of the trajectory to its initial data (a 0 , 0, i 0 ): Lemma 6.6 If i(z) ≥ i c , such that the system stays non-negative and converges to (0, 0, i +∞ ) as z → +∞, then i +∞ can be written as a function of a 0 and i 0 : The function i +∞ (a 0 , i 0 ) is decreasing in a 0 , for a 0 ∈ [0, 1 − i 0 ].
Proof We apply Proposition 3.2, and solve the resulting system for i +∞ . In the present case, since a 0 = 0, this results in a quadratic equation with two possible solutions.
Since i(z) is decreasing, it must be that i +∞ < 1, which uniquely determines (6.12). A short computations proves that ∂ ∂a 0 i +∞ (a 0 , i 0 ) ≤ 0. We look for values of a 0 that ensure i +∞ ≥ i c . Thus, we rearrange (6.12) for a 0 , set i +∞ = i c , and choose the only possible positive solution of the resulting quadratic equation: Lemma 6.7 Given i 0 and under the assumption that i +∞ = i c , the value of a 0 is uniquely determined by (6.13) Equation (6.13) can be restated as i +∞ α(i 0 ), i 0 = i c , but keep in mind that still have to prove convergence. It can easily be seen that α(i c ) = 0. Since we require that a 0 ∈ [0, 1 − i 0 ], such that a 0 ∈ T c (i 0 ), this leads to our Definition 6.8 (Upper bound for a 0 ) For fixed c > 0 and, we define (6.14) This will hold as sharp upper bound for a 0 , such that the trajectory stays nonnegative and converges. Before we state the corresponding theorem, we perform a last check that we are in the correct setup: Lemma 6.9 Let i 0 ∈ [i c , 1) and a 0 ∈ [0, a * (i 0 )]. If the system stays non-negative and converges to (0, 0, i +∞ ), then where i +∞ (a 0 , i 0 ) is given as in Lemma 6.6.

Proof It holds that
After all these preparations, we are now ready to prove that the system converges and stays non-negative. We state Theorem 6.10 (Attractor of S +∞ ) For r ≥ 0, c > 0, let i 0 ∈ [i c , 1) and a 0 ∈ 0, a * (i 0 ) . Let a(z), b(z), i(z) be the solution of Eq. (2.1) with initial data (a 0 , 0, i 0 ).

Proof Notation:
We fix i 0 ∈ (i c , 1) and change only a 0 . If i 0 = i c , we must choose a 0 = 0. For a compact notation, z (x) is the state of the system at phase-time z, starting in x = (a, b, i). If the limit of a trajectory exists, we denote (6.17) Step 1: starting interval For a 0 positive but small enough, Proposition 6.5 grants that for all z ≥ 0: where T c (i c ) is a bounded invariant region that contains only points such that a ≥ 0, b ≤ 0. Thus, a(z), i(z) → (0, i +∞ ) monotone. With the help of Lemma 6.6, we can explicitly calculate i +∞ as stated in Eq. 6.16, and our claim holds on some small non-empty interval a 0 ∈ [0, a u ).
Step 2: neighborhood of existing trajectories Pick some a 0 < a * (i 0 ) for which the statement is already proven. By choice of a 0 , it holds that i +∞ > i c . Thus, the limit +∞ (a 0 , 0, i 0 ) is Lyapunov stable by our previous analysis of the asymptotics, see Theorem 4.2: for every ∞ > 0, there exists a δ ∞ > 0, such that for all z ∈ [0, ∞). Choose ∞ ≤ i +∞ (a 0 , i 0 ) − i c and assure that 0 < δ ∞ ≤ ∞ . This grants i(z) ≥ i c after entering the δ ∞ -neighborhood. Within this attractor, also a(z) ≥ 0 in view of Theorem 6.1, since i(z) ≥ i c . Starting in (a 0 , 0, i 0 ), we follow the trajectory up to some finite time τ , where it has entered the δ ∞ -neighborhood: (6.20) The derivative of the system is locally Lipschitz and all trajectories are within a bounded domain. Thus, the trajectories z (x 0 ) are uniformly continuous in initial data x 0 on any finite time interval [0, T ], with respect to the maximum norm ||.|| [0,T ] . This is a classic result and can easily be proven via a Grönwall's inequality, we refer to chapter 2 of the textbook of Hsieh and Sibuya (1999). There exists some δ 0 > 0, s.t. for all x ∈ R 3 with ||x − (a 0 , 0, i 0 )|| < δ 0 : This implies for all such trajectories z : (6.23) In particular, τ (x) lies within the δ ∞ -neighborhood, so ultimately Again, Theorem 6.1 implies a(z), i(z) ∈ T c (i c ) for all z ≥ 0. As before, the system is integrable and converges as z → +∞, so we can explicitly calculate i +∞ (a 0 , i 0 ). Thus, our claim holds for all starting points a small open neighborhood of (a 0 , 0, i 0 ).
Step 3: limits of trajectories Assume that the claim holds for all a 0 ∈ [0, a u ). For all trajectories starting in (a 0 , 0, i 0 ), where a 0 ∈ [0, a u ), it holds that i(z) is monotone and bounded from below by i c , such that a(z), b(z) stay within T c (i c ). Fix any finite time-horizon [0, T ]. As mentioned before, the trajectories z (a 0 , 0, i 0 ) are uniformly continuous in initial data on finite time-intervals, and thus form a Cauchy-sequence on ||.|| [0,T ] as a 0 → a u . Since T can be chosen arbitrarily large and since the limits +∞ (a 0 , 0, i 0 ) are also continuous in a 0 , our claim holds for the trajectory that starts in a u .
Step 4: conclusion By step 1, the claim holds for a 0 in some small interval [0, a u ). By step 3, it then also holds for a 0 = a u . If now a u < a * , the claim holds for a 0 ∈ [0, a u + ) by step 2 for some > 0. Iterating these two steps, the claim ultimately holds for all a 0 ∈ 0, a * (i 0 ) . In particular, we have proven that the trajectories z (a 0 , 0, i 0 ) are uniformly continuous with respect to initial data on z ∈ [0, +∞]. This continuity allows us to finish the proof of Proposition 6.3. In the non-critical cases where c 2 /4 + i +∞ −1 > 0, the trajectories converge along a stable manifold with rate of convergence −c/2 + c 2 /4 + i +∞ − 1. As c 2 /4 + i +∞ − 1 → 0, the critical trajectory must converge along the limit of these manifolds.

The complete trajectory
We track the non-negative branch of the unstable manifold of S −∞ , see (2.8), and show that it stays positive and enters the attractor of S +∞ from the previous section, cf. Theorem 6.10.
Assumption We will use the following setup over the entire Sect. 7: Let i −∞ > 1 and let a(z), b(z), i(z) be the unique solution of the ODE-System (2.1) that emerges from (0, 0, i −∞ ) as z → −∞ and where a(z) > 0 asymptotically as z → −∞.
For all i −∞ > 1, existence and uniqueness of these trajectories have been proven in Sect. 4. Moreover, we know their asymptotic behavior: Lemma 7.1 The following holds as z → −∞: The first two inequalities are given by Theorem 4.1, which also yields b (z) > 0 asymptotically. Noticing that c(a + i) = −a(1 + r ) − b < 0 completes the proof.

The maximum of active particles
For connecting these trajectories with the attractor of S +∞ , we will prove Proposition 7.2 (The maximum of active particles) There exists a finite phase-time z 0 , such that b(z 0 ) = 0 for the first time.
We will prove that the sum a(z) + i(z) decreases below 1. Given this, the term cb(z) + b (z) = a(z) · [a(z) + i(z) − 1] becomes negative, so b must eventually reach 0.
Lemma 7.4 As long as b(s) > 0 for all s ∈ (−∞, z], it can not happen that a(z)+i(z) converges to some finite L > 0.
Proof By the previous lemma: (a + i) < 0 while b > 0. Assume that a(z) + i(z) converges to a finite value L > 0 from above, which we denote as a(z) + i(z) L. This implies that also b +i 0. By the Wave Equations (1.2), these two expressions are equivalent to −a(1 + r ) − b 0 and (7.9) cb − a(a + i + r ) 0. (7.10) The first convergence indicates that b ≤ δ < 0 after some time z δ , since a is strictly increasing and hence positive. The second statement is equivalent to cb−a ·(L +r ) 0. Thus, also b is increasing. But b (z) < 0 for all z ≥ z δ and while b > 0, a contradiction.
We can now show that there exists a finite phase-time z 0 such that b(z 0 ) = 0, finishing the Proof of Proposition 7.2 By the previous lemma, a(z) + i(z) decreases below every positive value as long as b(z) > 0. In particular, for some > 0: a(τ ) + i(τ ) ≤ 1 − after some phase-time τ . Then for all z ≥ τ , since a > 0: Either cb(z) < 0 and the system has already passed a first local maximum of a(z), or we may assume that b (z) ≤ −a(τ ) = −δ < 0. If now b (z) ≤ −δ, then b(z) reaches zero after a finite time z 0 , which can not be larger than τ + b(τ ) δ .
If we assume additionally that is decreasing, which was proven in Lemma 7.3. Since a(z) is strictly increasing up to z 0 , i(z) must be strictly decreasing, but not below i(z 0 ) > 0. Hence, the trajectory stays positive. The inequality a z 0 +i z 0 ≤ 1 implies that a z 0 < 1.
This allows us to connect the unstable manifold of (0, 0, i −∞ ) with the attractor of S +∞ : Proposition 7.6 (Reaching the attractor of S +∞ ) Let i −∞ ∈ (1, 2 − i c ]. The nonnegative branch of the unstable manifold of (0, 0, i −∞ ) reaches the point (a z 0 , 0, i z 0 ), where a z 0 ∈ (0, 1) and i z 0 ∈ (i c , 1). It then holds that 0 < a z 0 ≤ a * (i z 0 ), (7.15) for a * like in Definition 6.8. In view of Theorem 6.10, the trajectory that starts/continues in such a point (a z 0 , 0, i z 0 ) converges to S +∞ and stays non-negative.
Proof We have just shown that i z 0 , a z 0 > 0 and that a z 0 +i z 0 ≤ 1. Recall Definition 6.8: (7.16) We have already verified that a z 0 ≤ 1 − i z 0 , so proving a z 0 ≤ α(i z 0 ) suffices for proving a z 0 ≤ a * (i z 0 ). By (7.13), we know that (7.17) The two expressions (7.16) and (7.17) are very similar. After some elementary steps, the claim a z 0 ≤ α(i z 0 ) is equivalent to This is equivalent to i −∞ ≤ 2 − i c , since i −∞ > 1 and i c ≤ 1. But that is just how we have chosen i −∞ .

Concluding the proof of the main result
With the results from the previous sections, we complete the For all z < z 0 , it holds that a(z), b(z), i(z) > 0. Denote the state of the system at z 0 as (a z 0 , 0, i z 0 ). Lemma 7.5 states that i z 0 ∈ (i c , 1), Proposition 7.6 states that a z 0 ∈ (0, a * (i z 0 )], for a * as in Definition 6.8. By Theorem 6.10, we then know that (a z 0 , 0, i z 0 ) lies in a non-negative attractor of the set S +∞ . Thus, For any non-negative and bounded solution, the identity i −∞ + i +∞ = 2 holds by Proposition 3.3. For c > 0 and i −∞ ∈ (1, 2 − i c ], the previous paragraph proves existence and uniqueness of the claimed wave. For i −∞ = 1, the constant solution can be the only non-negative and bounded one.
We then consider an arbitrary non-constant, bounded and non-negative solution. By monotonicity of i(z), it must converge as z → ±∞. If we assume that i −∞ ∈ (1, 2 − i c ], it is one of the above solutions. If we assume that i −∞ > 2 − i c , then i +∞ < i c . In this case, the trajectory can not stay non-negative as z → +∞, which is stated by Proposition 4.3, contradicting the assumption.

FKPP-waves
We have given a description of all bounded and non-negative traveling waves of the Reaction-Diffusion System (1.1). For the most related systems, the FKPP-equation (Kolmogorov et al. 1937;Fisher 1937), the FitzHugh- Nagumo-equation (FitzHugh 1961;Nagumo et al. 1962) and combustion equations (Berestycki et al. 1985), no such continuum of traveling waves has yet been constructed.
Still, the non-negative traveling waves of System (1.1) are closely related to pulled FKPP-waves with only a single type of particles. The equation for such a wave w(z) reads 0 = cw + w + F(w). For the purpose of a simple comparison, we let F(w) = gw − w 2 , where g > 0 is the initial growth rate of the particles. In this case, Theorem 1.1 states that the convergence of System (1.1) as z → +∞ is identical to that of w, if g = 1 − i +∞ , see e.g. (Uchiyama 1977). In words, the asymptotic growth speed of traveling waves of System (1.1) coincides with that of simple FKPP-waves in presence of a constant density i +∞ of inhibiting particles. Moreover, Theorem 1.1 implies that i c = 0 for all c ≥ 2. Thus, the minimal speed of an invasive front, where i +∞ = 0, is given by c min = 2. Again, this coincides with the minimal wave speed of the associated FKPP-equation, i.e. in the absence of inactive particles. It is this critical front which can be interpreted as the most natural one, our simulations indicate that it always arises under compact initial data. If we assume convergence, a technique of Berestycki et al. (2018) yields an upper bound for the speed of the traveling front, just by ignoring the dampening influence of the inactive particles. For compact initial data, the system always chooses the smallest possible wave speed, as suggested.
The emergence of traveling fronts is known for many reaction-diffusion systems. We suggest the literature (Britton 2003;Volpert and Petrovskii 2009;Othmer et al. 2009) for more examples with a biological motivation. Rigorous proofs of these phenomena are rare. Often, only the form of the traveling waves is analyzed analytically. The FKPP-equation is one of the cases, where the convergence of the front of the PDE towards a traveling wave solution can be proved. The first rigorous proof was done by Kolmogorov, Petrovsky & Piscunov in 1937Kolmogorov et al. (1937. Extensions of this result to more general initial data and a more precise description of the speed of the front have been provided by Uchiyama (1977) and Bramson (1983). The approach of Kolmogorov et al. and Uchiyama seems to be restricted to systems with only a single type of particles, as it relies on a maximum principle and monotonicity of the front. The approach of Bramson relies on a relationship between the FKPP-equation and branched Brownian motion, which can not be applied in the present case since the inactive particles do not diffuse. A singular perturbation of System (1.1), which again introduces a small diffusion to the inactive particles will be subject to future investigations. Despite the fact that this system would be biologically interesting, since no tissue or population is entirely static, this would also rule out some difficulties when analyzing the stability of the traveling waves against perturbations.

Stability of the traveling waves
We give a brief introduction to the stability of traveling waves against small perturbations, in the spirit of the introduction in Ghazaryan et al. (2013). A good overview, where the following concepts are presented in greater depth, has been written by Sandstede (2002). . . . , d n ) with d i ≥ 0, and R a smooth reaction. In the moving frame z = x − ct, the System reads

Consider a reaction-diffusion system
A traveling wave w(z) with speed c is a constant solution of Eq. (9.2). The wave w is called non-linearly stable in a space X , if any solution of the PDE (9.2) which starts in Y 0 = w +Ỹ , whereỸ ∈ X is a sufficiently small perturbation, converges to a shift of w. This type of stability is often encoded in the spectrum of the operator L, that is obtained by linearizing the equation for the perturbationỸ in (9.2) around to the constant part w:Ỹ where J R is the Jacobian of the reaction R. Let L : X → X be the operator given byỸ → LỸ . We say that the wave w is spectrally stable in X if the spectrum of L is contained in the half-plane Re(γ ) < 0, except maybe a simple eigenvalue at 0 (that corresponds to the traveling wave itself, if w ∈ X ). For diffusive systems, where det(D) > 0, a quite general theory has been developed. If X is appropriately chosen, spectral stability implies non-linear stability, we refer to the literature (Sandstede 2002;Ghazaryan et al. 2013). Classical results are e.g. given for subspaces of X = H 1 , the L 2 -Sobolev space. As explained in the next paragraphs, we are not aware of any rigorous framework for studying the non-linear stability of System (1.1). Several technical problems arise, that so far have been treated only separately (Ghazaryan et al. 2013;Kirchgässner 1992).
Most importantly, the traveling waves of System (1.1) can not be stable against perturbations in the classical sense, since the inactive particles neither react nor diffuse. Any initial deviation remains for all times, as shown in Fig. 1. However, the actual front of the system does converge to a traveling wave. For capturing this idea, we introduce the weighted space X = H 1 α with norm || f || H 1 α = || f · e αz || H 1 for some α > 0. Non-linear stability in H α 1 is referred to as convective stability. Convergence of the PDE in the moving frame (9.2) in H 1 α means that the front of the system approaches the traveling wave, whereas any initial finite and local deviation is convected towards z = −∞ and vanishes due to the weighting. A first rigorous result regarding convective stability was obtained by Ghazaryan et al. (2013). They could show that in some cases, spectral stability in H 1 α implies convective stability against small perturbations in H 1 α ∩ H 1 . For their approach, the authors require that the weight α can be chosen such that all eigenvalues γ of L except zero fulfill Re(γ ) ≤ ν < 0 and such that the derivative w ∈ H 1 α of the traveling wave is an eigenfunction that corresponds to a simple eigenvalue at zero. Unfortunately, this setting is not suited for studying pulled FKPP-fronts: the assumption w ∈ H 1 α implies that the continuous spectrum of L touches the origin, see e.g. chapter 6 in the work of Sattinger (1976).
Another difficulty arises when studying critical pulled fronts (with minimal possible speed) whose tail as z → +∞ converges sub-exponentially, as in Theorem 1.1. In this case, the requirement w ∈ H 1 α is only fulfilled for rather small values of α, which do not suffice for shifting the continuous spectrum of L to the left half-plane. For diffusive systems, non-linear stability of this more delicate case was first treated rigorously by Kirchgässner (1992), a recent overview is given in Faye and Holzer (2018). After introducing a small diffusion to the inactive particles, we could apply this theory to the critical front.
For the most natural traveling wave solution, the critical one with speed c = 2 and i +∞ = i c = 0, we performed a numerical analysis that strongly indicates that this wave is spectrally stable in H 1 α , when we choose α = −μ +∞ = c/2. The details are presented in "Appendix B". Thus, based on the work of Ghazaryan et al. regarding convective stability (Ghazaryan et al. 2013) and the work of Kirchgässner regarding critical fronts (Kirchgässner 1992), we dare to make an educated guess: we expect that this traveling wave is convectively stable against small perturbations in H 1 c/2 ∩ H 1 .
such that dx/dt = f (x) is equivalent to: We require that the eigenvalues of A ∈ R k×k have zero real parts and those of B ∈ R l×l have nonzero real parts. Further, both functions g : R n → R k and h : R n → R l are smooth and vanish together with their first-order partial derivatives at the origin. Then, System (A1) is called the normal form of the system.
Proposition A.2 Let f : R n → R n be smooth. Let a dynamical system d x/dt = f (x), x ∈ R n have a fixed point x 0 ∈ R n , such that the eigenvectors of the Jacobian D f (x 0 ) span the entire R n . The system can be written in normal form as in Definition A.1.
The proof includes a change of coordinates into the system of eigenvectors of the Jacobian D f (x 0 ). This will be done explicitly in Section A.2. For the underlying theory and the more general case, we refer to the monograph of Kirchgraber and Palmer (1990).

Definition A.3 (Center manifold) Consider a dynamical system in normal form (A1).
Let φ : R k → R l be a smooth function such that φ(0) = 0 and also its derivative Dφ(0) = 0. Assume that the set is invariant under Dynamics (A1). The set CM is called a center manifold of the fixed point (due to its vanishing derivative at 0).
We will use a local version of the center manifold, which can be shown to exist in a neighborhood of the fixed point: Theorem A.4 (Local center manifold, cf. Theorem 4.1. in Kirchgraber and Palmer (1990)) Consider a smooth dynamical system in normal-form (A1), where dim(y) = k ≥ 1, such that the Jacobian at the fixed point has k eigenvalues with zero real part. Let c 1 + c 2 = dim(z), where the matrix B has c 1 eigenvalues with positive real part and c 2 eigenvalues with negative real part. Then locally, there exist a unique center manifold of dimension k, a unique unstable manifold of dimension c 1 and a unique stable manifold of dimension c 2 .
The center manifold can be written as (y, z) : z = φ(y) like in (A2). There exists a homeomorphism defined in an open neighborhood of the origin which takes solutions of d x/dt = f (x) onto solutions of dy dt = A · y + g y, φ(y) , Definition A.5 (Error of approximation of the center manifold) Consider a smooth dynamical system in normal-form (A1). For a smooth function T : R k → R l define the error of approximation of the normal form by Theorem A.6 (Approximating the center manifold, cf. Theorem 3 in Carr (1982)) Consider a smooth dynamical system in normal form (A1) with local center manifold {(y, z) : z = φ(y)} as in (A2). Let T : R k → R l be smooth with T (0) = 0 and DT (0) = 0. Suppose that as y → 0, for some q > 1: Then, as y → 0, also

A.2 Calculating the normal form and the center manifold
We analyze the flow of the ODE System (2.1) around its fixed points by applying the theory from the previous section. We therefore write the system into normal form, see Definition A.1. For a fixed point (a, b, i) = (0, 0, K ), we begin with the affine transformation and then decompose the resulting system into a linear part M and a non-linear part G.
To be concise with the notation from the previous section, which is adopted from the existing literature, we use a vectorial notation in coordinates ( j, a, b), such that the center manifold can be written as ( j, a, b) : (a, b) = φ( j) .
Definition A.7 Given c > 0, K ∈ R, introduce the matrix M as Further, define the non-linear functions g( j, a) := a 2 + a j and G : R 3 → R 3 : Definition A.9 For given c > 0, K ∈ R, we define the discriminant The eigenvalues and eigenvectors of M are then given by (cf. (2.5)) Technical difficulties arise when the eigenvectors no longer span the entire R 3 . We require K = 1 in view of (A13), which also eliminates the case λ + = 0. For similar reasons, we also exclude the case that λ + = λ − , so we require that = 0. This is Lemma A.11 (Dynamics in normal form) Let c > 0 and K / ∈ {1, 1−c 2 /4}. The eigenvectors e 0 , e + , e − of M form a basis of R 3 . We introduce the coordinates (u, v, w), such that any x ∈ R 3 can be written as x = u · e 0 + v · e + + w · e − . The System (2.1) in coordinates (u, v, w) follows dynamics where P is a polynomial such that P(u, 0, 0) = 0: Proof We change coordinates from u, v, w to j, a, b and back: Now recall the functions G and g, see (A9). Explicitely calculating the non-linear part Luckily, for evaluating g( j, a), we do not have to calculate the coordinate b (u, v, w).
Now we have all ingredients for computing the center manifold. We use the approximation argument from Theorem A.6.
α > 0 and analyze the spectral stability of the critical traveling wave in the weighted Fig. 7 Numerical evaluation of the Evans-function (B6) Ev(γ ) for r = 0 and γ on the boundary of the Domain S, defined in (B7). The Evans-function for r = 1 is very similar. The origin is marked with a small red cross. The graph does not enclose the origin and it can visually be seen that its winding number is equal to zero. We conclude that the Region S contains no zeros of Ev(γ ) We investigate the case α = c 2 , which is equal to the rate of convergence of the wave as z → +∞, up to a sub-exponential term, see Theorem 1.1. For α = c 2 , then within the region {Re(γ ) ≥ 0, γ = 0} the following holds: the dimension of the unstable space of M(−∞, γ ) + c 2 · 1 is given by k − = 2, and the dimension of the stable space of M(+∞, γ ) + c 2 · 1 is given by k + = 1. This can easily be deduced from the corresponding Eigenvalues (B3), (B4), which do not cross the imaginary axis. The values k − and k + add up to the dimension of the ODE (B2). We say that {Re(γ ) ≥ 0, γ = 0} is contained in the region of consistent splitting of the operator L. This implies that the non-negative part of the spectrum of L is contained in the point spectrum of the operator, which is a standard result (Sandstede 2002;Sattinger 1976). Within the region of consistent splitting, we can define the Evans-function E(γ ).
Given γ with Re(γ ) ≥ 0, γ = 0, we let X (z) be the unique (up to a shift) solution of Eq. (B5) that vanishes at z = +∞, and let Y 1 (z), Y 2 (z) span the two-dimensional space of solutions of Eq. (B5) that vanish at z = −∞. The Evans-function is defined as (B6) The Evans-function is not unique, but it holds that E(γ ) = 0 if and only if γ lies in the point spectrum of L. Moreover, E(γ ) is analytic if X , Y 1 , Y 2 are chosen such that they are analytic in γ (Sandstede 2002). Thus it suffices to calculate E(γ ) along the boundary of a domain: the winding number along this contour then corresponds to the number of zeros inside the domain. We use this to verify that there are no zeros of E(γ ) within the set S := γ ∈ C Re(γ ) ≥ 0, 10 −3 ≤ |γ | ≤ 1000 , where we keep a small distance from the origin to prevent that the eigenvalues of the matrices M(±∞, γ ) + c 2 · 1 touch the imaginary axis, and hope that there are no unexpectedly large eigenvalues. We want to remark that for systems with a degenerate diffusion, no general a priori upper bound for the size of the eigenvalues with non-negative real-part has been found yet, which would allow for a numerical proof of spectral stability. It may be possible to generalize the approach in Lattanzio and Zhelyazov (2021). Simple energy estimates exist for traveling waves of diffusive systems, see e.g. chapter 6 in Ozbag and Schecter (2018).
The various numerical challenges that arise when computing the Evans-function, as well as their solutions, are described in detail by Barker et al. (2018), who also suggest using their library STABLAB Barker et al. (2015). We gratefully followed this suggestion, and computed the left-adjoint Evans-function, a slight modification which is numerically advantageous in the present setting (Barker et al. 2018). The result is presented in Fig. 7 and yields a strong evidence that the critical wave is spectrally stable in H 1 c/2 . If c 2 4D + r S · i +∞ − r A = 0, convergence as z → +∞ is sub-exponentially fast and of order z · e − c 2D z . If c 2 4D + r S · i +∞ − r A > 0, convergence as z → +∞ is exponentially fast. Convergence as z → −∞ is exponentially fast in all cases. The corresponding rates are In particular, for any invading front, where i(z) → 0 as z → +∞, the remaining density of particles at the back of the wave is given by i −∞ = 2 · r A r S .