Shock-fronted travelling waves in a reaction-diffusion model with nonlinear forward-backward-forward diffusion

Reaction-diffusion equations (RDEs) are often derived as continuum limits of lattice-based discrete models. Recently, a discrete model which allows the rates of movement, proliferation and death to depend upon whether the agents are isolated has been proposed, and this approach gives various RDEs where the diffusion term is convex and can become negative (Johnston et al., Sci. Rep. 7, 2017), i.e. forward-backward-forward diffusion. Numerical simulations suggest these RDEs support shock-fronted travelling waves when the reaction term includes an Allee effect. In this work we formalise these preliminary numerical observations by analysing the shock-fronted travelling waves through embedding the RDE into a larger class of higher order partial differential equations (PDEs). Subsequently, we use geometric singular perturbation theory to study this larger class of equations and prove the existence of these shock-fronted travelling waves. Most notable, we show that different embeddings yield shock-fronted travelling waves with different properties.


Introduction
Reaction-diffusion equations (RDEs) are widely used to study population dynamics in cell biology and ecology [23]. Often, U (x, t) represents a  [16]. Pink discs represent isolated agents and blue discs represent grouped agents. During each time step of duration τ , isolated agents attempt to move to nearest neighbour lattice sites with a probability P i m , to proliferate to form new agents in neighbour sites with a probability P i p and to die with a probability P i d . Similarly, grouped agents attempt to move to neighbour sites with a probability P g m , to proliferate to form new agents in neighbour sites with a probability P g p and to die with a probability P g d . The attempts that would move to an occupied site or place an agent on an occupied site are aborted.
population density and provides a macroscopic description of individual behaviour. For RDEs established from the continuum limit of stochastic models, a solution of the RDE not only shows the macroscopic evolution of U (x, t), but it also reflects how microscopic behaviour of individuals influences the macroscopic outcomes [1,6,16,17,21,33]. Johnston et al. [16] introduced a lattice-based stochastic model to study how a population of individuals can undergo motility, proliferation and death events with the aim of studying biological and ecological invasion, see Figure 1. By considering different behaviours of isolated and grouped agents, including motility, proliferation and death events, an RDE with a nonlinear diffusivity function and a logistic or Allee type reaction term was derived as the continuum limit. In particular, where U (x, t) represents the total population density at position x ∈ R and time t ∈ R + .
The nonlinear diffusivity function is given by where D i ≥ 0 and D g ≥ 0 are diffusivities of the isolated and grouped agents, respectively. When D i > 4D g , D(U ) has two real roots, α and β, which are centred around 2/3, and are given by and D(U ) < 0 for U ∈ (α, β). While the negativity of a nonlinear diffusivity function is sometimes related to aggregation in the underlying discrete model [34], here it is actually a macroscopic effect of the isolated and the grouped motility of the agents, together with competition for space, that leads to a net aggregation effect [16]. The condition D i > 4D g implies that the motility rate of isolated agents is greater than the motility rate of grouped agents, which is consistent with the common biological observation that isolated leader cells are more motile than follower cells [28,32]. Note that D i and D g are related to P i m and P g m , respectively, in the lattice-based model in Figure 1. Full details of the discrete model and the continuum limit derivation are given in [16].
The reaction term, whose parameters are also related to parameters in the lattice-based model depicted in Figure 1, is given by where λ i ≥ 0 and λ g ≥ 0 are the proliferation rates of isolated and grouped agents, respectively; K i ≥ 0 and K g ≥ 0 are the death rates of isolated and grouped agents, respectively [16]. If the proliferation mechanism is the same for isolated and grouped agents and no death event occurs, that is, λ i = λ g and K i = K g = 0, then (4) simplifies to a logistic reaction term If the proliferation and death mechanisms are either competitive or cooperative, that is, λ i = λ g and K i = K g [35,38], then the reaction term takes the form of an Allee effect [38]. For simplicity, but without loss of generality, Figure 2: (a) The nonlinear diffusivity function D(U ) (2) centred around 2/3 (dashed line). (b) The reaction term R(u) corresponding to the logistic growth with λ g = 1 (red), the weak Allee effect with r = 0.8 and A = −0.9 (blue), and the strong Allee effect with r = 3 and A = 0.3 (green).
we assume K g = 0. 1 Subsequently, (4) simplifies to where r = K i − λ i + λ g is the intrinsic growth rate and A = 1 − λ g /r is the Allee parameter [16]. If r > λ g , which is equivalent to K i > λ i and thus implies that isolated agents have a higher death rate than proliferation rate, then 0 < A < 1 and R(U ) < 0 in (0, A) and R(U ) > 0 in (A, 1). This represents the strong Allee effect. Conversely, if 0 < r < λ g , which implies that isolated agents have a higher birth rate than death rate, then A < 0 and R(U ) > 0 in (0, 1). This is called the weak Allee effect [16]. See Figure 2 for the different potential forms of R(U ). For simplicity, we assume that A = α and A = β. Understanding travelling wave solutions, that is, solutions that propagate through space with a fixed shape and a constant speed, is important in the study of biological and ecological invasion processes [? 9, 14, 31]. In this work, we are interested in travelling wave solutions supported by (1) with D i > 4D g , such that the nonlinear diffusivity function is negative for U ∈ (α, β), see (3).
In this case, the nonlinear diffusivity function can be written as With the implicit finite difference method introduced in [16], numerical solutions of (1) with D(U ) as in (7) and with either logistic or weak Allee forms for R(U ) lead to smooth travelling wave solutions with positive speeds, while simulations of (1) with D(U ) as in (7) and strong Allee forms for R(U ) lead to shock-fronted travelling wave solutions with either positive or negative speeds [16,21], see Figure 3 for different travelling wave solutions at t = t 1 , t 2 , t 3 , with t 1 < t 2 < t 3 . To calculate the wave speed, we locate the front of the wave by looking for the left-most coordinate x l satisfying U (x l , t) < 10 −3 . Then, we estimate the speed from the distance the front of the wave has travelled from t 2 to t 3 . Interestingly, the speeds of the shock-fronted travelling wave solutions are much smaller than the speeds of the smooth travelling wave solutions, which potentially indicates that the mechanisms giving rise to shock-fronted travelling waves are fundamentally different to the mechanisms that give rise to smooth travelling waves. Note that with the nonlinear diffusivity function D(U ) centred around 2/3 given by (2) we only observe shock-fronted travelling wave solutions with the strong Allee effect. However, a forward-backward-forward nonlinear diffusivity function which is not centred around 2/3 may also lead to shock-fronted travelling wave solutions with logistic or weak Allee forms of R(U ), see Figures 9 and 10 in [21] for an example. Ferracuti et al. [12] showed that there exist smooth travelling wave solutions of (1) with logistic R(U ) for a range of positive wave speeds based on the comparison method [2]. Kuzmin and Ruggerini [20] provided necessary conditions for the existence of smooth travelling wave solutions of (1) with R(U ) that takes the form of an Allee effect and the speed of the wave can be either negative or positive according to the shape of D(U ) and R(U ). However, to the best of our knowledge, the existence of shock-fronted travelling wave solutions to (1) with D(U ) (7) and R(U ) taking the form of an Allee effect is an open question. The methods used in [20] can also be used to identify necessary conditions for the existence of shock-fronted travelling wave solutions. In particular, let U 1 and U 2 (with U 1 ≤ α < β ≤ U 2 ) denote the U -values at the endpoints of the shock, then a necessary condition for the existence of a monotonically decreasing shock-fronted left-travelling wave solution is for all U a ∈ (0, U 1 ). Similarly, a necessary condition for the existence of a monotonically decreasing shock-fronted right-travelling wave solution is for all U b ∈ (U 2 , 1). We refer to Appendix A for a derivation, inspired by [20], of these necessary conditions. Heuristically this means that for shock-fronted travelling wave solutions with a positive speed, the A value in R(U ), in order to satisfy (9), should not be too close to 1. Since A = 1 − λ g /r where r = K i − λ i + λ g , this implies that a relatively much higher death rate of isolated agents compared to the birth rate of isolated agents will not result in a successful invasion event. Similarly, a very small birth rate of grouped agents will also not result in a successful invasion event.
In [21], we derived, among other things, the same condition as in [12] for the existence of smooth travelling wave solutions of (1) with logistic R(U ) by using a geometric approach. Furthermore, geometric approaches have been used to study shock-fronted travelling wave solutions. For example, in [14,31], the authors studied shock-fronted travelling wave solutions in an advection-reaction-diffusion equation for malignant tumour invasion using geometric singular perturbation theory (GSPT) [11,15,18] and canard theory [37,39,40]. In this work, we use GSPT to further explore the existence of shock-fronted travelling wave solutions of (1) with lim x→−∞ U (x, t) = 1 and lim x→∞ U (x, t) = 0. 2 Moreover, we assume D i > 4D g -such that D(U ) < 0 for U ∈ (α, β) -and K g = 0 and r > λ g -such that we have a strong Allee effect type R(U ). To apply GSPT, we smooth out the shock and regularise (1) by adding a small higher order perturbation term. This embeds (1) into a 2 U ≡ 0 and U ≡ 1 are both constant solutions of (1) with the logistic R(U ) and the weak and strong Allee type R(U ). However, U ≡ 0 and U ≡ 1 are both PDE stable for the strong Allee type R(U ), while only U ≡ 1 is PDE stable for the logistic R(U ) and the weak Allee type R(U ). Therefore, it is no surprise that left-moving traveling wave solutions are only found for the strong Allee type R(U ), see Figure 3.
larger class of PDEs. Regularisation of RDEs is typically considered in one of two ways [25,26]. The first method of regularisation accounts for non-local effects by adding a small fourth-order spatial derivative term [26,41]. In particular, equation (1) becomes The second method of regularisation accounts for viscous relaxation by adding a small mixed derivative term [24,25,42]. In particular, equation (1) becomes It is important that the sign of the perturbation terms in (10) and (11) is such that setting ε > 0 generally leads to well-posed problems. However, see [21] and references therein, for a further discussion related to the well-posedness of (1). Also note that other types of regularisations have been used to smooth out shocks [4].
In §2, we study travelling wave solutions of (10) and first derive a higherdimensional slow-fast system of ordinary differential equations (ODEs). The related reduced singular limit ODE systems give useful information of underlying shock-fronted travelling wave solutions of (1) and (10) based on GSPT and Fenichel theory [11]. Because the reduced systems are algebraically intractable, we use a numerical ODE solver to determine the speed of the shock-fronted travelling wave solutions. In §3, we use a similar approach to establish a different higher-dimensional system of ODEs based on the viscous relaxation PDE (11) and find shock-fronted travelling wave solutions with different properties. Note that in this case, GSPT has to be extended since the critical manifold loses normal hyperbolicity near a fold point. Although (10) and (11) are the same in the singular limit ε = 0, they yield shock-fronted travelling wave solutions with different speeds and different shock sizes when ε > 0. Finally, we discuss various extensions of the current work including the relationship between the discrete model and the continuous description, the option of including different regularisation terms, the possibility of shockfronted travelling wave solutions with logistic R(U ) and the spectral stability of travelling wave solutions of (1).

Remark 1.
In the remainder of this article we will use nonlinear diffusivity functions D(U ) (2) and reaction terms R(U ) (4) that are larger than the D(U ) and R(U ) used in Figure 3 to generate larger speeds. As the model based on (1) is dimensionless, those larger parameters in D(U ) and R(U ) still correspond to the parameters introduced in the latticed-based model in Figure 1 upon rescaling space and/or time. However, note that the connection between the discrete and continuum models is only accurate when the rate of motility of both the grouped and isolated agents is much greater than rate of proliferation and death of both the grouped and isolated agents and this should be kept in mind when rescaling space and/or time. For more details, see [16].

Non-local regularisation
In this section, we look for shock-fronted travelling wave solutions of (10) connecting U = 1 to U = 0. We first introduce a travelling wave coordinate to transform (10) into a fourth-order ODE. Next, we use a dynamical system approach to transform the ODE into a four-dimensional singular perturbed slow-fast system. The four-dimensional system has two equivalent forms as ε = 0. However, these forms produce different lower-dimensional subsystems, called the reduced problem and the layer problem in the singular limit ε = 0. 3 The concatenation of solutions of each of the subsystems yields a solution of the four-dimensional system in the singular limit. We give an outline, and conclude based on GSPT, that it persists for ε sufficiently small in the full four-dimensional system. This solution corresponds to a travelling wave solution of (10).

Preliminary observations
A travelling wave solution of (10) is a solution of the form U (x, t) = u(x − ct) = u(z), where c ∈ R is the constant speed of the travelling wave solution and z = x − ct is the travelling wave coordinate. Writing (10) in its travelling wave coordinate leads to where F (U ) = D(U )dU and the reaction term R(U ) (6) is of strong Allee effect type. A travelling wave solution u(z) is a stationary solution to (12) that asymptotes to one as z → −∞ and to zero as z → ∞. Thus, it satisfies Upon defining (13) transforms into a four-dimensional singular perturbed slow-fast dynamical system Here, (u, w) ∈ R 2 are fast variables and (p, v) ∈ R 2 are slow variables. By using a stretched, or fast variable, ξ = z/ε [11], (15) is transformed into an equivalent fast system, provided ε = 0, The three fixed points 4 of the two equivalent systems (15) and (16) are and we are interested in heteroclinic orbits connecting P 1 ε with P 0 ε as these correspond to travelling wave solutions of (10) that asymptote to 1 as x → −∞ and to 0 as x → ∞. Note that due to the symmetry (w, p, z, c) → (−w, −p, −z, −c) of system (15), the existence of a heteroclinic orbit connecting P 1 ε with P 0 ε also implies the existence of a heteroclinic orbit connecting P 0 ε with P 1 ε and this latter orbit corresponds to a travelling wave solution of (10) that asymptotes to 0 as x → −∞ and to 1 as x → ∞ and moves in the opposite direction.
The characteristic equation of the Jacobian of (16) is given by where we used that F (u) = D(u) and observe that u = 0, 1 or A at a fixed point. Upon substituting a regular expansion τ = τ 0 + ετ 1 + O(ε 2 ) into (18), we obtain an expansion for the four eigenvalues of the Jacobian That is, both the stable and unstable manifolds P 0,1 ε are two-dimensional. At P A ε , the stable and unstable manifolds depend on the sign of c. If c > 0, the stable manifold of P A ε is three-dimensional and the unstable manifold of P A ε is one-dimensional, while the situation for the stable and unstable manifolds of P A ε is the opposite for c < 0. For c = 0, that is, for a standing wave, we again have that the stable and unstable manifold of P A ε are two-dimensional. While the slow system (15) and the fast system (16) are equivalent for ε = 0, they have different singular limits. The singular limit of the fast system, that is, the layer problem, describes the dynamics near the shock and the fast variables (u, w) will change significantly here while the slow variables (p, v) are to leading order constant. In contrast, the singular limit of the slow system, that is, the reduced problem, describes the dynamics away from the shock and here the fast variables will be slaved to the slow variables.

Layer problem
The layer problem is obtained by letting ε → 0 in the fast system (16), which gives as well as dp/dξ = 0 and dv/dξ = 0, that is, (p, v) ∈ R 2 are constants in (20). The union of fixed points of (20) forms a two-dimensional invariant manifold, which is the so-called critical manifold [18], see Figure 4. The Jacobian of (20) is Therefore, the manifold M 0 loses normal hyperbolicity when D(u) ≤ 0, that is, for u ∈ [α, β] the eigenvalues τ ± (22) are purely imaginary. As such, we split the critical manifold M 0 into two two-dimensional normally hyperbolic saddle-type branches Figure 4: (a) A projection of the four-dimensional phase plane of (15) and the critical manifold M 0 . A shock-fronted travelling wave solution of (10) as shown in (b) corresponds to a heteroclinic orbit (indicated in blue in (a)) that starts at P 1 ε on the normally hyperbolic branch M − 0 of M 0 and that follows the dynamics of the reduced problem (RP), whose projection on the (u, p)-plane is shown in (c), before it jumps to the other normally hyperbolic branch M + 0 according to the dynamics of the layer problem (LP). The projection of the layer dynamics on the (u, w)-plane, since p and v are constant, is shown in (d) and the two blue curves connecting u − and u + in (a) correspond to the heteroclinic orbits u 0,± h , w 0,± h in (d). On M + 0 , the heteroclinic orbit again follows the dynamics of the reduced problem and asymptotes to P 0 ε . A shock-fronted travelling wave solution is thus composed by orbits in the reduced problem and the layer problem as indicated in (b). a two-dimensional not normally hyperbolic centre-type branche and the two one-dimensional boundary sets The layer problem (20) describes the dynamics near the shock away from the critical manifold. It is a Hamiltonian system and, as such, we are looking for a heteroclinic orbit connecting M − 0 with M + 0 . The Hamiltonian of (20) is given by where G(u) = F (u) du. Any solution is confined to a level set of the Hamiltonian and we have that where u ± ∈ M ± 0 are the end-points of the heteroclinic orbit such that 0 < u + < α < β < u − < 1. This is equivalent to the integral equation which is the well-known equal area rule, see, for example, [41]. Recall that F (u) = D(u) du and F (u) thus has an integration constant. Therefore, for a specific F (u) the value of v satisfying the equal area rule (23) is unique.
In Appendix B we show that (20) supports two heteroclinic orbits connecting (u + , 0) and (u − , 0) and these heteroclinic orbits u 0,± h , w 0,± h are given by where a = (D i − D g )/2 and we recall that u − > u + by construction. See Figure 5.

Reduced problem
The reduced problem is obtained from (15) by letting ε → 0, which gives and the two algebraic constraints w = 0 and v + F (u) = 0. Hence, (24) simplifies to Morever, since w = 0 and F (u) = −v, (25) governs the flow on the critical manifold M 0 . The reduced problem is singular along the two lines u = α and u = β since D(α) = D(β) = 0. Therefore, we transform (25) into a desingularised system 5 by using a stretched variable dψ = dz/D(u) [3,21] It is important to note that, while the stretching changes the speed along a trajectory in a nonlinear fashion, the trajectories of the phase portraits of the reduced problem (25) and the desingularised problem (26) are the same. However, the orientation along a trajectory is reversed for u ∈ (α, β) as D(u) < 0. System (26)  However, the desingularised system is more amenable to analysis and we thus study the dynamics of this desingularised system.

The construction of the heteroclinic orbit in the singular limit
Since the fixed points P 0,1 ε (17) are on the normally hyperbolic branches M ± 0 of the critical manifold M 0 , a shock-fronted travelling wave solution to (10) corresponds to a heteroclinic orbit of (15) that, to leading order, starts on M − 0 , follows the dynamics of the reduced problem (25) before it jumps, according to the layer dynamics (20), to the other normally hyperbolic branch M + 0 on which it asymptotes to P 0 ε following the dynamics of (25) again. In particular, if we split the spatial domain z ∈ (−∞, ∞) into three parts: then the heteroclinic orbit is, to leading order, on M ± 0 and governed by the reduced problem (25) for z ∈ I ± s , while it is, to leading order, governed by the layer problem (20) for z ∈ I f , see Figure 4. Note that due to translation invariance of (12) we can, without loss of generality, set z * = 0 in (27).
Since w = 0 and F (u) = −v on the critical manifold, the fixed points (0, 0) and (1, −c) of the reduced problem (25) correspond to P 0 ε and P 1 ε , respectively. Furthermore, the analysis of the layer problem (20) -which is independent of the speed c -indicates there may exist shocks with endpoints u − (> β) and u + (< α). Consequently, if there exists a shock-fronted travelling wave solution of (10) with a shock from u − to u + , it relates to two trajectories in system (25), see also Figure 4. These, in turn, relate to two corresponding trajectories in the desingularised system (26). One is the unique trajectory γ + , for a given speed c, that starts on the line {(u + , p + ), p + ∈ R} and approaches (0, 0) as ψ → ∞, while the other one is the unique trajectory γ − that arrives at the line {(u − , p − ), p − ∈ R} and approaches (1, −c) as ψ → −∞. Note that these unique trajectories can intersect the lines {(u ± , p ± ), p ± ∈ R} multiple times, see, for instance, Figure 6a. However, only the first intersections may lead to monotone travelling wave solutions. We are mainly interested in monotone travelling wave solutions since nonmonotonic travelling wave solutions are often PDE unstable [? ]. Therefore, we only look for these first intersections. As p is a slow variable, it should, to leading order, hold constant at the endpoints of the shock (dp/dξ = 0 in the singular limit). Hence, we are interested in the speeds c 0 for which the p-value of the trajectory γ − at u − , say p − * , is the same as the p-value of the trajectory γ + at u + , say p + * , see Figure 6b. These c-values determine the actual speed of the shock-fronted travelling wave solution.
As the stable and unstable manifolds of (0, 0) and (1, −c) are algebraically too complicated to study analytically, we use numerical tools to detect the speeds leading to a feasible desingularised system (26). In particular, we use the function ode45 in MATLAB to obtain the phase plane of (26) and then calculate ∆p := p + * − p − * for different speeds c. Note that we locate the initial points of trajectories approaching (0, 0) or (1, −c) with a small step along their eigenvectors. We find the crossing point of the trajectory leaving from (1, −c) and the straight line u = u − as (u − , p − * ) and the crossing point of the trajectory arriving at (0, 0) and the straight line u = u + as (u + , p + * ). Finally, we calculate ∆p as function of c.
Due to the complexity of numerically simulating a singularly perturbed fourth-order PDE like (10), we simulate solutions of the perturbed ODE system (15) with Matlab's ODE solver 'bvp4c' and compare it with our analytical results from the singular limit. With the diffusivity function and reaction term as above, the numerical results and analytical results coincide (to leading order), see Figure 7. Furthermore, in Figure 7c we compare the numerical and analytical speeds for reaction terms of the form R(u) = 5u(1 − u)(u − A) with varying A. Again, the numerical and analytical speeds coincide (to leading order).

Persistence analysis
For c = c 0 , the orbit in the layer problem connecting u − to u + and the orbits in the reduced problem and desingularised problem connecting 1 to u − and connecting u + to 0 form a complete heteroclinic orbit connecting 1 to 0 in the singular limit ε → 0. Below we will argue that such solution persists in the four-dimensional system (15) for sufficiently small ε, i.e. 0 < ε 1. Note that we do not present the full proof for the persistence claim -which follows from geometric singular perturbation theory (GSPT) based on Fenichel's persistence theorems [11,15,18] since M ± 0 are normally hyperbolic -because this is rather standard, but quite technical, at this stage. Instead, we provide some heuristic arguments for the persistence.
The endpoints of the heteroclinic orbit in the full system (15) are P 0 ε and P 1 ε (17) and the heteroclinic orbit lies in the intersection of the two-dimensional stable manifold of P 0 ε , W s (P 0 ε ), and the two-dimensional unstable manifold of P 1 ε , W u (P 1 ε ), see (19). This intersection will generically not be transversal since the full system is four-dimensional, i.e. 2 + 2 − 1 < 4. Therefore, we extend the full system (15) to a five-dimensional system by appending it with an equation for the unknown speed {c = 0}. That is, we threat c as a variable and not as an unknown parameter. In the extended system the heteroclinic orbit now lies in the intersection of the three-dimensional centre stable manifold W cs (P 0 ε ) and the three-dimensional centre unstable manifold W u (P 1 ε ) and this intersection will generically be transversal since the full system is five-dimensional, i.e. 3 + 3 − 1 = 5. Typically, transversality follows from a Melnikov-type analysis [15,29,36]. We decided to omit this calculation, but its proof is numerically verified in Figure 6(b) and (d). As a result, and for sufficiently small ε, the heteroclinic orbit will persist with a nearby speed c(ε), with c(0) = c 0 , the speed found in the singular limit. Finally, recall that such a heteroclinic orbit corresponds to a shock-fronted travelling wave solution of (10).

Viscous relaxation
In this section, we study shock-fronted travelling wave solutions in (11) and we use similar mathematical techniques as in §2 to obtain a three-dimensional singular perturbed slow-fast system. The reduced problem is the same as in §2, however, it has different algebraic constraints. In contrast, the layer problem is different and only one-dimensional which leads to shocks with different characteristics. Since the methodology of the analysis is similar, we only present a succinct and brief derivation of the main results.

Preliminary observations
The travelling wave solution of (11) of interest here is a solution of that asymptotes to one as z → −∞ and to zero as z → ∞. Here, z := x − ct is again the travelling wave coordinate. Next, with some abuse of notation, we define and transform (28) into a three-dimensional singular perturbed slow-fast dynamical system where u ∈ R is fast variable and (p, v) ∈ R 2 are slow variables. By using a stretched variable ξ = z/ε, (30) is transformed into an equivalent fast system, provided ε = 0, The fixed points of the two equivalent systems (30) and (31) are and we are interested in heteroclinic orbits connecting Q 0 ε with Q 1 ε . The Jacobian of (31) has three eigenvalues with the expansion of ε At Q 0 ε , R (0) < 0, D(0) > 0, thus, τ + 1 (0) > 0, τ − 1 (0) < 0 and τ 2 (0) > 0. Similarly, at Q 1 ε , R (1) < 0, D(1) > 0, thus, τ + 1 (1) > 0, τ − 1 (1) < 0 and τ 2 (1) > 0. That is, the stable manifolds of Q 0,1 ε are one-dimensional and the unstable manifolds of Q 0,1 ε are two-dimensional. At Q A ε , for positive speeds, the stable manifold is two-dimensional and the unstable manifold is one-dimensional; for negative speeds, the stable manifold is one-dimensional and the unstable manifold is two-dimensional.

Layer problem
Letting ε → 0 in (31) gives the layer problem and dp/dξ = 0 and dv/dξ = 0. Thus, we have a two-dimensional critical manifoldM Upon recalling that F (u) = D(u), we observe that the critical manifold loses normal hyperbolicity along the one-dimensional set repelling manifolds for c > 0 and attracting manifolds for c < 0. Similarlŷ is an attracting manifold for c > 0 and an unstable manifold for c < 0, see Figure 8. Considering the stability of the different branches of critical manifold, there may exist connections betweenM ± 0 andM 0 0 and betweenM ± 0 and F ∓ . In contrast to the previous section, we are now interested in connections betweenM ± 0 andF ∓ since we are looking for travelling wave solutions that connect u = 0 and u = 1, and both of these points are onM ± 0 . There are two ways to establish these connections. If c > 0,M ± 0 are repelling and F (u) = −v has two non-repeating real roots β and u l (< α), or u = α and u r (> β). In this case, the related shocks are u r → α and u l → β, see Figure 8a. If c < 0,M ± 0 are attracting. and the related shocks are in the opposite direction α → u r and β → u l , see Figure 8b.

Reduced problem
The reduced problem of (30), obtained by letting ε → 0, is the same as the reduced problem (25) of the previous section and is given by Similarly, its desingularised system 6 is the same and given by However, note that the slow variable p is defined differently, see (14) and (29), and thus has a different meaning.
3.4. The construction of the heteroclinic orbit in the singular limit From the analysis of the layer problem (32), the shocks u r → α and u l → β have positive speeds, while the shocks in the opposite directions, α → u r and β → u l , have negative speed. The shocks u r → α and β → u l potentially relate, in the singular limit, to trajectories of (30) leaving from u = 1 and arriving at u = 0, that is, they have the asymptotic conditions lim z→−∞ u = 1 and lim z→∞ u = 0 we are interested in. In contrast, the shocks u l → β and α → u r correspond to trajectories with the opposite asymptotic conditions lim z→−∞ u = 0 and lim z→∞ u = 1. Thus, we are interested in positive speeds c for which there exist trajectories of the desingularised system (33) that connect (1, −c) with (u r , p * ) and (α, p * ) with (0, 0) (both in forward ψ). Similarly, we are interested in negative speeds c for which there exist trajectories of the The change of speed as function of A where the line with positive speed represents shocks u r → α and the line with negative speed represents shocks β → u l . We remark that we could not find an Allee type R(u) for which both types of travelling wave solutions exist simultaneously.
desingularised system (33) that connect (1, −c) with (β, p * ) and (u l , p * ) with (0, 0). Following the same procedure as in the previous section using ode45 in MATLAB, we can now construct orbits of the correspond to heteroclinic orbits in the singular limit of (30), and thus to shock-fronted travelling wave solutions of (11). See Figure 9 for two prototypical examples of these orbits. One corresponding to a shock-fronted travelling wave solution with positive speed and one with negative speed.

Persistence analysis
To show that these singular orbits indeed persist for ε sufficient small, and thus correspond to shock-fronted travelling wave solutions of (11), we have to proceed in a similar fashion as in the previous section and extend the full three-dimensional system (30) with an equation for the speed {c = 0} (since the stable and unstable manifold W s,u (Q 0,1 ε ) are respectively one and two-dimensional and 1 + 2 − 1 < 3) such that transversality is generically possible. Transversality again follows from a Melnikov-type argument, but we have to extend Fenichel theory near the regular fold pointF ± , where the critical manifold loses normal hyperbolicity -one of the necessary conditions for Fenichel's persistence theorems. This way, we can show that the orbits persist even though in the singular limit we leave, or arrive at, the critical manifold at a fold pointF ± . We decided to not go into the details of this analysis and refer to [5], and references therein, instead, for an outline how the persistence of these singular orbits can be shown. In the end, this shows the persistence of the heteroclinic orbit for sufficiently small ε and with nearby speed c(ε), with c(0) = c 0 , the speed found in the singular limit.

Summary, discussion and outlook
In this article, we studied shock-fronted travelling wave solutions supported by the RDE (1) with a convex nonlinear diffusivity function D(U ) (2) that is negative for U ∈ (α, β) (3), and with an Allee-type reaction-term R(U ) (6). This RDE with forward-backward diffusion was previously derived by [16] from a lattice-based stochastic model modelling a population of individuals and groups that can undergo movement, birth and death events to describe the its macroscopic behaviour. We studied the RDE by adding two different small regularisations; a non-local regularisation −ε 2 ∂ 4 U/∂x 4 , with ε small, see (10) and §2, and a viscous relaxation ε∂ 3 U/(∂x 2 ∂t), see (11) and §3. Note that in the singular limit ε → 0 both PDEs reduce to (1).
These two regularisations allowed us to use a dynamical systems approach to study the shock-fronted travelling wave solutions. In particular, for the nonlocal regularisation the PDE (10) could be reduced to a singularly perturbed four-dimensional system of ODEs (15). As the regularisation term is assumed to be small there is a scale separation in this system of ODEs. This allowed for a further reduction by investigating (15) singular limit in the fast and slow scaling. The singular limit in the fast scaling, called the layer problem (20), described the dynamics near the shock of a shock-fronted travelling wave solutions, and was a two-dimensional Hamiltonian system independent of the speed c, see Figure 5. The singular limit in the slow scaling, called the reduced problem (24), was a singular two-dimensional system of ODEs. It is constraint to the critical manifold M 0 (21) and described the dynamics away from the shock. Note that we use MATLAB to investigate the reduced problem as it is algebraically too involved to determine the sought after trajectories. A shock-fronted travelling wave solution can now be constructed, in the singular limit, upon concatenating the three parts of the solution, see Figures 4 and 6. Subsequently, GSPT can be used to show that the solution persists for sufficiently small ε. Note that the details of this final calculation were omitted, instead it was shown that the dynamics of the full ODE (15) agrees with the obtained results in the singular limit, see Figure 7.
For the viscous relaxation the PDE (10) could be reduced to a singularly perturbed three-dimensional system of ODEs (30). Whilst this ODE had the same reduced problem as with the non-local regularisation, it had a different layer problem (32). This difference can lead to shock-fronted travelling wave solutions with different characteristics for same nonlinear diffusivity function D(U ) (2) and reaction-term R(U ) (6), see Figure 10. In addition, as the shock-connection in the layer problem is at a point where the critical manifold loses normal hyperbolicity, GSPT has to be extended to prove the persistence of the singular orbit for sufficiently small ε. Again, details of this computation were omitted.

Regularisations and the lattice-based stochastic model
While the two regularised PDEs have the same singular limit (1), the different regularisations yielded shock-fronted travelling wave solutions with different characteristics. Therefore, we mainly compared the singular limit results of the two models with the travelling wave ODE systems, and not with the numerical results of (1). The reason for this is that the numerical schemes used to simulate (1) naturally introduce artificial regularisation (and error) terms and, as shown in this article, different regularisations yield shock-fronted travelling wave solutions with different characteristics. The connection between the numerical results of (1) and the analytical results therefore needs to be further explored.
In addition, (1) was derived from a lattice-based stochastic model and during this derivation of the continuous description small higher order terms were omitted. Including some of these small higher order terms would potentially result in a (differently) regularised version of (1), which in turn could lead to shock-fronted travelling wave solutions with different properties. Therefore, studying the connection between the lattice-based stochastic models and the regularisations is also an interesting topic.
For instance, a natural question to ask is what happens when we consider a linear combination of the non-local regularisation (considered in §2) and viscous regularisation (considered in §3)  where µ ∈ [0, 1] is a constant. Note that µ = 0 corresponds to the viscous regularisation (11) and µ = 1 corresponds to the non-local regularisation (10). The associated four-dimensional slow-fast system 7 is given by The corresponding layer problem, for µ = 0, is If v is such that (34) has three fixed points (u − , w − ), (u 0 , w 0 ) and (u + , w + ), where u + < α < u 0 < β < u − . Then, for µ = 1, (34) does not have heteroclinic orbits connecting (u − , w − ) with (u + , w + ). Hence, we do not expect shock-fronted travelling wave solutions in this case.

Generalisations
In this article, we concentrated on a specific quadratic nonlinear diffusivity function D(U ) (2) centred around 2/3 and a specific Allee-type reaction-term R(U ) (6) as these were derived from an underlying lattice-based stochastic model [16]. However, the techniques used in this article can in fact be easily extended to more general nonlinear diffusivity functions and reaction terms. For instance, if we change the reaction term from an Allee type (6) to a logistic type (5) (as studied in [21]), we can still construct the higher-dimensional systems based on the two regularisations (10) and (11). Since the two layer problems (20) and (32) only depend on F (u), the anti-derivative of D(u), and not on R(u), we obtain the same conditions for the shocks as for the Allee type reaction term. That is, for the non-local regularisation the shocks will have, to leading order, endpoints u − and u + , while the shocks will have, to leading order, endpoints u r and α or u l and β for the viscous relaxation. In other words, the size of the shock depends on the relaxation and the nonlinear diffusivity function D(U ), but not the reaction term R(U ). For both regularisations, the reduced desingularised problem has four fixed points which are determined by the roots of the product of the nonlinear diffusivity function D(U ) and the reaction term R(U ). In particular, the fixed points are (0, 0), (1, −c), (α, −cα) and (β, −cβ). In the desingularised system, the fixed point (0, 0) is a stable node or stable spiral for c > 0 and an unstable node or unstable spiral for c < 0. For shock-fronted travelling wave solutions with the asymptotic conditions lim z→−∞ U = 1 and lim z→∞ U = 0, we expect (0, 0) to be stable in the desingularised problem. Therefore, we expect those travelling wave solutions to have positive speeds. Hence, if the reaction term is logistic, we do not expect shock-fronted travelling wave solutions with negative speeds. However, using other boundary conditions may provide novel characteristics, see [9,10] for examples of moving boundary problems with logistic type reaction terms.

Stability
Another natural extension of this work is to analyse the stability of the constructed shock-fronted travelling wave solutions. This was partly done for smooth travelling wave solution supported by (1) with D(U ) as in (2) and logistic reaction term R(U ) (5) in [21]. In that article we studied the absolute spectrum of the associated desingularised stability problem and showed that for speeds above the minimal wave speed, the essential spectrum [19,30] of the desingularised system can always be weighted into the left-half plane, while this is not possible for speeds below the minimal wave speed [21]. This analysis can be repeated for the shock-fronted travelling wave solutions constructed in this article since the essential spectrum is related to the behaviour of the wave at infinity and thus only determined by the asymptotic end states of the shock-fronted travelling wave solution under consideration. For brevity we decided not to show this computation and instead refer to [21]. In short, the computation shows that the essential spectrum of the associated desingularised stability problems of (1), (10) and (11) are all fully contained in the left-half plane, see Figure 11, thus there are no absolute instabilities. However, what remains to be determined is the point spectrum, as well as the connection of the essential spectrum of the associated desingularised stability problem and (Λ) Figure 11: The essential spectrums (shaded green regions plus boundaries) of the desingularised stability problems associated to (1) (a), (10) (b) and (11) (c) with ε = 0.1. Note that the essential spectrums are to leading order the same and fully contained in the left-half plane.
the original stability problem, to complete the linear stability analysis. This is part of future work, see also the discussion in [21].
Similarly, integrating (A.1) between u b (> u 2 ) and 1 gives which, for c > 0, leads to the necessary condition

Appendix B. The heteroclinic orbits of the layer problem
We derive the analytic expressions for the heteroclinic orbits given in the layer problem supported by where v is a constant. Based on its Hamiltonian, we require H(u, w) = − 1 2 w 2 + G(u) + vu = 0, on the heteroclinic orbits u 0,± h , w 0,± h . Subsequently, we obtain w = ± 2 (G(u) + vu).
Note that G(u) has two integration constants. With specific integration constants, w(u) can become a second-order polynomial with specific roots. That is, we can write w as w(u) = ± 2 (G(u) + vu) = ± a 2 (u − B 1 ) 2 (u − B 2 ) 2 .