Pushed and pulled fronts in a discrete reaction–diffusion equation

We consider the propagation of wave fronts connecting unstable and stable uniform solutions to a discrete reaction–diffusion equation on a one-dimensional integer lattice. The dependence of the wavespeed on the coupling strength μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} between lattice points and on a detuning parameter (a) appearing in a nonlinear forcing is investigated thoroughly. Via asymptotic and numerical studies, the speed both of ‘pulled’ fronts (whereby the wavespeed can be characterised by the linear behaviour at the leading edge of the wave) and of ‘pushed’ fronts (for which the nonlinear dynamics of the entire front determine the wavespeed) is investigated in detail. The asymptotic and numerical techniques employed complement each other in highlighting the transition between pushed and pulled fronts under variations of μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} and a.


Introduction
Mathematical analyses of a wide variety of physical and biological systems lead to investigation of spatially discrete nonlinear reaction-diffusion equations of the form on a discrete integer lattice with lattice points j ∈ Z at which u j = u j (t), the parameter μ > 0 dictating the coupling strength between lattice points while the constant a parameterises the nonlinearity (see below). Such models arise in, for example, the study of crystal growth [1] and binary alloy evolution [2] in material science and of neural networks for image processing and pattern recognition [3,4]; in addition, analogous formulations emerge naturally in the study of cell population growth [5] and cell signalling [6][7][8][9]. Lastly, such systems of course occur in the numerical solution by spatial discretisation of partial differential equations. We remark, however, that the formulations investigated in the above studies, and the equation analysed herein, are to be viewed as truly spatially discrete, and not as discretised versions of PDEs; indeed, the results that we present serve to illustrate how the behaviour of the continuous analogue of (1) relates to that of the discrete system, representing one of the simplest systems in which such a discrete-to-continuous transition can be explored. The study of spatially discrete diffusion equations has a long history, especially in the area of material sciences (see, e.g. Cook et al. [10] and references therein). An important aspect of such studies is in the choice of the nonlinearity f (u j ; a). For example, discrete diffusion equations with a bistable nonlinearity have been widely studied: see Keener [11], Elmer and van Vleck [12], Cahn et al. [13], Chow et al. [14], Fáth [15] and King and Chapman [16], and references therein, in which propagation failure was considered in detail, highlighting some key differences between discrete models and their continuous counterparts. Here, however, our focus is on the propagation of heteroclinic connections between unstable and stable states (rather than stable-stable connections), for which travelling-wave solutions to the continuous Fisher [or Kolmogorov-Petrovskii-Piskunov (KPP)] equation, in particular, have been extremely widely studied. In the discrete case, Zinner et al. [17] considered the existence of travelling-wave solutions to the discrete Fisher-KPP equation on a one-dimensional lattice, obtaining a constraint on the coupling strength for which strictly increasing travelling wavefronts of speed c > 0 exist; more recent work on related systems includes that of Guo and Morita [18] and Hakberg [19].
In this study, we consider for definiteness (and because of its widespread adoption as a model nonlinearity in the PDE case) equation (1) with the following cubic nonlinearity: and investigate the dependence on the parameters μ and a of the wave propagation speed through a discrete lattice. We emphasise that the specific form (2) is adopted for illustrative purposes; the features that we analyse below are applicable to a broad class of nonlinearities and much of the analysis applies to more general cases of smooth f with (1) scaled such that f (u; a) = u + o(1) as u → 0, f (1; a) = 0, f (u; a) > 0 for0 < u < 1.
The parameter a > 0 is often termed the 'detuning parameter' in the literature (e.g. Cahn et al. [13]). The uniform steady-state solutions of (1), (2) are u j ≡ 0 (unstable) and u j ≡ 1, −a (stable for a > 0). We remark that (1), (2) contains as limit cases both the (monostable) discrete Fisher equation (a → ∞), with steady states u j ≡ 0, 1, and the (bistable) discrete Newell-Whitehead-Segel equation with stable states u j ≡ ±1 (a = 1). Our interest here will be in fronts in which the stable state u j = 1 overruns the unstable one u j = 0. In the continuous counterpart of equation (1), μ corresponds to the diffusion coefficient; the study of travelling wave fronts in such continuous diffusion equations with a range of choices for the nonlinearity f (u; a) has been the subject of extensive investigation, in a wide variety of contexts, including (but not limited to) mathematical biology, plasma dynamics and pattern formation. We do not give a thorough review here, since our main focus is on the discrete equation; instead, we include a summary of the results of most relevance to the current work in Sect. 2. The emphasis of this work is to characterise, for the first time, the transition that occurs between the so-called 'pulled' (where the propagation speed, denoted c * , is characterised by the leading-edge behaviour) and 'pushed' (the entire nonlinear wave profile determines the wavespeed, c † ) fronts as the parameters μ and a are varied. While equivalent results are available for the continuum version of (1), (2), we are not aware of previous detailed analysis in discrete equations. Our results may be summarised briefly as follows. For μ = O(1), c * (μ) is straightforward to obtain; however, neither the transition a T (μ) nor c † (a, μ) are available analytically. In the limit a → 0 + , we obtain a new estimate for c † (a, μ). For μ 1, we isolate the transition as a T (μ) ∼ 1/ ln(1/μ) and provide new results for c † (a, μ) (a = O(μ)) and c † and c * at transition. The applicability of these results is illustrated by comparison with numerical simulation, and with the more well-known results that pertain for μ 1.
From the point of view of asymptotic methods, we additionally highlight the following features of the current work. A general remark is that discrete models have been far less studied from the point of view of formal asymptotics than their continuous counterparts. In contrast to the discrete bistable case (in which wave pinning can occur), the continuum limit μ → +∞ (addressed in Sect. 2) here does not lead to phenomena qualitatively absent in the continuous (PDE) case; asymptotic analysis of the weak-coupling limit μ → 0 + , analysed in Sect. 5 is, however, notably challenging and is correspondingly delicate from a numerical point of view (e.g. in that the transition from pushed to pulled fronts is particularly difficult to capture), making the asymptotic study particularly helpful in characterising the regimes in which both pulled-and pushed-front behaviours (see Sect. reftravellingspeeds) occur. The combination of Liouville-Green (JWKB) and matched-asymptotic approaches needed here is likely also to be applicable in a number of related contexts.
There is an issue that needs setting to one side before proceeding with the detailed analysis: we concern ourselves only with initial data that decay sufficiently rapidly such that the minimum wavespeed (denoted henceforth by c min ) of strictly positive travelling waves is realised. To characterise what 'sufficiently rapidly' means here, we note that (1), (3) linearised about u = 0, has travelling-wave solutions u j = exp(−λ( j − ct)), with Hence c(λ) takes its minimum value when and (4) both hold. When the initial data behave as exp(−λj) as j → +∞, where λ is less than the value determined by (4), (5) and the resulting c given by (4) has c > c min , then that value of c is realised as the large-time behaviour; for more rapidly decaying initial data, we are in the territory analysed in the remainder of the paper. The paper is organised as follows. In Sect. 2, we outline the two distinct classes (pulled and pushed) of waves of relevance here, and analyse the continuum (strong-coupling) limit μ → +∞. In Sect. 3, we set up the corresponding discrete travelling-wave problem and in Sect. 4, an asymptotic analysis of the transition between pulled and pushed waves in (1) is performed. Section 5 considers the weak-coupling limit (μ → 0 + ) in detail. In Sect. 6, we present numerical simulations indicating the key features of travelling-wave behaviour, and illustrating the applicability of the asymptotic results derived in previous sections. Section 7 includes a summary of the main results of this paper. Appendix 1 addresses the limit a → 0 + , Appendix 2 treats a linear differential-difference equation crucial to the results of Sect. 5, Appendix 3 gives an analysis of both stable-unstable and stable-stable connections for the case of a specific (explicitly solvable) nonlinear coupling between lattice points and Appendix 4 generalises aspects of the analysis of Sects. 2 and 4 to a much broader setting.

Pushed and pulled fronts and the continuum limit µ → +∞
This paper focusses on the propagation speed of large-time (travelling-wave) solutions to Eqs. (1), (2) in a onedimensional integer lattice j ∈ Z on which u j = u j (t). We consider in particular how these speeds differ from those obtained in the corresponding continuous reaction-diffusion equation.
The type of behaviour with which we are concerned is illustrated in Fig. 1, which shows a numerical solution indicating the evolution to two travelling waves (one propagating in each direction) from initial data given by (see Sect. 6 for details of the numerical procedures adopted) Our analysis will concern the rightward travelling wave, the other following from an obvious symmetry argument. Stable-unstable connections in systems such as (1) may be categorised into two distinct classes: in the case for which the propagation speed is determined by the leading edge of the wave, the front is termed 'pulled'; conversely, in the case of 'pushed' waves, the wavespeed is determined by the whole of the wavefront 1 . The pulled speed can be obtained by linearising the system around the homogeneous state (and is therefore also frequently referred to as the 'linearly selected' wavespeed); in this regime, the unstable homogeneous state is invaded by the stable one at a speed which tends to a constant c * as t → +∞; when nonlinear effects are dominant in a sense that will be clarified below, invading waves of speed c † > c * are obtained, and in such circumstances the linear analysis provides a lower bound for the pushed speed, but does not give the realised speed. Determining the class of wavefronts exhibited by a given discrete system is, in general, difficult (see, e.g. Plahte and Øyehaug [8]); in the continuous case (with a monostable nonlinearity), approaches to achieve this include that of Lucia et al. [22], but such methods do not seem likely to be applicable in the discrete case.
In the limit μ → +∞ in (1), neighbouring lattice points have (after an initial transient t = O(1/μ), unless the initial data are suitably prepared) almost equal values of u j , i.e. the solution is slowly varying with respect to j. The appropriate spatial variable is x = j/ √ μ; under this rescaling, with u j (t) ∼ u(x, t), we obtain Since the seminal studies of Fisher [23] and Kolmogorov et al. [24], travelling-wave solutions to (7) have been widely investigated; see, e.g. the review of van Saarloos [25].
Travelling-wave solutions of (7) take the form u(x, t) ∼ U (z), z = x − S(t), S(t) ∼ ct for constant 2 c and satisfying These describe the heteroclinic connections between the stable state u = 1 and the unstable one u = 0; note that for given c > 0, (8), (9) should be viewed as an initial value problem from z = −∞, the solution then being determined up to translations in z. Linearising (8) as z → +∞ gives (in view of (3)) so that Positivity (required by the comparison theorem for non-negative initial data) thus requires c ≥ 2. In the phase plane of (8), the origin is a degenerate stable node if c = 2, in which case, (10) implies and, when A > 0 holds for the solution to (8), (9), c = 2 is the wavespeed associated with a pulled front (i.e. c * = 2). As a decreases (in the case of (2), and more generally in (1), (3) if a is suitably defined), A reaches zero at a = a T , say, and for smaller a with 2 < c < c † for some c † (a), U (z) in (10) becomes negative at some z, approaching zero from below as z → +∞ (as can readily be demonstrated by a phase-plane analysis, a method that is not of course available in the discrete case); such non-monotonic waves are unstable and are in any case again precluded for non-negative initial data by the comparison theorem. In this case, pushed fronts occur, and we now further classify the two cases (summarising well-known results for (7)), whereby c min = c * for a > a T (pulled front) and c min = c † (a) > c * for 0 < a < a T (pushed front).
(a) Pulled fronts: a > a T , c min = 2. The wavespeed here follows either via (12) from the linearisation (10) of the travelling-wave ODE or from an application of the Liouville-Green (JWKB) method to the linearised PDE: (cf. Cuesta and King [26] and references therein), whereby gives (since φ x x does not contribute at leading order) where η = x/t. We therefore obtain giving the point at which linearisation becomes inapplicable (i.e. at which the solution (14) fails to be exponentially small) as η = 2 (implying S(t) ∼ 2t), thereby identifying the nonlinear wavefront as having c = 2.
Here the wavespeed is determined by the solution to (8) having (11), i.e. fast decay into the origin (the phase-plane analysis demonstrates that for a < a T this is the borderline between those c for which U (z) remains positive and those for which it crosses zero, whereas for a > a T this borderline is given by the degenerate-stable-node case (12)). For a < a T , the wavespeed c min can therefore be viewed as an eigenvalue (or second-kind similarity exponent), as is the wavespeed in the case of stable-stable connections; such an interpretation is not appropriate for pulled fronts. For a = a T , we have c min = 2, the value of a T being determined via the requirement that A = 0, B > 0 in (12), i.e. in this transition case a rather than c can be viewed as an eigenvalue. A specific reason for the choice of the nonlinearity (2) is that c † can be determined explicitly. Setting U = 1/V in (8) gives the quadratically nonlinear ODE: and established approaches to such problems suggest seeking a solution of the form: (the solution seems to have been identified for the first time by Hadeler and Rothe [20], without explicitly exploiting the quadratically nonlinear form). Equation (17) leads to two algebraic equations, thereby prescribing c as well as λ, as follows: Thus, in (11) we have and conversely for the negative square root. Hence (18), (20) corresponds to fast decay for a < 1/2 but slow decay for a > 1/2 (the latter will pertain for initial data that decay as exp(−x/ √ 2a) as x → +∞). We thus conclude that and we note that implying that c min (a) behaves rather smoothly as a passes through a T (cf. Sect. 4). In summary, for (7) we have (the latter in fact also holds for a < −1, but we shall limit ourselves to a > 0). For more general nonlinearities (3), c min = 2 for a ≥ a T remains true but c † (a) and a T each needs determining numerically.   (24) for 0 < a ≤ 1/2, the pushed speed for (7), showing good agreement for μ large and a lying in the range 10 −1 a 1, c, d c min given by (24) for a ≥ 1/2, the pulled speed for (7), showing good agreement for μ large and a sufficiently large, e, f c min given by (4), (5), the pulled speed for (1), showing good agreement for a sufficiently large and arbitrary μ. In each case, the range of agreement is as anticipated The surface and contour plots in Fig. 2 illustrate the range of applicability (by comparison with numerical simulations of (1); see Sect. 6 for details of the numerical procedures employed) of the above results in (μ, a) parameter space, namely the continuum wavespeeds (24) (rescaled by √ μ, as required for comparison with the discrete problem), appropriate for large μ (each in its own range of a-see Sect. 7), and the pulled speed determined by (4), (5), appropriate for a > a T (μ).

Travelling-wave analysis of the discrete system
We now set 3 in (1) to give the ordinary-differential-difference equation We note that while (1) holds for j ∈ Z, (26) applies for z ∈ R, in keeping with t ∈ R + applying in (1). Unlike the continuous case, this cannot be treated as an initial value problem from z = −∞; a boundary-condition count (respecting the invariance of (26) under translations of z) implies that (for given c) the required degrees of freedom as z → +∞ (of the form exp(−λz)) involve all the roots of (4) having positive real part. Pulled fronts have c * (μ) and λ * (μ) given by (4), (5) with λ * real, (5) being the repeated-root condition (cf. Zinner et al. [17], equation (6), the supremum in which is related to the minimum wavespeed statement that led to (5) above); (4) has two real roots for c > c * > 0 (these also being its roots with the smallest positive real part) and none for c < c * . Pushed fronts are those for which f (u; a) is such that there is a c = c † > c * for which the exponential corresponding to the smaller positive real root of (4) is absent as z → +∞ in the solution to (26), (27) (i.e. fast decay occurs). From (4), (5), we have that λ * (μ), c * (μ) are given by This implies the leading orders in which are consistent with the results of the previous section; for μ → 0 + , we have λ * → +∞ and with only terms exponentially smaller in λ * neglected, so that as μ → 0 + , These logarithmic dependencies provide advanced warning of the asymptotic challenges that will arise in Sect. 5.4.
In the discrete case, we have been unable to identify any f (u; a) for which c † can be determined explicitly for a pushed front (a somewhat abortive attempt to generalise (17) in this direction, exploiting a corresponding quadratically nonlinear form, is given in Appendix 3); explicit solutions can be constructed by specifying a suitable specific profile for U (z) and then determining the nonlinearity f from (26), but the resulting f will in general then depend on μ, so this procedure is not suitable for the current purposes. For example, proceeding in this way, we find that for the ansatz again applies with λ and c given by however, the dependence of (32) upon μ makes it of limited value to the current study.

Asymptotic analysis of the transition between pushed and pulled fronts
In this section, we analyse the behaviour close to the transition between pushed and pulled fronts-while we shall do so here in the framework of the ordinary-differential-difference equation (26), the criteria in question are of much more general applicability: see Appendix 4. In the current context, the analysis of this section will provide a diagnostic that is helpful in the analysis of the weak-coupling limit μ → 0 + . We require the solution to (26), defined uniquely for given wavespeed c up to translations in z -an indeterminacy that we shall eliminate by specifying U (0) = 1/2 -and (as usual) determine c through the requirements that U > 0 for z ∈ (−∞, ∞) and that U (z) exhibits a maximal decay rate as z → +∞.
For pulled fronts, having c = c * (μ), we denote U = U * (z; a). We reiterate that we denote by a = a T the value of a at which a transition between pulled fronts (for which a > a T without loss of generality 4 ) and pushed fronts (a < a T ) occurs. Then for a = a T we have (cf. (12)) 5 U * (z; a) ∼ (A(a)z + B(a))e −λ * z as z → +∞, the z prefactor being associated with λ = λ * being a repeated root. In (37), we have A > 0 for a > a T and A < 0 for a < a T . While only the former is realisable for non-negative initial data, U * will also play a role in the sequel for a < a T . At the transition, we have (cf. Cuesta and King [26]) with A † + > 0, and where U † (z; a) denotes the pushed-front profile, 7 having propagation speed c † . The key implications of the above analysis for the current section are as follows. Setting a = a T + δ with |δ| 1 we have from (37) that for δ > 0 while for δ < 0 it follows from (40) that, 8 retaining terms up to O(δ), as λ * → λ + , z → +∞ (as we shall see in (46), λ + − λ * = O(δ)). Now, returning to (26), we set as δ → 0, and it follows that in both the pulled (δ > 0) and pushed (δ < 0) cases. In view of (38), the matching condition (41) is automatically satisfied, while in (42) we infer from (35)-(36) that 6 Note that c † (a T , μ) = c * (μ). 7 The far-field expansion (39) will also in general contain a contribution of the form (40) with A † + replaced by A + (a, c), with A † + (a) = A + (a, c † ). 8 Significantly, the coefficient of the pre-exponential term linear in z is positive in (41) and negative in (42).
expresses the rather smooth transition of c min (a, μ) through a T . 5 The weak-coupling limit, µ → 0 +

Naïve time-dependent analysis
For reasons of notational convention, we set μ = ε in this section, and consider the behaviour of propagating fronts in (26), (27) for 0 < ε 1. It should be clear how the nonlinearity f (u; a) could be generalised; we shall again limit ourselves to the cubic form defined by (2).
Returning to the notation u j (t), consider, for definiteness, the following initial data: Then u j ∼ 1 for j ≤ 0 holds as ε → 0 for all time, and setting u j ∼ ε j v j for j > 0, t = O(1) gives and so Thus with leading order balance 9 with u 1 ∼ e t 1 as t 1 → −∞, and hence u 1 (t 1 ) is given by In (50), the final term in the summation dominates for t j, so application of Stirling's formula gives Setting j = ct, it is required for (54) to be of O(1) that which corresponds exactly to eliminating λ * in (30), as might be expected. Subsequent lattice points ( j > 1) are activated on timescales from which the wavespeed will be inferred via c(ε) = j/T j (ε) ln(1/ε); then and 10 u j ∼ e t j as t j → −∞, so that u j = u 1 (t j ). That the initial value problem (57) is identical to (52) corresponds to the solutions approaching a waveform of fixed profile and wavespeed Thus, once activated, each lattice point in turn activates its neighbour after a fixed timestep t ∼ ln(1/ε), corresponding to a simple cellular automaton. The analysis above gives little insight into whether the mechanism of wavespeed selection is of pulled or pushed type (even though (55) has the same asymptotic behaviour as c * (ε)). For example, (49) is independent of the nonlinearity while (57) is not; more substantially, the use above of (49) (the linearised behaviour ahead of the front) is consistent with pulled behaviour (and we shall find that this indeed occurs for a = O(1)), whereas the discrete-time and two-state (u j = 0 or u j = 1) cellular automaton just referred to would have u j ≡ 0 ahead of the wavefront, which might be interpreted as requiring pushed behaviour (which will arise below for a logarithmically small in ε). The transition between pulled and pushed fronts proves delicate to analyse and is therefore treated in some detail below (see Sect. 5.4).

Travelling-wave preliminaries
We now consider in more detail the small-μ asymptotics of travelling waves in the system (26), and their propagation speed. From (4), with μ replaced by ε, it follows that for c = O(1) the real roots λ + > λ − have so the faster decaying case, at least, exhibits very rapid decay in z in the limit ε → 0. It will prove useful in the sequel to introduce the slowness s = 1/c and the distinguished limit that includes the repeated-root (pulled-front) case corresponds to σ = O(1), where i.e.
Since λ ± are both large in this regime, (4) can be approximated (to all orders in ln(1/ε)) by i.e. setting λ = s + Λ, the repeated-root case (5) having Λ ∼ 1, σ ∼ 1/e. We now discuss the travelling-wave behaviour in the two regimes implicitly identified above, namely (i) a pushed case, c † = O(1) and (ii) σ = O(1), the former being significantly simpler but of less import for our purposes, in the sense that it resides firmly in the pushed-front regime.

a = O(ε)
Pushed fronts with c † = O(1) have a = O(ε). Here we first analyse an inner problem, which involves setting a = ε/α and z = εc † ξ to give at leading order for ξ = O(1) with U 0 → 1 as ξ → −∞, U 0 → 0 as ξ → +∞, and U 0 (0) = 1/2. On straightforward integration, the solution may be obtained in implicit form as with limiting behaviour Such algebraic decay is usually, but not here, associated for such nonlinearities with slowly decaying initial data generating waves travelling faster than the minimal speed. The outer scaling sets z = c † ζ , U = εV whereby 11 dV 0 dζ and, as ζ → 0 + , V 0 ∼ 1/αζ . The solution to (67) decays as e −ζ as ζ → ∞ for almost all s, corresponding to first of (59). To obtain a pushed front, we thus require that s be determined such that the leading-order solution satisfies Hence so that thereby providing explicitly the pushed speed (via c = 1/s) in terms of the detuning parameter in this regime (recall, a = ε/α).
As α → +∞ we have this limit corresponding to the case in which the nonlinearity is quadratic, rather than linear, as u → 0, in which case it is trivially clear that the front must be nonlinearly selected (i.e. pushed); see also Appendix 1. and Because V = o(1) for ζ > s, describing the transition to the exponential decay associated with λ + in (59) requires significantly more involved asymptotics, which we shall not pursue. In summary, we have shown here that a pushed front occurs for μ → 0 + , a = O(ε); setting a = ε/α the leading-order expressions for c † (a, μ) are A comparison of these results, and of the pulled ones (4), (5), with the speed observed in numerical simulations of (1), (2) is shown in Fig. 3; the agreement is excellent other than in an intermediate range of a that we address in the next subsection. We shall also in fact cover the regime a = O(1/ ln(1/ε)), which proves to be the most important of all, in this subsection. We first consider the pulled case with σ = O(1) in (60). Setting z = c * ζ , with (see (61)) gives Hence, setting we have so that, for the nonlinearity (2), U 0 is given in implicit form by with far-field behaviour where Moreover, so that U 1 ∼ −1 as ζ → +∞. For z = O(1) with 0 < z < 1, the contributions for which we need to account follow from the balance (which we continue to write in terms of ζ to make the matching into ζ = O(1) transparent) In the remaining regions, the inequalities U (ζ + s) provides the relevant balance, and we set 12 to give We can now appeal directly to the analysis of Appendix 2, which is devoted to (87) with z = ξ + 1. The initial data in 0 < ξ ≤ 1 follow 13 from (84); adopting (111) with ν = 0 and ν = s, and using (114), we then find that 14 It follows from (110) that λ − < 1, while s 1, so the sign of the coefficient in (88) is determined by that of 15 (using (60) and (110)) It remains to use (89), with K (a) defined by (81), to infer the value of σ . For given λ − (and hence σ ), if κ > 0 one has slow decay associated with the initial data having behaviour proportional to exp(−λ − j) as j → ∞ (with λ − here given by (4)). It is clear that κ is a decreasing function of λ − , the latter being maximal in the repeated-root 12 That the other terms from the central difference operator are negligible in (85) follows from the factors e −s arising from the transformation from U to V . 13 Because ζ s in (78), this fully nonlinear region makes no leading-order contribution to the integral in (109), though the contribution it does make is only logarithmically smaller: self-consistency checks incorporating such correction terms have been undertaken to confirm that they do not lead to the expansions derived in this section becoming invalid in the regimes considered. 14 The notation λ − is that of Appendix 2, not that above; the two usages of λ ± are in correspondence, however. case, i.e. at λ − = 1 . If κ(1; a) > 0, we infer that the front is pulled, with (114) being modified in the obvious way (cf. (37)) due to ρ = −1 being a second-order pole when σ = 1/e; however, if κ(1; a) < 0 the minimum-speed non-negative wave will be that for which corresponding to the wave being of pushed class. Since s ∼ ln(1/ε), the condition (90) requires that a 1, and hence, in view of (81) , K (a) = a + O(a 2 ) and 16 a = O(1/ ln(1/ε)). Hence, at transition, we have and the pushed fronts that occur for a < a T have (using (110)) and hence, from (61), so that It is clear from (91) and (93) that, to these orders, ds/da = 0 at a = a T , consistent with the analysis of Sect. 4 (recalling that c = 1/s). Moreover, formally setting a 1/ ln(1/ε) in (93) gives thereby matching successfully with the regime analysed in Sect. 5.3 (see (73)).
Much of the analysis of pushed fronts (notably in Appendix 2) here relies, unusually, on linearisation ahead of the wavefront; the nonlinear profile from (79) feeds through only in the small-a behaviour of K (a) in (81). We observe that the propagation speed is rather insensitive here to the values of ε and of a, and hence even to whether the front is pulled or pushed (compare (31) and (94)), making the numerical identification of the transition particularly challenging.

Numerical results
In this section, we outline our numerical approach and present further simulations complementing, and comparing the observed wave propagation speed in (1), (2) to, the asymptotic results provided above, thereby investigating a wide range of parameter values. We consider a domain j = 1 . . . N , and the system of N ODEs is solved using the initial value problem solver ode15s in MATLAB with Neumann-type boundary conditions; i.e. u 0 = u 2 , u N +1 = u N −1 . Initial conditions comprise a small region of the stable state u j (0) = 1, j ∈ [1, 50], adjacent to the unstable (trivial) steady state in the remainder of the domain. We define the wavefront position to be the lattice position j given by where t * j is the time at which each node attains the wavefront value, u j = 0.5. We remark that in order to obtain accurate wavespeed estimates, we employ a relatively large domain (we choose a lattice of 5000 nodes in our numerical simulations); however, the above numerical approach is significantly more accurate than the naïve method of calculating the speed directly from the position of the wavefront, which is adversely affected by the fact that the wave then moves in discrete jumps. Figure 4 illustrates the propagation speed c of a travelling wave observed in numerical simulations of (1), (2), under variation of the parameters a and μ. Figure 4a illustrates the dependence for various a. More instructively, Fig. 4b indicates two distinct regimes: for a sufficiently large, the wavespeed is insensitive to the value of a, in accord with expectations; for smaller a the wavespeed varies with a in correspondence with pushed-front behaviour.
Figures 5, 6, 7 and 9 explore further how the observed wavespeed compares with the asymptotic predictions. Figure 5 illustrates the applicability of the large-μ results. The deviation of the numerics from the pushed-front expression (24) for a small is noteworthy and is a consequence of (24) ceasing to be valid for μ = O(1/a), as captured by the analysis of Appendix 1. Figure 6 complements Fig. 3 in exploring in more detail the applicability of the asymptotic results of Sect. 5.3 for small μ, emphasising in particular, through the use of a linear scale, the relative portions of parameter space in which the different estimates apply. Figure 7 applies the results of Sect. Figs. 3 and 7 show deviation between numerics and asymptotics; given the logarithmic dependence of the asymptotics it is unsurprising that the quantitative agreement here is not good (numerically, reducing μ further would represent a significant challenge). The qualitative agreement is, however, reasonable given those issues and that the asymptotic expressions are derived on the basis that a = O(1/ ln(1/μ)): in Fig. 7 the pulled speed is given by the asymptotic expression (31) rather than the exact result used in Fig. 3, explaining the relatively poor agreement in the pulled regime-if both pulled and pushed asymptotic results are translated upwards in accordance with the exact pulled speed, the quantitative agreement is significantly improved. Figure 8 demonstrates the behaviour for μ = O(1), giving the comparison between numerical and analytical results (the latter are available when 0 < a < a T only for 0 < a 1).  pairs (a, μ) for which the numerically obtained speed deviates from that given by (4), (5) by less than 1 insensitivity of the wavespeed at small μ with a μ to whether the front is pulled or pushed, the quality of the agreement for small μ is unsurprising, but that for large μ is encouraging.  (1) with coupling strength and nonlinearity given by (115) and (2), respectively. Also shown are the results corresponding to the strong (dashed line) and weak-coupling (dotted line) limits (127) was on characterising, for the first time, the transition between so-called 'pulled' (where the propagation speed is characterised by the leading-edge behaviour) and 'pushed' (the entire nonlinear wave profile determines the wavespeed) fronts, under variation of the coupling strength μ and of the detuning parameter a that appears in the nonlinear term. While results characterising the transition between pushed and pulled fronts in the continuum version of the equation studied herein are well established (see e.g. Hadeler and Rothe [20], Stokes [21] and Rothe [27]), and the linearly selected speed of travelling waves in discrete systems is straightforward to obtain, we are not aware of previous detailed results for pushed waves in discrete equations.

in the intermediate regime in which
In summary, we have the following. For μ = O(1), neither a T (μ) nor c † (a, μ) are available analytically. However, the following general results hold.
The asymptotic results in terms of μ can be summarised as follows. Two distinguished limits arise when μ 1, in which case a T (μ) ∼ 1/ ln(1/μ).
See Figs. 3 and 6 for comparison of these results and of (97) to the numerical ones.
the latter expression for c * of course (because pulled speeds are independent of a) also holds for a ≥ O(1). Figure 7 is devoted to a comparison in this regime.
For μ 1, we instead have the following, with a T (μ) ∼ 1/2: whereĉ min is determined as in Appendix 1 via (104). (B) a = O(1) (with c scaled as in (26)); see Fig. 5 for a comparison. We note that in the weak-coupling limit analytical expressions are thus available for the wavespeed in both pulled and pushed regimes and, unlike those for the continuum limit, the latter can relatively easily be generalised to other forms of nonlinearity. It is important to stress that the pulled speed is available through the transcendental expression in (a) for arbitrary μ and is independent of a, depending only on the problem linearised about u = 0; by contrast, the pushed speed is much more challenging to obtain and the transition point a T (μ) is extremely difficult to get at, even numerically (particularly given that the transition between pulled and pushed occurs rather smoothly, as described in Sect. 4).
A number of conjectures warranting rigorous investigation, such as that a T (μ) is an increasing function of μ (contrast Appendix 3) arise naturally from this work. Similarly, obvious potential extensions come to mind, the higher-dimensional generalisation being the subject of current study.
The poles ofv thus occur at ρ = −λ with (cf. (63)), so for σ < 1/e there are two (real) roots for λ that we denote here by λ − and λ + , with 0 < λ − < λ + , and it is easy to show that the complex roots all have real part larger than λ + . Equation (107) requires initial data for 0 ≤ ξ ≤ 1 and for our purposes in Sect. 5 it suffices to consider the case v(ξ ) = e νξ for 0 ≤ ξ ≤ 1, for constant ν, and it then follows from (109) that If ν = −λ, with λ satisfying (110), we havê as is clear beforehand. More significantly for our purposes, the far-field behaviour follows from (112) as a residue contribution; suppressing such a slowly decaying term is a central ingredient in the selection mechanism for a pushed front in §5.4.
has μ ∼ μ in (115), the pulled wavespeed c * (μ) is again given by (4), (5). II Pushed waves Wave propagation speeds determined by the whole nonlinear wavefront may be obtained by the ansatz (which motivated (115) in the first place): We thereby obtain and hence with a pushed front occurring when λ † > λ * ; λ † (a, μ), and c † (a, μ) can be determined explicitly in the form: and the transition relationships λ † = λ * , c † = c * imply that a T (μ) is given by Equation (116) is evidently not of the class discussed in Appendix 3; moreover, setting in (116) with c = c * and W (z) slowly varying yields the dominant balance so the far-field behaviour differs from (37), and the behaviour outlined in Appendix 3 might not be expected to pertain. Nevertheless, it is easy to see that ∂c † (a T , μ)/∂a = 0 also holds in this case. Equation (121)  as a → 0 + , c † ∼ 2aμ as a → +∞, the former being consistent with the scaling argument embodied by (105) and the latter corresponding to the continuum limit. More importantly, and it follows that a T ∼ ln 1 2μ as μ → 0 + , a T ∼ 1 2 as μ → +∞; (127) the latter is as expected from the continuum limit, but the former implies a T → +∞ as μ → 0 + , i.e. the opposite behaviour from that in Sect. 5.4. Thus the analysis of (119) is counterproductive in terms of gaining insight into (26)-it does, however, reinforce the point that intuition into whether pushed or pulled behaviour is to be expected is hard to come by in the weakly coupled case. Figure 10 shows the curve (122) separating pushed and pulled waves, together with the weak-and strong-coupling limits (127). The offset of the latter for μ 1 illustrates further the implications of logarithmic terms in the associated asymptotic expansions (here through closed-form expressions, in contrast to §5.4).

Stable-stable connections
For completeness, we exploit the analytic tractability of (115) to address this case also. For a connection between the stable states U = 1, −a (stable for a > 0), representing the propagation of the state U = −a overrunning U = 1 (and vice versa, denotedŨ ), we set where V (z) is defined in (118). In each case, the wavespeed and decay rate may be obtained as in Sect. 3.1. The resulting wavespeeds are where c + corresponds to the wave U (z) and c − toŨ , and the decay rate λ(μ, a) is given in each case by λ(μ, a) = cosh −1 1 + (a + 1) 2 4aμ ; indeed the ostensibly distinct solution ansätze (128) in fact represent the same propagating front, moving in opposite directions.

Appendix 4: The generic pushed/pulled transition
Here we revisit briefly the analysis of Sect. 4 in a much more general setting. We consider the travelling-wave problem for a pseudo-differential operator with symbol P, whereby P d dz e −λz = P(−λ)e −λz (132) and the analysis henceforth will for the most part depend on P only in the form of the function P(−λ); for the discrete problem (1), P(−λ) is P(−λ) = 2μ(cosh(λ) − 1).
We again let f (u; a) = u + o(u) as u → 0 and assume P and f are such that for any given c > 0 a unique connection satisfying U → 1 as z → −∞, U → 0 as z → +∞, U (0) = 1 2 (135) exists. As z → +∞ we will in general have (39), where λ is a root of with smallest real part. The repeated-root case has c * and λ * determined by in which case U satisfies (37). The transition value a = a T is given by A(a T ) = 0 with, generically, B(a T ) > 0, and we again take the dependence of f upon a to be such that a pulled front arises for a < a T and a pushed one for a > a T . In the Liouville-Green approach (14), we have on the envelope solution to the second of these (i.e. the Clairaut equation) holds, so by (137) it follows that F = 0 on η = c * , dF/dη = λ * , as is to be anticipated. Correspondingly, the expansion-fan solution of the first of (138), parameterised by q = ∂φ/∂x (which is constant on rays) is and determining where φ = 0 gives a further derivation of the same result.