Slow travelling wave solutions of the nonlocal Fisher-KPP equation

We study travelling wave solutions, , of the nonlocal Fisher-KPP equation in one spatial dimension, with and , where is the spatial convolution of the population density, , with a continuous, symmetric, strictly positive kernel, , which is decreasing for x  >  0 and has a finite derivative as , normalized so that . In addition, we restrict our attention to kernels for which the spatially-uniform steady state u  =  1 is stable, so that travelling wave solutions have as and as for c  >  0. We use the formal method of matched asymptotic expansions and numerical methods to solve the travelling wave equation for various kernels, , when . The most interesting feature of the leading order solution behind the wavefront is a sequence of tall, narrow spikes with weight, separated by regions where U is exponentially small. The regularity of at x  =  0 is a key factor in determining the number and spacing of the spikes, and the spatial extent of the region where spikes exist.


Introduction
In this paper we study travelling wave solutions of the nonlocal Fisher-KPP equation in one spatial dimension, where D is a positive constant and the kernel, φ(x), is continuous, symmetric, strictly positive and decreasing for x > 0, with finite derivative as x → 0 + , normalized so that We are interested in permanent form travelling wave (TW) solutions of (1). Numerical simulations suggest that minimum speed travelling waves are formed in the solution of initial value problems with sufficiently localized initial conditions. We focus on the case D 1, which we will see is the limit in which interesting, highly nonlinear behaviour can be found.
TW solutions with speed c are of the form, u(x, t) ≡ U(z), where z = x − ct, and satisfy the nonlocal ordinary differential equation We will confine our attention to TW solutions that satisfy with c > 0 by studying kernels for which it is straightforward to show that the uniform solutions u = 0 and u = 1 of (1) are respectively unstable and stable. TW solutions with these boundary conditions are right-propagating and connect the uniform steady states, U = 1 and U = 0. Note that discontinuous kernels, [1], and kernels that become negative, [2], can lead to instability of the steady state u = 1, as can some smooth positive kernels such as φ ∝ e −x 4 and φ ∝ (1 + x 4 ) −1 . Instability of u = 1 leads to the existence of periodic travelling wave solutions, which are not the subject of the present work (see, for example, [3]). The nonlocal Fisher-KPP equation, (1), is relevant to many different scientific areas, [3][4][5], and has been studied extensively, but mainly through proofs of the existence of travelling wave solutions, [6][7][8], and in the limit of fast travelling waves ( D 1). Fast travelling wave solutions are a small perturbation of the travelling wave solution of the local Fisher-KPP equation, which is equivalent to (1) with φ(x) = δ(x), [9]. When D 1, (1) can be characterised as an equation with short range activation and long range inhibition, which is known to lead to the generation of localized spikes in other reaction-diffusion systems, [10].
It is helpful to work in terms of the variable L(z) ≡ log U(z), so that (3) and (4) become subject to with since the TW solution must connect with the stable manifold of the steady state U = 0 as z → ∞. This transformation was also used in [11], which considered the formation of spikes, represented by delta functions, in steady and unsteady solutions of (1) on a finite domain using a kernel with compact support. Equation (7) suggests, as was proved rigorously in [6], that a TW solution with U > 0 only exists for c 2 √ D. It is therefore convenient to write D = D 0 c 2 with D 0 1 4 and consider the solutions of (5) when c 1. We will find below that the diffusive terms in (5) are only significant ahead of the wavefront when c 1, and we will later assume that D 0 = 0 in order to simplify parts of the asymptotic analysis of the solution behind the wavefront.
We begin in section 2 with a preliminary numerical investigation of the solution of (5) for the two most obvious choices of kernel, This motivates a discussion of how to investigate the very different forms of TW solution in these two cases through choices of kernel that are intermediate between the two, which is presented in section 3. In section 4 we construct the asymptotic solution when φ(x) has a discontinuous derivative at x = 0. We study the kernel φ = Φ 1 in detail, since it allows (5) to be written as a single ordinary differential equation. We then show how to generalize this analysis to other kernels with discontinuous derivative at x = 0 by writing down an ansatz for the behaviour of U for all z, which includes a delta function spike behind the wavefront. We then discuss the asymptotic solution of (5) for other kernels, in particular a family of kernels (discussed in section 3) that allows us to unfold the bifurcation from a single spike solution when φ = Φ 1 to a solution with infinitely-many spikes that exist in a finite region z ∞ < z < 0 when φ = Ψ 2 ≡ 1 4 (1 + |x|) e −|x| , a kernel that is three-, but not four-, times differentiable at x = 0. In section 5 we show that kernels that are twice-, but not infinitely-, differentiable have a similar asymptotic structure to the solution with φ = Ψ 2 , and discuss how the region z ∞ < z < 0, where the spikes exist, grows as n increases for a family of kernels φ = Φ n . We conclude in section 6 that infinitely-many spikes exist in any semi-infinite region z < z 1 for all z 1 < 0 when φ = Φ ∞ and for other infinitely-differentiable kernels.

Numerical solutions of the Fisher-KPP travelling wave equation for two typical kernels
Natural choices of kernel, φ(x), are either the exponential kernel, φ = Φ 1 , or the Gaussian kernel, φ = Φ ∞ , (8). Although these kernels are most obviously distinguished by the more rapid rate of decay of the latter as x → ±∞, we find that it is only their behaviour as x → 0, where the former has a discontinuous derivative and the latter is infinitely-differentiable, that has a significant effect on the qualitative nature of TW solutions of (5) as c → 0. In order to gain some insight into TW solutions of (5), we begin by examining numerical solutions for these two kernels. After truncating to a finite domain, we use central-differences to evaluate the spatial derivatives in (5) and two-point Gaussian quadrature to compute the integral, taking a linear variation in L(z) across each element, collocating at the midpoint. We took the translational invariance of (5) into account by fixing the area under the solution in the truncated domain of solution. We solved the resulting system of nonlinear algebraic equations using fsolve in MATLAB, which uses the trust region dogleg method, providing the analytical Jacobian to this routine to increase its execution speed. We proceeded by continuation from c = 0.1 down to a small value of c, using adaptive regridding to resolve the steep gradients of the solution, with the length scale estimated from U/ |U zz | at each local maximum.
The TW solution with kernel Φ 1 (x), wavespeed c = 10 −12 and diffusion coefficient D 0 = 1 4 is shown in figure 1. The top panel shows the form of the TW solution on an O (1) lengthscale. There appears to be a single, large, narrow spike at the wavefront. Note that the U-axis has been restricted; the actual maximum value of the spike is close to 3 × 10 5 as can be seen in the bottom right hand panel, where the z-axis has been restricted so that the full structure of the spike can be seen. Note the asymmetry of the spike, with U decaying more rapidly as z increases than it does as z decreases. The spike contains an approximately unit area, which leads us to believe (correctly, as we shall see) that the leading order solution can be written as U(z) ∼ H(−z) + δ(z) as c → 0, up to a translation in z, with H the Heaviside step function. There is, however, some small scale structure to the solution, as can be seen in the bottom left hand panel of figure 1. There is an additional, smaller asymmetric spike to the left of the large spike, along with some further oscillations as z decreases. Figure 2 shows that the heights of the spikes scale with c − 1 2 and c − 1 8 , whilst their separation scales with c 1 4 . We will analyse this structure in detail in section 4.1.
We now turn our attention to TW solutions with the kernel Φ ∞ . The solution for c = 10 −10 (note that the existence of so many large spikes means that we are unable to resolve the solution adequately with a reasonable number of grid points for smaller values of c) and D 0 = 1 4 is shown in figure 3, where it is immediately clear that the structure of the solution is very different to that shown in figure 1. There is a sequence of spikes behind the wavefront with O(1) spacings. The bottom panels of figure 3 suggest that each is locally Gaussian, in contrast to the asymmetric spikes shown in figure 1 when φ = Φ 1 . Figure 4 indicates that the height of each spike scales with c −1/2 as c → 0, but decreases with each successive spike (numbered starting from the largest spike, which lies at the wavefront, and proceeding in the direction of decreasing z). The final panels of figure 4 show how the weight, w m , (computed as the integral of U(z) between successive minima) and spacing, ∆z m , between the spikes behaves with spike number, m. It appears that w m ∼ ∆z m as m → ∞, but it is not clear whether ∆z m approaches a finite value or tends to zero very slowly as m → ∞. In either case, this suggests that the spikes fill the region behind the wavefront as c → 0. We will discuss this point further for Φ ∞ and other infinitely-differentiable kernels in section 6.
These numerical solutions of (5) suggest that the structure of the TW solution is very different for these two typical kernels. In the rest of this paper, we will try to understand what controls these differences. We begin by presenting some families of kernels that allow us to unfold the bifurcation structure of the TW solution as the kernel changes from Φ 1 to Φ ∞ .

Some families of kernels
A fact that we will use extensively in section 4 is that Φ 1 (x) ≡ 1 2 e −|x| is a Green's function, since it satisfies The smooth kernel is not a Green's function. However, consider the bounded function Φ n (x) that satisfies For positive integer values of n, this is easy to solve analytically, for example, In general, Φ n has 2n − 2 continuous derivatives at x = 0. The Fourier transform of (10) is and hence This shows that it is consistent to , and that the sequence of kernels defined by (10) connects Φ 1 to Φ ∞ as n increases. Moreover, it is also possible to invert (11) for non-integer values of n to obtain the general result, where K n− 1 2 is a modified Bessel function of the second kind. This gives us a means of varying a parameter, n, that transforms the kernel Φ 1 to Φ ∞ as n → ∞ whilst smoothly changing the regularity of the function at the origin, since the Taylor series of Φ n has a non-analytic part of O(|x| 2n−1 ) as |x| → 0. The kernel Φ n (x) is shown in figure 5 for various values of n.
We will see below that it is also of interest to consider kernels whose slope at the origin varies smoothly without a change of regularity. An obvious candidate is φ(x) = (1 − k) Φ 1 (x) + kΦ 2 (x) with 0 k 1. However, by defining Ψ n (x) to be the bounded solution of and usingΦ a kernel that reduces (5) to a very simple ordinary differential equation, we are able to investigate analytically the effect of changing the slope of φ(x) at the origin. We will also consider the kernelΦ along with a couple of others, to illustrate additional possible features of TW solutions. Note that The results that we obtain in this paper are summarized in table 1.

Kernels that are not twice-differentiable
We begin our investigation of kernels that are not twice-differentiable at the origin by considering φ(x) = Φ 1 (x) in section 4.1. We then generalize this analysis to similar kernels in section 4.2. In sections 4.3 and 4.4 we study the kernel φ(x) =Φ k (x), and show that, as k increases, a second spike appears when k = 3 4 and that the number of spikes tends to infinity as k → 1.    we can use the fact that Φ 1 satisfies (9) to rewrite (5) as and hence, by eliminating v, subject to (6). Being able to reformulate (5) as an ordinary differential equation is extremely helpful, as we can use formal asymptotic methods without having to postulate a global form for U(z) ≡ e L(z) to substitute into the convolution integral. In its current form, (18) gives L = 0 when c = 0, which reproduces the boundary condition (6). Two obvious rescalings that bring in more terms at leading order are L = O(1/c), with L < 0, z = O(1), which we will find gives the leading order problem ahead of the wavefront, and L = O(1), z = O(c 1/3 ), which controls the solution behind the wavefront, and will be the starting point for our analysis. We define At leading order, When L 1, this gives Lzzz ∼ −L, so that with l 0 and z ∞ constants. Note that when L is large and negative, Lzzz ∼ 1, so that L increases like z 3 . When L is large and positive, Lzzz is large and negative and L rapidly decreases. These qualitative observations suggest that the oscillations in the solution that are present

Countably infinite
Semi-infinite, z < 0 as z → −∞ persist as L moves away from equilibrium. A typical solution, obtained using ode15s in MATLAB, is shown in figure 6 and confirms that the initial oscillations grow as z increases. Although the solution cannot have L → −∞ as z → ∞, the rapid growth of −L shown in figure 6 means that we are unable to find a numerical solution beyond the value of z shown. In order to understand how to construct the asymptotic solution of the full problem, we need to understand the structure of solutions such as that shown in figure 6 as L varies between a large negative local minimum, a large positive local maximum and back to a large negative minimum. We can do this by considering the asymptotic solution of (21) subject to the initial conditionsL with Â = O(1) as → 0, so that we have, without loss of generality, placed a large negative minimum in L with L = − −1 at z = 0. Note that this approach was used in [12, appendix A] to study a similar third order ordinary differential equation. As we shall see, the solution in this case can be constructed as a sequence of asymptotic regions centred on either a local minimum (Region A m ) or a local maximum in L (Region B m ).

Region A m+1 .
In region A m+1 , we rescale using subject toL Figure 6. A typical solution of (21). Although we know that the solution cannot have L → −∞ as z → ∞, it changes so rapidly that the numerical method used, ode15s in MATLAB, is unable to compute accurate solutions for L large and negative.
To all algebraic orders, this giveŝ which has a local maximum at . We shall see below that matching to the preceding Region B m+1 forces L to be zero at this maximum, and hence Â ≡ (3/2) 1/3 ≈ 1.145. This means that the remaining root of (26) is simple, and lies at ẑ =Â. The solution is shown in figure 7. As ẑ →Â, L ∼Â 5 (ẑ −Â), and L becomes positive for ẑ >Â. We must therefore rescale into a new asymptotic region, which we denote Region B m .

Region B m .
By matching with Region A m+1 , we find that appropriate scalings arê In terms of these scaled variables, (24) becomes In this region we need to develop a two term asymptotic expansion, so we writē At leading order, we havē subject to the matching condition Note that (29) is invariant under the transformation Although this does not allow us to solve (29) analytically, we can solve (29) numerically subject to the generic boundary condition to obtain the solution L 0 (z) =L(z), which is shown in figure 8. This has maximum value L max ≈ −0.867. For L large and negative, eL is exponentially small, so L ∼ − 1 2 az 2 as z → ∞, and we find that a ≈ 1.302. It is this asymmetry in the behaviour of L in Region B m that leads to the asymmetry of the spike in the TW solution, shown in figure 1. From the boundary condition (30) and the invariance (31), we can see that , and we must rescale into a new asymptotic region, which we label Region A m . At leading order

Region A m . We define new scaled variables
to be solved subject to the matching condition This has solution with a minimum of L min = − 2 3ā at z = 2ā. By matching the solutions in Regions A m and A m+1 , where L is cubic at leading order, through Region B m , we have been able to obtain the leading order asymptotic values of the maximum and minimum in L that follow a large, negative minimum in L . After rewriting the solutions in terms of the original variable, L , we find that after the first minimum, the following maximum is and the following minimum This shows that a deep minimum is followed by a moderate maximum, and then a much deeper minimum, as can be seen in figure 6. By solving (19) subject to (22) for various values of l 0 and z 0 , we can plot the values of the maxima and minima of L in figure 9 and find excellent agreement between the numerical solution and the leading order asymptotic solutions, (34) to (36). Note the huge difference in size between the maxima and minima. For example, It is this disparity in scales that leads to the sequence of spikes of rapidly decreasing height that we discuss below.  (19), which are neglected in the analysis above, become active in Region A 1 . By considering the size of these terms after rescaling, we find that we need = O(c 1/4 ) as c → 0, so without loss of generality we take = c 1 4 . The leading order problem in Region A 1 is then which, using the fact that it comes from (17), can be written as Note that this has L ∼ − 1 2z 2 as z → 0. The factor of ā has been scaled back into the solution in the other asymptotic regions, but we omit the details here (see section 4.2 for the general case). In addition, L ∼ k +z as z → ∞, where the solution remains uniform. The constant k + is defined in (7), and it is here, ahead of the wavefront, that the minimum wave speed condition, c 2 √ D, and hence D 0 1 4 , is determined. We now know the asymptotic scalings for the successive Regions A m+1 , B m and A m , along with the fact that = c 1/4 in Region A 1 . The scalings of z in the Regions A m give the distances between the sequence of spikes, whilst the scalings of z, L and hence U, in the Regions B m give the widths and heights of the spikes. By comparing successive Regions A m , we find that, in terms of z = c −1/3 z, the (large) width decreases by a fourth root, as does the size of L as m increases to m + 1. The width of the spike Region B m scales with the inverse square root of the width of Region A m , whilst the size of U in Region B m scales with the square root of the size of L in Region A m . When we put all of this together in terms of the original variables, z and U = e L , we find that the asymptotic scalings are:  figure 1. However, the scalings for the mth spike are only asymptotic provided that c is small enough that the height of the spike is large, so that the number of spikes in the solution is For c = 10 −12 , the value used in figure 1, m 0 (c) ≈ 2.4, so we see just two spikes. For z < z m0 , the solution has L 1, and hence U ∼ 1, at leading order. To obtain, for example, m 0 ≈ 3, we would need c ≈ 10 −28 so, for all practical purposes, there are no more than two clearly defined spikes. Figure 10(a) shows the numerically-calculated heights of the first two spikes in U and their leading order asymptotic approximations. The agreement is very good, noting that the second spike has height of O(c −1/8 ) and is far slower to converge to the leading order approximation than the first spike. Figure 10(b) shows the numerically-calculated spacing between the first two spikes, which is also in excellent agreement with the leading order asymptotic solution.
The distance between the first two spikes. The broken line shows the asymptotic approximation

Generalization to other non-differentiable kernels
We have been able to construct a complete asymptotic solution when φ = Φ 1 ≡ 1 2 e −|x| because this is the Green's function for a simple, constant coefficient, ordinary differential equation, so that (5) can be written as a single ordinary differential equation. If φ is not twice-differentiable at the origin but not a Green's function, and we can formally write we can show that the asymptotic structure of the solution is qualitatively the same as for φ = Φ 1 (note that for φ = Φ 1 , κ = 1). This means that after differentiating (5) twice we obtain We can now proceed by using the ansatz In other words, we assume that the leading order structure of the solution is the same as that for φ = Φ 1 . If this gives a consistent asymptotic solution, we have a post hoc justification for using (43). In (43), the function U ∞ (z) is the leading order solution for −z c 1/3 . When φ = Φ 1 , U ∞ = 1, but for other kernels, U ∞ need not be constant.
On substituting (43) into (5) and (42), we obtain and (45) It now becomes clear after using the scalings (37) and (38) that, apart from some differences in constant terms, we recover the same leading order problems and structure that define the solution when φ = Φ 1 . Specifically, in Region A 1 , where L = O(c −1 ) and z = O(1), the solution of (44) matches with the solution in Region B 1 provided that L z (0) = 0, so that In the far field, for −z (47) Note that when φ = Φ 1 , (47) has solution U ∞ = 1, Ū = 0. Indeed, (47) shows that the only kernel that leads to the constant solution U ∞ = 1 is φ = Φ 1 and scalings thereof.
To illustrate how this works, we solved (47) numerically, using the trapezium rule to evaluate the integrals, for the kernel The asymptotic solution, along with the numerical solution for c = 10 −12 , is shown in figure 11, and the agreement is excellent. The numerical solution of the full problem also shows the spike at z = 0 which appears as a delta function when z = O(1) as c → 0. We use the kernel (48) to illustrate, firstly that this asymptotic solution structure correctly describes the solution for a kernel that does not allow (5) to be reduced to an ordinary differential equation, and secondly that, since this kernel decays like a Gaussian as |x| → ∞, the far field behaviour of the kernel is not its most important feature; the behaviour in the neighbourhood of the origin determines the structure of the solution as c → 0. In Regions B m , where the spikes exist, (45) shows that the leading order equation is After scaling z to remove the factor of κ, this is the leading order equation that we found for φ = Φ 1 . In Regions A m with m > 1, (45) gives at leading order as c → 0. Note that for φ = Φ 1 , A = 1 and that for the kernel given by (48), A ≈ 0.255. The solution is again cubic in z, and the whole asymptotic structure is the same. By matching together these asymptotic regions we are able to obtain the leading order asymptotic solution and in particular show that, if we define Region B m lies at z = z m , where as c → 0. Figure 10 shows the excellent agreement between these leading order asymptotic expressions for φ = Φ 1 . Figures 12 shows the same comparisons for the kernel given by (48) which are also in good agreement, although the rate of convergence is somewhat slower for the second spike.

Numerical solution for φ =Φ k
Although it appears at first sight that we now have a complete asymptotic theory, there are other possibilities, which we can investigate by considering the kernel Φ k (x) defined by (14). The slope of this kernel changes by (1 − k) at x = 0, and Ψ k = Ψ 2 when k = 1. Figure 13 shows the travelling wave solution for c = 10 −12 and various values of k. For k < k 1 (we will demonstrate below that k 1 = 3 4 ), the structure of the TW solution is of the form (43), with a single spike of O(1) weight separating regions where U > 0 and U = 0. However, as k → k 1 , U ∞ (0) → 0 and as k increases past k 1 a new spike of O(1) weight is formed. The original spike remains, isolated at z = z 1 , whilst the new spike, at z = z 2 , has the sequence of small spikes of decreasing height to its left at z = z m with m > 2.
These numerical solutions suggest that there exists a strictly increasing sequence, k M , with k 0 = 0 and k M → 1 as M → ∞, and, for k M−1 < k < k M , a strictly decreasing sequence z m (k) for m = 1, 2, . . . M, with z M → z ∞ > −∞ as M → ∞, such that the leading order asymptotic solution for c 1 is In addition, the sequence of small spikes, with small weight as c → 0, lies at z = z m with m > M, to the left of the final spike of O(1) weight at z = z M . As this notation suggests, it seems natural to think of the solution as consisting of an infinite sequence of spikes at z = z m , with the functional form of the kernel controlling which of these have weight of O(1) as c → 0. The solid lines in figure 14 show the weights, w m (k), and spacings, ∆z m (k) = z m − z m+1 , for the travelling wave solution with c = 10 −12 as a function of k, which illustrates the emergence of seven distinct spikes as k → 1. For smaller values of c, which we were unable to access numerically, we would expect further members of the infinite sequence of spikes to be visible in the numerical solution.

Asymptotic solution for φ =Φ k as c → 0
We can provide evidence for our assertion that an infinite sequence of spikes of O(1) weight emerges as k → 1 when φ =Φ k by constructing the asymptotic solution as c → 0. Note that in this case the TW equation (5) can be written as At leading order as c → 0 with z = O (1) and U, L = O(1), we find that the solution that is bounded as z → −∞ is On substituting the ansatz (52) into (3), for z z M , we obtain, at leading order as c → 0, since U > 0, For this simple kernel, we can evaluate the integral, which gives the product of a linear factor and an exponential, and hence find two relationships between A and w m , including which we will need later (see appendix for more details). For 0 < k < k 1 , since M = 1, evaluating (55) is sufficient to find the two unknowns, A and w 1 , (translational invariance allows us to take z 1 ≡ 0) as  (57).
Since (54) shows that U ∞ (z M ; k) = 1 + A, we require that A −1 for this single spike solution structure to be available. Using (57), this means that we need k 3 4 , and hence that k 1 = 3 4 , consistent with the numerical solutions shown in figure 13. Note that the weight of the spike, w 1 (k), increases monotonically from one to two as k increases from zero to k 1 . These asymptotic expressions for w 1 and A for 0 k k 1 are in excellent agreement with numerical solutions of (5) (see figures 14 and 15).
For k > k 1 we expect multiple spikes, as given by (52), and we need to consider the asymptotic solution in the regions between the spikes, where U 1 when c 1. Specifically, we expect that L = O(c −1 ) with L < 0, so that U is exponentially small in the regions between the spikes. If we define L = cL, with L = O(1) as c → 0, (note that this is a new definition, and not the notation used in section 4.1) we find that, at leading order, (5) gives with F(z; k), since the ansatz (52) remains the same, given by (55), which determines L z in terms of F. In particular, matching to the spike regions requires L z (z m ) = 0, and hence It is helpful to treat these equations as a system that determines w m for given values of z m . Matching also requires that L = 0 at z = z m , and hence that This can be treated as a system of equations that determines z m for given values of w m . For D 0 > 0, L z is related to F(z; k) nonlinearly through the solution of the quadratic equation (58). However, we have found that, even for D 0 = 1 4 , the maximum value of D 0 , which corresponds to the minimum speed TW solution, D 0Lz is numerically very small (typically less than 0.01) in all the cases that we have studied in this paper. In other words, the solution with D 0 = 0 is adequate to describe, with great accuracy, the solution for 0 < D 0 1 4 . As we saw in section 4.1, the term D 0L 2 z is crucial in the region z > 0, ahead of the wavefront, and determines the wavespeed, but between the spikes, this diffusive term is never relevant. We therefore focus our attention on the case D 0 = 0, since (58) which can be evaluated analytically.
In order to determine the positions and weights of the spikes, we need to solve (56), (59) and (61). Note that (56) and (59) are linear in A and w m , so that these variables can be eliminated. For a solution with M spikes, (61) then provides M − 1 nonlinear equations for the M − 1 distances between the spikes, ∆z m . We solve these in MATLAB using continuation in k, increasing the number of spikes in the solution by one as A + 1 passes through zero. Note that the linear system (59) is extremely poorly conditioned even for moderately large values of M, with the condition number of the corresponding matrix increasing exponentially with M. This is because the coefficients in (59) become very close to each other when z m − z M is small. We were able to compute solutions with up to 18 spikes by using variable precision arithmetic in MATLAB. The asymptotic solution that we calculate is shown as broken lines in figure 14. This is in excellent agreement with the numerical solution, shown as solid lines. We are able to calculate many more spikes asymptotically than we can using the full numerical solution with c = 10 −12 . The function A(k) is shown in figure 15, along with the analytical expression (57).
Before we move on to consider the TW solution for differentiable kernels, note that the solution for φ =Φ k described above, with a new spike being born next to the leftmost spike as k increases, and a monotonically-decreasing solution U ∞ (z) behind this spike, is actually a special case. We can illustrate this by briefly considering the solution with φ =Φ k , given by (15). In this case, (5) can be written as the equivalent system of ordinary differential equations Proceeding as we did for φ =Φ k , as c → 0, the leading order equation In contrast to (54), the solution of (63) is not monotonic, so there is no guarantee that if it reaches zero (indicating the birth of a new spike) this will be at the position z = z M of the leftmost spike. It is possible to determine the bounded solution of (63) analytically when there is just one spike. We do not give the details here, but note that we find k 1 ≈ 0.907 94, with the point where U ∞ reaches zero occuring away from the leftmost spike. This is consistent with the numerical solution of (5) for c = 10 −12 shown in figure 16, where it is clear that U ∞ is not monotonic. Subsequent spikes are also born via this mechanism.

Kernels that are twice-, but not infinitely-, differentiable
We begin by discussing the governing equation in the spike regions when φ(x) is twice-, but not infinitely-, differentiable. If the first discontinuous derivative of φ is the p th with p > 2 (we will discuss the case of infinitely-differentiable kernels in section 6), and if, following (40) and (41), we can write then, after differentiating (5) p times we obtain, neglecting the diffusive terms as c → 0, In the narrow spike regions, the only possible leading order balance is between the two terms on the left hand side of (66). However, for a spike with weight of O(1), so that e L = O(dL/dz), the leading order equation for p > 2, is simply d p+1 L/dz p+1 = 0, and the first p derivatives of the outer solution match across the spike region. In particular, apart from a constant correction of O(log c), the leading order solution simply reproduces the quadratic term in the inner limit of the outer expansion for L ≡ cL, and the spike is therefore Gaussian with width of O(c 1/2 ), weight of O (1) and hence height of O(c −1/2 ), all determined from matching to the outer solution, consistent with numerical solutions of the full TW problem with c 1. The spike regions are entirely passive in terms of the asymptotic structure of the solution.
The results of section 4 suggest that when the kernel has a continuous derivative at x = 0, an infinite number of spikes of O(1) weight exist behind the wavefront, so we use an ansatz of the form We will confine our attention to the sequence of kernels φ = Φ n , for which we can reduce the TW equation to an ordinary differential equation and U ∞ = 1, and only consider the case of n a positive integer, even though the kernel (12) can be used in (5) for non-integer n. It would be interesting to use fractional calculus as a framework to investigate the non-integer case in detail (see (83) below), but we have not attempted this here.
Guided by the analysis of the case φ =Φ k given in section 4.4, we postulate that the asymptotic solution of (5) for c 1 when φ = Φ n takes the form (67) with U ∞ = 1 (the unique solution when L = O(1) as c → 0). Proceeding as before, we find that the positions and weights of the spikes are determined by the infinite system of equations where In order to solve (68) and (69) numerically, we need to know the functional form of the far field solution, w m and ∆z m ≡ z m − z m+1 for m 1. We can gain some insight into this by considering the case φ = Ψ 2 . Note that Ψ 2 =Φ k when k = 1. Figure 17(a) shows w m and ∆z m for the numerical solution of (59) and (61) when 1 − k = 10 −10 , which we expect to be close to the solution when k = 1, but with a finite number of spikes. This strongly suggests that both w m and ∆z m decay exponentially as m → ∞, and that w m , ∆z m = O(δ m ) for some 0 < δ < 1 is an excellent approximation even for moderate values of m. We also note from figure 17(b) that w m /∆z m appears to asymptote to a constant ratio as m → ∞. We can confirm this analytically for φ = Ψ 2 ≡ 1 4 (1 + |z|) e −|z| and D 0 = 0, in which case If we define F m (z) to be F(z) for z m+1 z z m, then, noting that eL /c 1 for L = O(1) < 0, we find by solving the leading order equation By applying the three conditions F m (z m+1 ) = F m (z m ) = 0 and zm zm+1 F m (z)dz = 0, we can determine the constants a m , b m , c m and d m in terms of the as yet unknown slope F m (z m+1 ), using computer algebra in Mathematica. Now noting that, from the definition (70), F(z) is as smooth as φ(z − z m ) at z = z m , we can see that F and F are continuous at z = z m . This leads to the first equation being simply a continuity condition and the second coming from (72) with the constants determined by computer algebra, which gives After substituting the definition of φ = Ψ 2 into (70), evaluating the integral and comparing to (72) we find that By comparing these expressions in the limit m → ∞ with those that we obtained from Mathematica, we find that there are two constraints that the far field behaviour of F m (z m+1 ) must satisfy, namely and Along with (73) and (74), these expressions strongly suggest that w m = O(∆z m ) as m → ∞. If we write and substitute into (73), (75) and (76), we find, after rearranging and summing various geometrical series, that (78) On simplifying, we find that δ is a solution of The only root with 0 < δ < 1 is The lower panel shows the value of δ calculated in the same way, either from curve fitting to w m or ∆z m . As indicated on the figure, either method of estimating δ gives an almost identical result. In each case, the squares are the values obtained from the far field asymptotic analysis described in section 5.
and hence The broken lines in figure 17 show that these constants are in excellent agreement with the numerically-calculated far field solution for m 5. This far field analysis can also be carried out for the kernels φ = Ψ n , making use of the continuity of higher derivatives of F n . This has to be done entirely using computer algebra, and, since the polynomial that determines δ is of order greater than four for n > 2, no analytical expressions for the solutions are available. We find that the kernel φ = Φ n gives n − 1 possible values for both δ and W ∞ /∆Z ∞ . The values corresponding to the largest value of δ in each case are plotted for n 7 in figure 18. These are shown in comparison to the values obtained from the numerical solution of (5) with c = 10 −12 and φ = Φ n by curve fitting. These are inevitably rather crude estimates because of the limited number of spikes that emerge in each case, even for such a small value of c. However, the agreement is good, and indicates that the far field asymptotic analysis captures the behaviour of the solution. Note that figure 18 includes some results for 1 < n < 2 that are consistent with δ → 0 and z ∞ → 0 as n → 1. For values of n sufficiently close to one it is difficult to resolve the tightly packed spikes that exist behind the wavefront in the full numerical solution. As n → ∞, a key observation, see figure 19, is that the far field and numerical solutions strongly suggest that δ → 1 and W ∞ /∆Z ∞ → 1, and . This is therefore consistent with our earlier assertion that for the Gaussian kernel, φ = Φ ∞ , the TW solution consists of an infinite sequence of spikes that fills the region z < 0 behind the wavefront.
Although we were able to solve (68) and (69) numerically to find w m and ∆z m when n = 2, for larger values of n, the system is so ill-posed that we could not obtain a numerical solution. We were, however, able to find a numerical solution for n = 2, 3 and 4, but not for larger n, by taking a different approach that does not involve w m directly, and which also allows us to explain why the largest possible value of δ is appropriate in the far field solution. Instead of using (70) directly, we write where, for φ = Φ n , At leading order as c → 0, in the regions, z m+1 < z < z m = z m+1 + ∆z m between the spikes, we therefore havẽ In addition, (70) shows that L z is as smooth as φ = Φ n , so that We must also havẽ so that the solution matches with the spike regions in the neighbourhood of z = z m . This gives 2n + 2 conditions for the 2n + 2 unknowns, A m , a m k , b m k and ∆z m in each region between the spikes.
For z > z 1 = 0, we havẽ which satisfies eL 0 → 0 as z → ∞, and L 0 = dL 0 /dz = 0 at z = 0, and involves n − 1 unknown constants. In order to close this system we therefore need n − 1 conditions as m → ∞, so we must consider the far field solution.
As m → ∞, ∆z m → 0, so we rescale using so that (83) with c = 0 becomes At leading order as ∆z m → 0, Since L m = dL m /dz = 0 at z = 0 and z = 1, the solution is the polynomial It is also straightforward, using computer algebra, to linearise (91) and (96) about the equilibrium solution and find that each of the n − 1 equilibrium solutions has a stable manifold of a different dimension. In order to solve numerically, we truncate at large, finite m and need n − 1 boundary conditions, one of which is to specify ∆z m , so we need the solution with an (n − 2)-dimensional stable manifold. This, consistent with numerical solutions of the TW equation, is the solution with the largest value of δ. This gives us enough information to be able to solve the nonlinear system of equations given by (84) to (87) using fsolve in MATLAB. We used the Advanpix multiprecision computing toolbox, [13], to increase the precision of the computation, up to 66 digits of accuracy for n = 4. For n > 4 the system was too stiff for convergence to be achieved. The results, which are in excellent agreement with full numerical solutions of (5), are shown in figure 20, along with the far field behaviour indicated by ∆z m , w m = O(δ m ). We note that neither of the two numerical approaches described above is entirely satisfactory, and further work is required to devise a method better suited to this challenging system of equations. We have included these brief details of both methods here, since the former is more intuitive and leads naturally to the analytical results noted above, whilst the latter is more numerically stable.
Finally, we note that the case φ = Φ n is unusual since it has U ∞ = 1. A more typical solution, with a kernel that is a linear combination of the Gaussian, Φ ∞ , and Ψ 2 , is shown in figure 21. For z < z ∞ ≈ 17, U is oscillatory. This solution would be very hard to understand without the context provided by the preceding analysis.

Infinitely-differentiable kernels
When the kernel, φ, is infinitely-differentiable, the results of section 5, specifically that δ → 1 as n → ∞, indicate that z ∞ → −∞ as n → ∞, and that for φ = Φ ∞ and other infinitelydifferentiable kernels, the whole of the region behind the wavefront is filled with spikes in the limit c → 0. There is no obvious way to determine the large n asymptotic limit of the TW solution when φ = Φ n , specifically the solution of (89) when n −1 ∆z m 1, and the limit n → ∞ for (83) may not be well-defined, as discussed in [14].
Taking all of this into account, an appropriate ansatz for infinitely-differentiable kernels is and hence (98) After integrating (98), it is convenient to write the solution as with A a constant and Since I(z − z m ) → 0 as z → −∞, this means that w m ∼ ∆z m as m → ∞ and A = ∞ m=1 (w m − ∆z m ), so that Note that this is consistent with the result W ∞ ∼ ∆Z ∞ as n → ∞ (see figure 19). The conditions L =L z = 0 at each spike then lead to the infinite system of nonlinear equations for ∆z m and w m . Since there is no ordinary differential equation that determines L , solving (101) and (102) is the only way to investigate the behaviour of the spikes for infinitelydifferentiable kernels. In order to do this numerically, we need to quantify the behaviour of the far field solution. Although we have seen that ∆z m ∼ w m as m → ∞, it is not clear how to determine the individual behaviour of ∆z m and w m . We can, however, show that ∆z m cannot asymptote to a constant as m → ∞. For p 1, (101) becomes at leading order However, ∞ −∞ φ(z)dz = 1, so, for p 1, (103) is equivalent to a simple discretization of this integral, and if ∆z m → ∆z ∞ , a constant, as m → ∞, this is the trapezium rule, which is spectrally-accurate for integrals of analytic functions along the real axis, [15]. For example, A similar result holds for the smooth exponentially-and algebraically-decaying kernels that we study below, with the exponential rate at which G tends to zero controlled by the distance of the poles of φ from the real axis, [15]. Numerical calculation for these three kernels shows that in each case G is strictly positive, so that a far field solution of (101) cannot have a constant spacing of the spikes. We will see below that there is strong numerical evidence that ∆z m decays algebraically fast as m → ∞, but we cannot prove this. Fortunately, ∆z m ∼ w m as m → ∞ is enough information about the far field to allow us to truncate the infinite set of equations (101) and (102) and solve them numerically, although we will see that the crude approximation to the far field that we use introduces a significant error into the solution. We first reinstate the constant A defined in (99) and write (102) as We now note that Using a similar approach for (105), we truncate the infinite set of equations as (y)dy, at z = z p for p = 1, 2, . . . , M.
(107) This is a set of 2M equations in the 2M unknowns A, w m for m = 1, 2, . . . , M and ∆z m for m = 1, 2, . . . , M − 1. We solve the nonlinear system for A and ∆z m numerically in MATLAB using fsolve, solving the linear system for A and w m using linsolve. We begin with M = 2, for which a crude inital guess converges rapidly to a solution, and then successively increase M by one, with an initial guess of the solution formed by extending the vector ∆z m obtained for the previous value of M. Figure 22 shows the spacing and weight of the spikes for the smooth kernels Φ ∞ , Φ sech and Φ alg , calculated both from the full numerical solution of (5) and from the asymptotic solution, calculated numerically from (106) and (107). We find that the asymptotic solution reproduces the numerical solution, but that the inaccurate far field approximation discussed above leads to large errors for m > 1 2 M . In each case, the spacings that we are able to compute are not inconsistent with the region behind the wavefront being completely filled with spikes as c → 0, since it appears that 1 ∆z m m −1 as m → ∞. In addition, there is numerical evidence that for the Gaussian kernel, φ = Φ ∞ , an entire function, ∆z m = O(m −1/4 ) as m → ∞, whilst for the other two kernels, which have poles, at x = ±iπ for Φ exp and at x = ±i for Φ alg , ∆z m = O(m −1/2 ) as m → ∞. Moreover, we speculate that it is the position of the poles nearest to the real axis that controls this rate of decay, drawing on the results of [15], and as suggested by figure 22.

Conclusion
In this paper we have studied the structure of slow travelling wave solutions of the nonlocal Fisher equation, (1), when the diffusivity, D, and hence minimum wavespeed, c = c min ≡ 2 √ D, is small. For all the kernels that we used, the obvious common feature is the existence of one or more large, narrow spikes, with height of O(c −1/2 ) and width of O(c 1/2 ). The number and spacing of these spikes depend crucially on the behaviour of the kernel, φ(x), in the neighbourhood of x = 0. Kernels with a discontinuous derivative at x = 0 lead to travelling wave solutions with a finite number of spikes. Kernels that are differentiable, but not infinitely-differentiable, lead to an infinite number of spikes that fill a finite region behind the wavefront. Infinitelydifferentiable kernels lead to infinitely-many spikes that extend to infinity behind the wavefront. Although we were able to obtain detailed information about the solution for kernels that are not too smooth, the leading order problem for c 1 when the order of the discontinuous derivative of the kernel is greater than around eight is too poorly conditioned for the numerical methods that we developed to converge to a solution. For infinitely-differentiable kernels, we have been unable to determine analytically from the leading order equations any information about the far field behaviour of the spikes beyond w m ∼ ∆z m as m → ∞ (the weight and spacing of the spikes become identical in the far field), although numerical solutions of both the full travelling wave equation and the asymptotic solution for c 1 suggest that the weight and spacings of the spikes decay algebraically-fast as z → −∞, and we speculate that the position of the poles of the kernel in the complex plane may be the controlling factor. The results that we have obtained in this paper are summarized in table 1.
One interpretation of the solutions that we have studied here when c 1, using both numerical and asymptotic methods, is as nontrivial, steady state solutions of the nonlocal Fisher-KPP equation, (1), for D = 0, since they satisfy u (x) 1 − ∞ −∞ φ (x − y) u (y) dy = 0. The limit D → 0 is, as we have demonstrated, singular, but it is worth considering whether such steady state solutions could be relevant in physical situations for which the nonlocal Fisher-KPP is appropriate.
We now consider possible futher work related to the nonlocal Fisher-KPP equation that is suggested by the results presented in this paper. The use of fractional calculus to study travelling wave solutions when φ = Φ n for non-integer n, for example in (83), would be interesting, and could shed some light on the issue of convergence to a Gaussian kernel as n → ∞. Other extensions of this analysis could include kernels φ(x) that are non-monotonic, non-positive and/or discontinuous in x > 0; such kernels are used extensively in theoretical neuroscience, [16]. The extension of our analysis to two spatial dimensions is challenging, but would also be of interest.
Throughout the paper we have used asymptotic and numerical methods to gain insight into the travelling wave solutions. Our analysis lacks rigour, not only because of a lack of expertise, but also because the limit D → 0 is far from that where rigorous results are currently available, namely D → ∞, where the nonlocal Fisher-KPP equation asymptotes to the . The approximate algebraic rate of decay of ∆z m ∼ w m as m → ∞ suggested by these results in each case is shown by broken lines. standard, local Fisher-KPP equation, [9]. It seems clear that subtle features of the kernel control the gross features of the solution when D 1, and that the methods of modern nonlinear functional analysis could provide a route to a deeper understanding of this relationship.