Nonlinear stability of pulse solutions for the discrete FitzHugh-Nagumo equation with infinite-range interactions

We establish the existence and nonlinear stability of travelling pulse solutions for the discrete FitzHugh-Nagumo equation with infinite-range interactions close to the continuum limit. For the verification of the spectral properties, we need to study a functional differential equation of mixed type (MFDE) with unbounded shifts. We avoid the use of exponential dichotomies and phase spaces, by building on a technique developed by Bates, Chen and Chmaj for the discrete Nagumo equation. This allows us to transfer several crucial Fredholm properties from the PDE setting to our discrete setting.


Introduction
The FitzHugh-Nagumo partial differential equation (PDE) is given by where g(·; r 0 ) is the cubic bistable nonlinearity given by g(u; r 0 ) = u(1 − u)(u − r 0 ) (1.2) and ρ, γ are positive constants. This PDE is commonly used as a simplification of the Hodgkin-Huxley equations, which describe the propagation of signals through nerve fibres. The spatially homogeneous version of this equation was first stated by FitzHugh in 1961 [25] in order to describe the potential felt at a single point along a nerve axon as a signal travels by. A few years later [26], the diffusion term in (1.1) was added to describe the dynamics of the full nerve axon instead of just a single point. As early as 1968 [27], FitzHugh released a computer animation based on numerical simulations of (1.1). This video clip clearly shows that (1.1) admits isolated pulse solutions resembling the spike signals that were measured by Hodgkin and Huxley in the nerve fibres of giant squids [32]. As a consequence of this rich behaviour and the relative simplicity of its structure, (1.1) has served as a prototype for several similar systems. For example, memory devices have been designed using a planar version of (1.1), which supports stable stationary radially symmetric spot patterns [43]. In addition, gas discharges have been described using a three-component FitzHugh-Nagumo system [49,54], for which it is possible to find stable travelling spots [59].
Mathematically, it turned out to be a major challenge to control the interplay between the excitation and recovery dynamics and rigorously construct the travelling pulses visualized by FitzHugh in [27]. Such pulse solutions have the form (u, w)(x, t) = (u 0 , w 0 )(x + c 0 t), (1.3) in which c 0 is the wavespeed and the wave profile (u 0 , w 0 ) satisfies the limits lim |ξ|→∞ (u 0 , w 0 )(ξ) = 0. (1.4) Plugging this Ansatz into (1.1) and writing ξ = x + c 0 t, we see that the profiles are homoclinic solutions to the travelling wave ordinary differential equation (ODE) c 0 u 0 (ξ) = u 0 (ξ) + g(u 0 (ξ); r 0 ) − w 0 (ξ) c 0 w 0 (ξ) = ρ u 0 (ξ) − γw 0 (ξ) . (1.5) The analysis of this equation in the singular limit ρ ↓ 0 led to the birth of many techniques in geometric singular perturbation theory, see for example [40] for an interesting overview. Indeed, the early works [10,31,41,42] used geometric techniques such as the Conley index, exchange lemmas and differential forms to construct pulses and analyse their stability. A more analytic approach was later developed in [46], where Lin's method was used in the r 0 ≈ 1 2 regime to connect a branch of so-called slow-pulse solutions to (1.5) to a branch of fast-pulse solutions. This equation is still under active investigation, see for example [11,12], where the birth of oscillating tails for the pulse solutions is described as the unstable root r 0 of the nonlinearity g moves towards the stable root at zero.
Many physical, chemical and biological systems have an inherent discrete structure that strongly influences their dynamical behaviour. In such settings lattice differential equations (LDEs), i.e. differential equations where the spatial variable can only take values on a lattice such as Z n , are the natural replacements for PDEs, see for example [1,36,48]. Although mathematically it is a relatively young field of interest, LDEs have already appeared frequently in the more applied literature. For example, they have been used to describe phase transitions in Ising models [1], crystal growth in materials [9] and phase mixing in martensitic structures [58].
To illustrate these points, let us return to the nerve axon described above and reconsider the propagation of electrical signals through nerve fibres. It is well known that electrical signals can only travel at adequate speeds if the nerve fibre is insulated by a myeline coating. This coating admits regularly spaced gaps at the so-called nodes of Ranvier [51]. Through a process called saltatory conduction, it turns out that excitations of nerves effectively jump from one node to the next [47]. Exploiting this fact, it is possible [45] to model this jumping process with the discrete FitzHugh-Nagumo LDEu j = 1 h 2 (u j+1 + u j−1 − 2u j ) + g(u j ; r 0 ) − w j w j = ρ[u j − γw j ]. (1.6) The variable u j now represents the potential at the j th node, while the variable w j denotes a recovery component. The nonlinearity g describes the ionic interactions. Note that this equation arises directly from the FitzHugh-Nagumo PDE upon taking the nearest-neighbour discretisation of the Laplacian on a grid with spacing h > 0. Inspired by the procedure for partial differential equations, one can substitute a travelling pulse Ansatz (u j , w j )(t) = (u h , w h )(hj + c h t) (1.7) into (1.6). Instead of an ODE, we obtain the system (1. 8) in which ξ = hj + c h t. Such equations are called functional differential equations of mixed type (MFDEs).
In [35,36], Hupkes and Sandstede studied (1.6) and showed that for small values of ρ and r 0 sufficiently far from 1 2 , there exists a locally unique travelling pulse solution of this system and that it is asymptotically stable with an asymptotic phase shift. No restrictions were required on the discretisation distance h, but the results relied heavily on the existence of exponential dichotomies for MFDEs. As a consequence, the techniques developed in [35,36] can only be used if the discretisation involves finitely many neighbours. Such discretisation schemes are said to have finite range.
Recently, an active interest has arisen in non-local equations that feature infinite-range interactions. For example, Ising models have been used to describe the infinite-range interactions between magnetic spins arranged on a grid [1]. In addition, many physical systems, such as amorphous semiconductors [29] and liquid crystals [17], feature non-standard diffusion processes, which are generated by fractional Laplacians. Such operators are intrinsically non-local and hence often require infinite-range discretisation schemes [16].
Our primary interest here, however, comes from so-called neural field models, which aim to describe the dynamics of large networks of neurons. These neurons interact with each other by exchanging signals across long distances through their interconnecting nerve axons [7,8,50,57]. It is of course a major challenge to find effective equations to describe such complex interactions. One model that has been proposed [7,Eq. (3.31)] features a FitzHugh-Nagumo type system with infinite-range interactions.
Motivated by the above, we consider a class of infinite-range FitzHugh-Nagumo LDEs that includes the prototypė e −k 2 [u j+k + u j−k − 2u j ] + g(u j ; r 0 ) − w j w j = ρ[u j − γw j ], (1.9) in which κ > 0 is a normalisation constant. In [22], Faye and Scheel studied equations such as (1.9) for discretisations with infinite-range interactions featuring exponential decay in the coupling strength.
They circumvented the need to use a state space as in [35], which enabled them to construct pulses to (1.9) for arbitrary discretisation distance h. Very recently [23], they developed a center manifold approach that allows bifurcation results to be obtained for neural field equations.
In this paper, we also construct pulse solutions to equations such as (1.9), but under weaker assumptions on the decay rate of the couplings. Moreover, we will establish the nonlinear stability of these pulse solutions, provided the coupling strength decays exponentially. However, both results do require the discretisation distance h to be very small.
In particular, we will be working in the continuum limit. The pulses we construct can be seen as perturbations of the travelling pulse solution of the FitzHugh-Nagumo PDE. However, we will see that the travelling wave equations are highly singular perturbations of (1.5), which poses a significant mathematical challenge. On the other hand, we do not need to use exponential dichotomies directly in our non-local setting as in [36]. Instead we are able to exploit the detailed knowledge that has been obtained using these techniques for the pulses in the PDE setting.
Our approach to tackle the difficulties arising from this singular perturbation is strongly inspired by the work of Bates, Chen and Chmaj. Indeed, in their excellent paper [1], they study a class of systems that includes the infinite-range discrete Nagumo equatioṅ e −k 2 [u j+k + u j−k − 2u j ] + g(u j ; r 0 ), (1.10) in which κ > 0 is a normalisation constant. This equation can be seen as a discretisation of the Nagumo PDE u t = u xx + g(u; r 0 ). (1.11) The authors show that, under some natural assumptions, these systems admit travelling pulse solutions for h small enough.
In the remainder of this introduction we outline their approach and discuss our modifications, which significantly broaden the application range of these methods. We discuss these modifications for the prototype (1.9), but naturally they can be applied to a broad class of systems.
The key contribution in [1] is that the authors fix a constant δ > 0 and use the invertibility of L 0;u0:sc;c0:sc + δ to show that also L h;u0:sc;c0:sc + δ is invertible. In particular, they consider weaklyconverging sequences {v n } and {w n } with (L h;u0:sc;c0:sc + δ)v n = w n and try to find a uniform (in δ and h) upper bound for the L 2 -norm of v n in terms of the L 2 -norm of w n . Such a bound is required to rule out the limitless transfer of energy into oscillatory modes, a key complication when taking weak limits. To obtain this bound, the authors exploit the bistable structure of the nonlinearity g to control the behaviour at ±∞. This allows the local L 2 -norm of v n on a compact set to be uniformly bounded away from zero. Since the operator L h;u0:sc;c0:sc + δ is not self-adjoint, this procedure must be repeated for the adjoint operator.
Transfer of Fredholm properties: system case. Plugging the travelling pulse Ansatz (u, w) j (t) = (u h , w h )(hj + c h t) (1.15) into (1.9) and writing ξ = hj + c h t, we see that the profiles are homoclinic solutions to the equation e −k 2 u h (ξ + kh) + u h (ξ − kh) − 2u h (ξ) + g(u h (ξ); r 0 ) − w h (ξ) (1. 16) We start by considering the linearised operator K h;u0;c0 of the system (1.16) around the pulse solution (u 0 , w 0 ) of the FitzHugh-Nagumo PDE with wavespeed c 0 . This operator is given by K h;u0;c0 v w (ξ) = L h;u0;c0 v(ξ) + w(ξ) c 0 w (ξ) − ρv(ξ) + ργw(ξ) , (1.17) where L h;u0;c0 is given by equation (1.12), but with u 0:sc replaced by u 0 and c 0:sc by c 0 . In §3 we use a Fredholm alternative as described above to establish the invertibility of K h;u0;c0 +δ for fixed δ > 0. However, the transition from a scalar equation to a system is far from trivial. Indeed, when transferring the Fredholm properties there are multiple cross terms that need to be controlled. We are aided here by the relative simplicity of the terms in the equation that involve w. In particular, three of the four matrix-elements of the linearisation (1.17) have constant coefficients. We emphasize that it is not sufficient to merely assume that the limiting state (0, 0) is a stable equilibrium of (1.9). In [56] we explore a number of structural conditions that allow these types of arguments to be extended to general multi-component systems.
Construction of pulses. Using these results for K h;u0;c0 a standard fixed point argument can be used to show that the system (1.9) has a locally unique travelling pulse solution (U h (t)) j = (u h , w h )(hj +c h t) for h small enough, which converges to a travelling pulse solution of the FitzHugh-Nagumo PDE as h ↓ 0. Indeed, one can mimic the approach developed in [1, §4], which in turn closely follows the lines of a standard proof of the implicit function theorem.
Spectral stability. The natural next step is to study the linear operator K h;u h ;c h that arises after linearising the system (1.9) around its new-found pulse solution. This operator is given by where L h;u h is given by equation (1.12), but with u 0:sc replaced by u h and c 0:sc by c h . The procedure above can be repeated to show that for fixed δ > 0, it also holds that K h;u h ;c h + δ is invertible for h small enough. However, to understand the spectral stability of the pulse, we need to consider the eigenvalue problem for fixed values of h and λ ranging throughout a half-plane. Switching between these two points of view turns out to be a delicate task. We start in §4 by showing that K h;u h ;c h and its adjoint K * h;u h ;c h are Fredholm operators with one-dimensional kernels. This is achieved by explicitly constructing a kernel element for K * h;u h ;c h that converges to a kernel element of the adjoint of the operator corresponding to the linearised PDE. An abstract perturbation argument then yields the result.
In particular, we see that λ = 0 is a simple eigenvalue of K h;u h ;c h . In §5 we establish that in a suitable half-plane, the spectrum of this operator consists precisely of the points {k2πic h 1 h : k ∈ Z}, which are all simple eigenvalues. We do this by first showing that the spectrum is invariant under the operation λ → λ + 2πic h h , which allows us to restrict ourselves to values of λ with imaginary part in between − πc h h and πc h h . Note that the period of the spectrum is dependent on h and grows to infinity as h ↓ 0. This is not too surprising, since the spectrum of the linearisation of the PDE around its pulse solution is not periodic. However, this means that we cannot restrict ourselves to a fixed compact subset of the complex plane for all values of h at the same time. In fact, it takes quite some effort to keep the part of the spectrum with large imaginary part under control.
It turns out to be convenient to partition our 'half-strip' into four parts and to calculate the spectrum in each part using different methods. Values close to 0 are analysed using the Fredholm properties of K h;u h ;c h exploiting many of the results from §4; values with a large real part are considered using standard norm estimates, but values with a large imaginary part are treated using a Fourier transform. The final set to consider is a compact set that is independent of h and bounded away from the origin. This allows us to apply a modified version of the procedure described above that exploits the absence of spectrum in this region for the FitzHugh-Nagumo PDE.
Let us emphasize that our arguments here for bounded values of the spectral parameter λ strongly use the fact that the PDE pulse is spectrally stable. The main complication to establish the latter fact is the presence of a secondary eigenvalue that is O(ρ)-close to the origin. Intuitvely, this eigenvalue arises as a consequence of the interaction between the front and back solution to the Nagumo equation that are both part of the singular pulse that arises in the ρ ↓ 0 limit. In the PDE case, Jones [41] and Yanagida [60] essentially used shooting arguments to construct and analyze an Evans function E(λ) that vanishes precisely at eigenvalues. In particular, they computed the sign of E (0) and used counting arguments to show that the secondary eigenvalue discussed above lies to the left of the origin. Currently, a program is underway to build a general framework in this spirit based on the Maslov index [2,13,33], which also works in multi-dimensional spatial settings. In [18,19] this framework was applied to an equal-diffusion version of the FitzHugh-Nagumo PDE.
An alternative approach involving Lin's method and exponential dichotomies was pioneered in [46]. Based upon these ideas stability results have been obtained for the LDE (1.6) [36] and the PDE (1.1) [11] in the non-hyperbolic regime r 0 ∼ 0. The first major advantage of this approach is that explicit bifurcation equations can be formulated that allow asymptotic expansions to be developed for the location of the interaction eigenvalue discussed above. The second major advantage is that it allows us to avoid the use of the Evans function, which cannot easily be defined in discrete settings because MFDEs are ill-posed as initial value problems [52]. We believe that a direct approach along these lines should also be possible for the infinite range system (1.9) as soon as exponential dichotomies are available in this setting.
Nonlinear stability. The final step in our program is to leverage the spectral stability results to obtain a nonlinear stability result. To do so, we follow [36] and derive a formula that links the pointwise Green's function of our general problem (1.9) to resolvents of the operator K h;u h ;c h in §6. Since we have already analysed the latter operator in detail, we readily obtain a spectral decomposition of this Green's function into an explicit neutral part and a residual that decays exponentially in time and space. Therefore we obtain detailed estimates on the decay rate of the Green's function for the general problem. These Green's functions allow us to establish the nonlinear stability of the family of travelling pulse solutions U h . To be more precise, for each initial condition close to U h (0), we show that the solution with that initial condition converges at an exponential rate to the solution U h (· +θ) for a small (and unique) phase shiftθ.
We emphasize that by now there are several techniques available to obtain nonlinear stability results in the relatively simple spectral setting encountered in this paper. If a comparison principle is available, which is not the case for the FitzHugh-Nagumo system, one can follow the classic approach developed by Fife and McLeod [24] to show that travelling waves have a large basin of attraction. Indeed, one can construct explicit sub-and super-solutions that trade-off additive perturbations at t = 0 to phase-shifts at t = ∞. In fact, one can actually use this type of argument to establish the existence of travelling waves by letting an appropriate initial condition evolve and tracking its asymptotic behaviour [14,37]. For systems that can be written as gradient flows, which is also not the case here, the existence and stability of travelling waves can be obtained by using an elegant variational technique that was developed by Gallay and Risler [28].
In the spatially continuous setting, it is possible to freeze a travelling wave by passing to a comoving frame. In our setting, one can achieve this by simply adding a convective term −c 0 ∂ x (u, w) to the right hand side of (1.1). The main advantage is that one can immediately use the semigroup exp[tL 0 ] to describe the evolution of the linearised system in this co-moving frame, which is temporally autonomous. Here L 0 is the standard linear operator associated to the linearisation of (1.1) around (u 0 , w 0 ); see (2.10). For each ϑ ∈ R one can subsequently construct the stable manifold of u 0 (· + ϑ), w 0 (· + ϑ) by applying a fixed point argument to Duhamel's formula. Upon varying ϑ, these stable manifolds span a tubular neighbourhood of the family (u 0 , w 0 )(· + R). This readily leads to the desired stability result; see e.g. [44, §4]. We remark here that these stable manifolds are all related to each other via spatial shifts.
In the spatially discrete setting the wave can no longer be frozen. In particular, the linearisation of (1.6) around the pulse (1.7) leads to an equation that is temporally shift-periodic. In [15] the authors attack this problem head-on by developing a shift-periodic version of Floquet theory that leads to a nonlinear stability result in ∞ . However, they delicately exploit the geometric structure of ∞ and it is not clear how more degenerate spectral pictures can be fitted into the framework. These issues are explained in detail in [36, §2].
In [5] the authors found a way to express the Green's function of the temporally shift-periodic linear discrete equation in terms of resolvents of the linear operator L h associated to the pulse (1.7). Based on this procedure, it is possible to follow the spirit of the powerful pointwise Green's function techniques pioneered by Zumbrun and Howard [62]. Indeed, in [3] a stability result is obtained in the setting of discrete conservation laws, where one encounters curves of essential spectrum that touch the imaginary axis. Using exponential dichotomies in a setting with extended state-spaces This allowed the techniques from [4] to be transferred from the continuous to the discrete setting. A slightly more streamlined approach was developed in [36], which does not need the extended statespace and avoids the use of a variation-of-constants formula. However, exponential dichotomies are still used at certain key points.
In our paper we follow the spirit of the latter approach and extend it to the present setting with infinite-range interactions. In particular, we show how the use of exponential dichotomies can be eliminated all together, which is a delicate task. In addition, we need to be very careful in many computations since integrals and sums over shifts as in (1.16) can no longer be freely exchanged. We emphasize here that our techniques do not depend on the specific LDE that we are analyzing. All that is required is the spectral setting described above and the fact that the shifts appearing in the problem are all rationally related.
Let us mention that it is also possible to bypass the construction of the stable manifolds altogether and employ a direct phase-tracking approach along the lines of [61]. In particular, one can couple the system with an extra equation for the phase. To close the system, one chooses this extra equation in such a way that the resulting nonlinear terms never encounter the non-decaying part of the relevant semigroup. Such an approach has been used in the current spectral setting to show that travelling waves remain stable under the influence of a small stochastic noise term [30].
Acknowledgements. Both authors acknowledge support from the Netherlands Organization for Scientific Research (NWO) (grant 639.032.612).

Main results
We consider the following system of equationṡ which we refer to as the (spatially) discrete FitzHugh-Nagumo equation with infinite-range interactions. Often, for example in [35,36], it is assumed that only finitely many of these coefficients α k are non-zero. However, we will impose the following much weaker conditions here.
We note that (2.4) is automatically satisfied if α 1 > 0 and α k ≥ 0 for all k ∈ Z >1 . The conditions in (Hα1) ensure that for φ ∈ L ∞ (R) with φ ∈ L 2 (R), we have the limit see Lemma 3.5. In particular, we can see (2.1) as the discretisation of the FitzHugh-Nagumo PDE (1.1) on a grid with distance h. Additional remarks concerning the assumption (Hα1) can be found in [1, §1]. Throughout this paper, we impose the following standard assumptions on the remaining parameters in (2.1). The last condition on γ in (HS) ensures that the origin is the only j-independent equilibrium of (2.1).
Without explicitly mentioning it, we will allow all constants in this work to depend on r 0 , ρ and γ. Dependence on h will always be mentioned explicitly. We will mainly work on the Sobolev spaces with their standard norms (2.7) Our goal is to construct pulse solutions of (2.1) as small perturbations to the travelling pulse solutions of the FitzHugh-Nagumo PDE. These latter pulses satisfy the system with the boundary conditions If (u 0 , w 0 ) is a solution of (2.8) with wavespeed c 0 , then the linearisation of (2.8) around this solution is characterized by the operator L 0 : (2.10) The existence of such pulse solutions for the case when ρ is close to 0 is established in [40, §5.3].
Here, we do not require ρ > 0 to be small, but we simply impose the following conditions. Assumption (HP1). There exists a solution (u 0 , w 0 ) of (2.8) that satisfies the conditions (2.9) and has wavespeed c 0 = 0. Furthermore the operator L 0 is Fredholm with index zero and it has a simple eigenvalue in zero.
Recall that an eigenvalue λ of a Fredholm operator L is said to be simple if the kernel of L − λ is spanned by one vector v and the equation (L − λ)w = v does not have a solution w. Note that if L has a formal adjoint L * , this is equivalent to the condition that v, w = 0 for all non-trivial w ∈ ker(L * − λ).
We note that the conditions on L 0 formulated in (HP1) were established in [41] for small ρ > 0. In addition, these conditions imply that u 0 and w 0 decay exponentially.
In order to find travelling pulse solutions of (2.1), we substitute the Ansatz into (2.1) to obtain the system (2.12) in which ξ = hj + c h t. The boundary conditions are given by The existence of such solutions is established in our first main theorem.
Note that this result is very similar to [22, Corollary 2.1]. However, Faye and Scheel impose an extra assumption, similar to (Hα2) below, which we do not need in our proof. This is a direct consequence of the strength of the method from [1] that we described in §1.
Building on the existence of the travelling pulse solution, the natural next step is to show that our new-found pulse is asymptotically stable. However, we now do need to impose an extra condition on the coefficients {α k } k>0 .
Note that the prototype equation (1.9) indeed satisfies both assumptions (Hα1) and (Hα2). An example of a system which satisfies (Hα1), but not (Hα2) is given bẏ (2.16) in which κ = 6 π 2 is the normalisation constant. Moreover we need to impose an extra condition on the operator L 0 given by (2.10). This spectral stability condition is established in [20,Theorem 2] together with [60, Theorem 3.1] for the case where ρ is close to 0. Assumption (HP2). There exists a constant λ * > 0 such that for each λ ∈ C with Re λ ≥ −λ * and λ = 0, the operator is invertible.
To determine if the pulse solution described in Theorem 2.1 is nonlinearly stable, we must first linearise (2.12) around this pulse and determine the spectral stability. The linearised operator now takes the form Here the operator ∆ h is given by As usual, we define the spectrum, σ(L), of a bounded linear operator L : Our second main theorem describes the spectrum of this operator L h , or rather of −L h , in a suitable half-plane. We emphasize that λ 3 does not depend on h. The translational invariance of (2.12) guarantees that λ = 0 is an eigenvalue of −L h . In Lemma 5.1 we show that the spectrum of the operator L h is periodic with period 2πic h 1 h , which means that the eigenvalues k2πic h 1 h for k ∈ Z all have the same properties as the zero eigenvalue.
Our final result concerns the nonlinear stability of our pulse solution, which we represent with the shorthand The perturbations are measured in the spaces p , which are defined by . Assume that (HP1),(HP2), (HS), (Hα1) and (Hα2) are satisfied. Fix 0 < h ≤ h * * and 1 ≤ p ≤ ∞. Then there exist constants δ > 0, C > 0 and β > 0, which may depend on h but not on p, such that for all initial conditions U 0 ∈ p with U 0 − U h (0) p < δ, there exists an asymptotic phase shiftθ ∈ R such that the solution U = (u, w) of (2.1) with U (0) = U 0 satisfies the bound for all t > 0.

The singular perturbation
The main difficulty in analysing the travelling wave MFDE (2.12) is that it is a singular perturbation of the ODE (2.8). Indeed, the second derivative in (2.8) is replaced by the linear operator ∆ h : We will see in Lemma 3.5 that for all φ ∈ L ∞ (R) with φ ∈ L 2 (R), we have that lim Hence the bounded operator ∆ h converges pointwise on a dense subset of H 1 (R) to an unbounded operator on that same dense subset. In particular, the norm of the operator ∆ h grows to infinity as h ↓ 0.
Since there are no second derivatives involved in (2.12), we have to view it as an equation posed on the space H 1 (R) × H 1 (R), while the ODE (2.8) is posed on the space H 2 (R) × H 1 (R). From now on we write (3. 2) The main results in this section will be used in several different settings. In order to accommodate this, we introduce the following conditions.
Assumption (hM). The set M ⊂ C is compact with 0 / ∈ M . In addition, recalling the constant λ * appearing in (HP2), we have Re z ≥ −λ * for all z ∈ M .
In the proof of Theorem 2.1 we choose (ũ h ,w h ) andc 0 to be (u 0 , w 0 ) and c 0 for all values of h. However, in §4 we let (ũ h ,w h ) be the travelling pulse (u h , w h ) from Theorem 2.1 and we letc h be its wave speed c h .
If (hFam) is satisfied, then for δ > 0 and h > 0 we define the operators These operators are bounded linear functions from H 1 to L 2 . We see that L − h,δ is the adjoint operator of L + h,δ , in the sense that holds for all (φ, ψ), (θ, χ) ∈ L 2 . Here we have introduced the notation for (φ, ψ), (θ, χ) ∈ L 2 .
Since at some point we want to consider complex-valued functions, we also work in the spaces H 2 C (R), H 1 C (R) and L 2 C (R), which are given by (3.7) These spaces are equipped with the inner product As before, we write (3.9) Each operator L from H 1 to L 2 can be extended to an operator from H 1 C to L 2 C by writing It is well-known that this complexification preserves adjoints, invertibility, inverses, injectivity, surjectivity and boundedness, see for example [53]. If λ ∈ C then the operators L ± h,λ are defined analogously to their real counterparts, but now we view them as operators from H 1 Whenever it is clear that we are working in the complex setting we drop the subscript C from the spaces H 1 C and L 2 C and simply write H 1 and L 2 . We also introduce the operators L ± 0 : H 2 (R) × H 1 (R) → L 2 (R) × L 2 (R), that act as and These operators can be viewed as the formal h ↓ 0 limits of the operators L ± h,0 . Upon introducing the notation we see that L + 0 (φ + 0 , ψ + 0 ) = 0 by differentiating (2.8). To set the stage, we summarize several basic properties of L ± 0 . For completeness, we provide a proof of these facts later on in this section.

For every
4. There exists a positive constant C 1 such that There exists a positive constant C 2 and a small constant δ 0 > 0 such that for all 0 < δ < δ 0 we have 6. If (HP2) is also satisfied, then for each M ⊂ C that satisfies (hM), there exists a constant C 3 > 0 such that the uniform bound The main goal of this section is to prove the following two propositions, which transfer parts (5) and (6) of Lemma 3.1 to the discrete setting.
Proposition 3.2. Assume that (hFam), (HP1), (HS) and (Hα1) are satisfied. There exists a positive constant C 0 and a positive function h 0 (·) : R + → R + , depending only on the choice of (ũ h ,w h ) and c h , such that for every 0 < δ < δ 0 and every h ∈ (0, h 0 (δ)), the operators L ± h,δ are homeomorphisms from H 1 to L 2 that satisfy the bounds Our techniques here are inspired strongly by the approach developed in [1, §2-4]. Indeed, Proposition 3.2 and Proposition 3.4 are the equivalents of [1, Theorem 4] and [1, Lemma 6] respectively. The difference between our results and those in [1] is that Bates, Chen and Chmaj study the discrete Nagumo equation, which can be seen as the one-dimensional fast component of the FitzHugh-Nagumo equation by setting ρ = 0 in (2.1). In addition, the results in [1] are restricted to λ ∈ R, while we allow λ ∈ C in Proposition 3.3. These differences play a crucial role in the proof of Lemma 3.10 below.
Recall the constant δ 0 > 0 appearing in Lemma 3.1. For 0 < δ < δ 0 and h > 0 we define the quantities Similarly for M ⊂ C that satisfies (hM) and h > 0 we define together with The key ingredients that we need to establish Propositions 3.2 and 3.3 are lower bounds on the quantities Λ ± (δ) and Λ ± (M ). These are provided in the result below, which we consider to be the technical heart of this section.
Finally we prove part (6). Let M ⊂ C satisfy (hM). The invertibility of L ± 0 + λ for λ ∈ M follows directly from (HP2). Since M is compact and the map is continuous, we immediately see that there exists C 3 > 0 such that and all λ ∈ M . We now establish some basic facts concerning the operator ∆ h . In particular, we extend the realvalued results from [1] to complex-valued functions. We emphasize that the inequalities in Lemma 3.6 in general do not hold for the imaginary parts of these inner products.

For all
Lemma 3.6. Assume that (Hα1) is satisfied and pick f ∈ H 1 C (R). Then the following properties hold.

For all
Similarly we have For λ ∈ C we may compute (3.36) Taking λ = 1 gives the third property. We now set out to prove Proposition 3.4. In Lemmas 3.7 and 3.8 we construct weakly converging sequences that realize the infima in (3.18)-(3.21). In Lemmas 3.9-3.11 we exploit the structure of our operators (3.3) and (3.4) to recover bounds on the derivatives of these sequences that are typically lost when taking weak limits. Recall the constant C 2 > 0 defined in Lemma 3.1, which does not depend on δ > 0.
In particular, we can choose h > 0 and M > 0 in such a way that |ũ h | ≤ M and |u 0 | ≤ M for In particular, we see that (φ, ψ) is a weak solution to (L we get φ ∈ L 2 and hence φ ∈ H 2 . Since we already know that ψ ∈ H 1 , we may apply Lemma 3.1 and (3.39) to obtain Recall the constant C 3 > 0 from Lemma 3.1, which only depends on the choice of the set M ⊂ C that satisfies (hM).
where the constant C 3 is defined in Lemma 3.1.
The same statements hold upon replacing L The proof of this result is almost identical to the one of Lemma 3.7 and will be omitted. In our arguments below, we often consider the sequences {(h j , φ j , ψ j )} and {(λ j , h j , φ j , ψ j )} defined in Lemmas 3.7 and 3.8 in a similar fashion. To streamline our notation, we simply write {(λ j , h j , φ j , ψ j )} for all these sequences, with the understanding that λ j = δ when referring to Lemma 3.7.
As argued in the proof of Lemma 3.7, it is possible to choose h > 0 in such a way that By taking a subsequence if necessary, we assume from now on that h j < h for all j. It remains to find a positive lower bound for (φ, ψ) L 2 .
Proof. We first consider the sequence for Λ + . Using L + hj ,λj (φ j , ψ j ) = (θ j , χ j ) and Re ∆ hj φ j , φ j = 0 = Re φ j , φ j = Re ψ j , ψ j , which follow from Lemma 3.6, we obtain We write λ max = δ 0 in the setting of Lemma 3.7 or λ max = max{|z| : z ∈ M } in the setting of Lemma 3.8. We write (3.54) Using the Cauchy-Schwarz inequality we now obtain (3.56) Squaring this equation and using the standard inequality 2µω ≤ µ 2 + ω 2 , this implies that (3.57) In particular, we see We now look at the sequence for Λ − . Using L − hj ,λj (φ j , ψ j ) = (θ j , χ j ) and Re ∆ hj φ j , φ j = 0 = Re φ j , φ j = Re ψ j , ψ j , which follow from Lemma 3.6, we obtain Using the Cauchy-Schwarz inequality we now obtain This is the same equation that we derived for Λ + . Hence we again obtain We would like to point out that in the following lemma we heavily exploit the bistable structure of the non-linearity g. Moreover, we are aided by the fact that the off-diagonal elements are constant, which allows us to keep the cross-terms under control. In fact, one might be tempted to think that it is sufficient to note that the eigenvalues of the matrix −g u (0) 1 −ρ γρ all have positive real part, as then one would be able to find a basis in which this matrix is positive definite. However, passing over to another basis destroys the structure of the diffusion terms and therefore does not give any insight. (3.64) Here we write λ min = 0 in the setting of Lemma 3.7 or λ min = min{Re λ : λ ∈ M } in the setting of Lemma 3.8, together with Proof. Again we first look at the sequence for Λ By the bistable nature of our non-linearity g, we can choose m to be a positive constant such that for Here r 0 is the constant appearing in the choice of our function g in (HS). Then we obtain, using Re φ j , φ j = Re ψ j , ψ j = 0 and Re −∆ hj φ j , φ j ≥ 0, which we know from Lemma 3.6, that (3.67) We assumed that 0 < ρ < 1 so we see that 1−ρ −ρ < 0. We set (3.68) Now we obtain (3.69) Note that the denominator 4( ρ 1−ρ 1 2 γρ + γρ + Re λ j ) is never zero since we can assume that λ * is small enough to have Re λ j ≥ −λ * > −γρ. Using the identity χ j = −ρφ j +c hj ψ j + γρψ j + λ j ψ j (3.70) and the fact that Re ψ j , ψ j = 0, we also have Hence we must have that (3.72) Combining this bound with (3.67) yields the estimate (3.73) We now look at the sequence for Λ − . Let m and a be as before. Then we obtain, using L (3.75) Arguing as in (3.69) with different constants, we obtain (3.76) Note that the denominator 4(a + Re λ j ) is never zero since we can assume that λ * is small enough to have Re λ j ≥ −λ * > −a. Using the identity and the fact that Re φ j , φ j = 0, we also have Hence we must have that Combining this with the estimate (3.74) and noting that 1−ρ we note that β + j ≤ β and β − j ≤ β for all j since ρ < 1 and since β + j and β − j are maximal for Re λ = λ min . Therefore, in both cases, we obtain (3.82) and thus, again using the inequality 2µω ≤ µ 2 + ω 2 for µ, ω ∈ R, Proof. Without loss of generality we assume that 1 2 min{a, 1 2 ργ} + λ min > 0. Write This allows us to compute Proof of Proposition 3.4. We first choose 0 < δ < δ 0 and consider the setting of Lemma 3.7. Sending j → ∞ in (3.84), Lemma 3.7 implies (3.90) Solving this quadratic inequality, we obtain The analogous computation in the setting of Lemma 3.8 yields For the remainder of this proof, we only consider the operators L + h,δ , noting that their counterparts L − h,δ can be treated in an identical fashion. Seeking a contradiction, let us assume that L + h,δ (H 1 ) = L 2 , which implies that there exists a non-zero (θ, χ) ∈ L 2 orthogonal to L + h,δ (H 1 ). For any φ ∈ C ∞ c (R), we hence obtain By definition this implies that θ has a weak derivative and thatc h θ = −∆ h θ − g u (ũ h )θ + δθ − ρχ ∈ L 2 (R). In particular, we see that θ ∈ H 1 (R).
For any ψ ∈ C ∞ c (R) a similar computation yields Again, this means that χ has a weak derivative and in factc h χ = θ + (γρ + δ)χ. In particular it follows that χ ∈ H 1 (R).
We therefore conclude that holds for all (φ, ψ) ∈ H 1 . Since H 1 is dense in L 2 this implies that L − h,δ (θ, χ) = 0. Since we already know that L − h,δ is injective, this means that (θ, χ) = 0, which gives a contradiction. Hence we must have L

The point and essential spectrum
In this section we discuss several properties of the operator that arises after linearising the travelling pulse MFDE (2.12) around our wave solution (u h , w h ). The main goals are to determine the Fredholm properties of this operator. In particular we show that both the linearised operator and its adjoint have Fredholm index 0 and that they both have a one-dimensional kernel. Moreover, we construct a family of kernel elements of the adjoint operator that converges to (φ − 0 , ψ − 0 ), the kernel element of the operator L − 0 . Pick 0 < h < min{h * , h}, where h * is given in Theorem 2.1 and h is characterized by (3.51). We recall the operator L h : H 1 → L 2 , introduced in §2, which acts as In addition, we write L * h : H 1 → L 2 for the formal adjoint of L h , which is given by We emphasize that L h and L * h correspond to the operators L for the family featuring in (hFam). Finally, we introduce the notation The results of this section should be seen as a bridge between the singular perturbation theory developed in §3 and the spectral analysis preformed in §5. Indeed, one might be tempted to think that most of the work required for the spectral analysis of the operator L h can already be found in Proposition 3.2 and Proposition 3.3. However, the problem is that we have no control over the δ-dependence of the interval (0, h 0 (δ)), which contains all values of h for which L h + δ = L + h,δ is invertible. In particular, for fixed h > 0 we cannot directly conclude that L + h,δ is invertible for all δ in a subset of the positive real axis.
Our main task in this section is therefore to remove the δ-dependence and study L h and L * h directly. The main conclusions are summarized in the results below.  and

Upon introducing the spaces
and the operator L h : X h → Y h is invertible and there exists a constant C unif > 0 such that for each 0 < h < h * * we have the uniform bound A direct consequence of these results is that the zero eigenvalue of L h is simple. In addition, these results allow us to construct a quasi-inverse for L h that we use heavily in §5 and §6.
, which together with (4.5) completes the proof.  such that for all Θ ∈ L 2 and each 0 < h < h * * the pair is the unique solution to the problem that satisfies the normalisation condition Proof. Fix 0 < h < h * * and Θ ∈ L 2 . Upon defining we see that Θ + γ h [Θ]Φ + h ∈ Y h . In particular, Proposition 4.2 implies that the problem has a unique solution Ψ ∈ X h , which we refer to as L qinv h Θ.
The results in [21,48] allow us to read off the Fredholm properties of L h from the behaviour of this operator in the limits ξ → ±∞. In particular, we let L h,∞ be the operator defined by (4.16) This system has constant coefficients. For λ ∈ C we introduce the notation We show that for λ in a suitable right half-plane the operator L h,∞;λ is hyperbolic in the sense of [21,48], i.e. we write and show that det(∆ L h,∞;λ (iy)) = 0 for all y ∈ R. In the terminology of [21,48], this means that L h + λ is asymptotically hyperbolic. This allows us to compute the Fredholm index of L h + λ.
Remark 4.5. From this section onward we assume that (Hα2) is satisfied. This is done for technical reasons, allowing us to apply the results from [21]. In particular, this condition implies that the function ∆ L h,∞;λ (z) defined in (4.18) is well-defined in a vertical strip |Re (z)| < ν .
Before we consider the Fredholm properties of L h + λ, we establish a technical estimate for the function ∆ L h,∞;λ , which we need in §6 later on.
Lemma 4.7. Assume that (HP1), (HS), (Hα1) and (Hα2) are satisfied. Fix 0 < h < min{h * , h} and S ⊂ C compact in such a way that Re λ > −λ for all λ ∈ S. Then there exist constants κ > 0 and Γ > 0, possibly depending on h and S, such that for all z = x + iy ∈ C with |x| ≤ κ and all λ ∈ S we have the bound Proof. Using assumption (Hα2) we can pick κ 1 > 0 and Γ 1 > 0 in such a way that the bound holds for all z = x + iy ∈ C with |x| ≤ κ 1 .

(4.25)
Since S is compact we can find Y > 0 such that for all z = x + iy ∈ C with |y| ≥ Y and |x| ≤ k 1 and all λ ∈ S we have Re det(∆ L h,∞;λ (z)) (4.26) Seeking a contradiction, let us assume that for each 0 < κ ≤ κ 1 and each Γ > 0 there exist λ ∈ S and z = x + iy ∈ C with |x| ≤ κ and |y| ≤ Y for which Then we can construct a sequence {κ n , z n , λ n } with 0 < κ n ≤ κ 1 for each n, κ n → 0, λ n ∈ S for each n and z n = x n + iy n ∈ C with |x n | ≤ κ n and |y n | ≤ Y in such a way that | det(∆ L h,∞;λn (z n ))| < 1 n for each n. By taking a subsequence if necessary we see that λ n → λ for some λ ∈ S and z n → iy for some y ∈ R with |y| ≤ Y . Since det(∆ L h,∞;λ (z)) is continuous as a function of λ and z it follows that det(∆ L h,∞;λ (iy)) = lim n→∞ det(∆ L h,∞;λn (z n )) = 0, (4.28) which contradicts Lemma 4.6. Hence we can find κ > 0 and Γ > 0 as desired.
Proof of Proposition 4.1. We have already seen in Lemma 4.6 that L h +λ is asymptotically hyperbolic in the sense of [21,48]. Now according to [21,Theorem 1.6], we obtain that L h + λ is a Fredholm operator and that the following identities hold is the Fredholm index of L h + λ.
We follow the proof of [48,Theorem B]. For 0 ≤ ϑ ≤ 1, we let the operator L ϑ (h) be defined by Note that the operator L ϑ (h) is asymptotically hyperbolic for each ϑ and thus [21,Theorem 1.6] implies that these operators L ϑ (h) are Fredholm. Moreover, the family L ϑ (h) varies continuously with ϑ in B(H 1 , L 2 ), which means the Fredholm index is constant. In particular, we see that where the last equality follows from [21, Theorem 1.7].
We can now concentrate on the kernel of L h . The goal is to exclude kernel elements other than Φ + h . In order to accomplish this, we construct a quasi-inverse for L h by mimicking the approach of [38, Proposition 3.2]. As a preparation, we obtain the following technical result.
and C 0 is defined in Proposition 3.2. The factor 4 in the definition is for technical reasons in a later proof. We assume from now on that 0 < h ≤ h * 1 . Using L h Φ + h = 0 we readily see Recall that Φ − 0 L 2 = 1. Since 1 < λ −1 , we may use Proposition 3.2 to obtain (4.37) Remembering Φ − 0 , Φ + h > 0 and using Cauchy-Schwarz, we see that (4.38) Hence we must have  defined for all 0 < h < h * * , such that for all Θ ∈ L 2 the pair is the unique solution to the problem that satisfies the normalisation condition In addition, there exists C > 0 such that for all 0 < h < h * * and all Θ ∈ L 2 we have the bound Proof. Fix 0 < λ < min{ 1 2 , δ 0 } and let 0 < h ≤ min{h * , h, h 0 (λ)} be given, where h 0 (λ) is defined in Proposition 3.2. We define the set Pick Θ ∈ L 2 . We look for a solution (γ, Ψ) ∈ R × Z 1 of the problem which is the unique value for γ for which Recall the constant C unif from (4.35). With Proposition 3.2 we obtain for some C 1 that is independent of h, λ. Here we used that λ < 1 and thus 1 + 1 λ < 2 λ . Exploiting λ < 1 2 and applying Lemma 4.8, we see that (4.50) Here we used that and (4.53) For Ψ 1 , Ψ 2 ∈ Z 1 , a second application of Proposition 3.2 yields (4.54) Applying Proposition 3.2 for the final time, we see (4.55) In view of these bounds, we pick λ to be small enough to have C 3 λ < 1 2 and C 5 λ < 1 2 . In addition, we write h * * = min{h * 1 , h 0 (λ)} and pick 0 < h < h * * . Then T : Z 1 → Z 1 is a contraction, so the fixed point theorem implies that there is a uniqueL qinv h (Θ) ∈ Z 1 for which (4.57)  where L * h is the formal adjoint of L h . Proof. By differentiating the differential equation (2.12) we see that L h Φ + h = 0. We know that (u h , w h ) − (u 0 , w 0 ) → 0 ∈ H 1 . Since (u 0 , w 0 ) decays exponentially, we get (u 0 , w 0 ) ∈ L 2 . Hence we can assume that h * * is small enough such that Φ + h ∈ L 2 for all 0 < h < h * * . Since L h Φ + h = 0 we obtain from the differential equation that also (Φ + h ) ∈ L 2 . In particular, we see that Φ + h ∈ H 1 and hence Φ + h ∈ ker(L h ). Lemma 4.11. Assume that (HP1), (HS), (Hα1) and (Hα2) are satisfied. Let 0 < h < h * * be given. Then we have where L * h is the formal adjoint of L h . Proof. We show that dim(ker(L h )) = 1. Since Φ + h ∈ ker(L h ), we assume that there exists Ψ ∈ ker(L h ) in such a way that Ψ is not a scalar multiple of Φ + h . Suppose first that Ψ, Φ − 0 = 0. Then Lemma 4.9 yieds by linearity ofL qinv which gives a contradiction. Hence we suppose that Ψ, Φ − 0 = 0. In the proof of Lemma 4.8 we saw that Φ + h , Φ − 0 = 0. As such, we can pick a, b ∈ R \ {0} in such a way that Again it follows from Lemma 4.9 that aΦ + ε + bΨ = 0 which gives a contradiction. Therefore such kernel element Ψ does not exist.
Since we already know that Φ + h ∈ ker(L h ), we must have dim ker(L h ) = 1, which completes the proof.
The remaining major goal of this section is to find a family of elements To establish this, we repeat part of the process above for the adjoint operator L * h . The key difference is that we must construct the family Φ − h by hand. This requires a significant refinement of the argument used above to characterize ker(L * h ). First we need a technical result, similar to Lemma 4.8. (4.65) Recall the constant C unif from (4.35). Proposition 3.2 yields (4.66) Using Lemma 3.5 and the fact that c h converges to c 0 and g u (u h ) to g u (u 0 ), it follows that as h ↓ 0. Possibly after decreasing h 0 (λ) > 0, we hence see that holds for all 0 < h < min{h * * , h 0 (λ)}. is the unique solution to the problem that satisfies the normalisation condition Furthermore, there exists C * > 0, such that for all 0 < h < h * * and all Θ ∈ L 2 we have the bound Proof. We define the set Pick Θ ∈ L 2 . We look for a solution (γ, Ψ) ∈ R × Z 1 of the problem which is the unique value for γ for which From now on the proof is identical to that of Lemma 4.9, so we omit it. which leads to a contradiction. Hence, we can define the kernel element Φ − h of L * h as follows: we obtain, upon defining Using Lemma 4.13, we can estimate Proof. Fix 0 < h < h * * . Clearly L h : X h → Y h is a bounded bijective linear map, so the Banach isomorphism theorem implies that L −1 h : Y h → X h is bounded. Now let δ > 0 be a small constant such that δC unif < 1. Without loss of generality we assume that 0 < h * * ≤ h 0 (δ) and that This is possible by Lemma 4.14. Pick any Ψ ∈ X h . Remembering that Ψ, Φ − h = 0 and L h Ψ, Φ − h = 0, we obtain the estimate (4.86) Applying Proposition 3.2, we hence see (4.87) We therefore get the bound

The resolvent set
In this section we prove Theorem 2.2 by explicitly determining the spectrum of the operator −L h defined in (2.18). Our approach hinges on the periodicity of this spectrum, which we describe in our first result.  From this section onward we need to assume that (HP2) is satisfied. Indeed, this allows us to lift the invertibility of L 0 + λ to L h + λ simultaneously for all λ in appropriate compact sets. For any λ ∈ C, Ψ = (φ, ψ) ∈ H 1 and x ∈ R we obtain since plh ∈ 2πiZ for all l > 0. In particular, we can compute Since e p and e −p are invertible operators on H 1 and L 2 respectively, we know that the spectrum of L h equals that of e −p L h e p and thus that of L h + pc h .
Since L h has a simple eigenvalue at zero, it is relatively straightforward to construct a small neighbourhood around the origin that contains no other part of the spectrum. Exploiting the results from §4, it is possible to control the size of this neighbourhood as h ↓ 0.
Proposition 5.2. Assume that (HP1),(HP2), (HS), (Hα1) and (Hα2) are satisfied. There exists a constant λ 0 > 0 such that for all 0 < h < h * * the operator L h + λ : Proof. Fix 0 < h < h * * and Θ ∈ L 2 . We recall the notation (γ h [Θ], L qinv h Θ) from Corollary 4.4 for the unique solution (γ, Ψ) of the equation Now for λ ∈ C with |λ| small enough, but λ = 0, we want to solve the equation withΨ ∈ X h , we see that In particular, we must find a solutionΨ ∈ X h for the equation which we can rewrite as We choose λ 0 in such a way that 0 < λ 0 C unif < 1. Then it is well-known that I − λL −1 h is invertible as an operator on X h for 0 < |λ| < λ 0 . Since λL −1 h L qinv h Θ ∈ X h , we see that (5.10) indeed has a unique solutionΨ ∈ X h . Hence the equation (L h −λ)Ψ = Θ always has a unique solution. Proposition 4.1 states that L h −λ is Fredholm with index 0, which now implies that L h −λ is indeed invertible.

Region R 2 .
We now show that in an appropriate right half-plane, which can be chosen independently of h, the spectrum of −L h is empty. The proof proceeds via a relatively direct estimate that is strongly inspired by [1, Lemma 5].

Region R 3 .
This region is the most delicate to handle on account of the periodicity of the spectrum. Indeed, one cannot simply take Im λ → ±∞ in a fashion that is uniform in h. We pursue a direct approach here, using the Fourier transform to isolate the problematic part of L h + λ, which has constant coefficients. The corresponding portion of the resolvent can be estimated in a controlled way by rescaling the imaginary part of λ. We remark that an alternative approach could be to factor out the periodicity in a more operator-theoretic setting, but we do not pursue such an argument here.
Pick λ ∈ C with λ 0 < | Im λ| ≤ |c h | h π and write Introducing the new variable τ = Im λξ, we can write the eigenvalue problem (L h + λ)(v, w) = 0 in the form

(5.17)
Our computations below show that the leading order terms in the appropriate |λ im | → ∞ limit are encoded by the 'homogeneous operator' H h,λ that acts as Writing H h,λ for the Fourier symbol associated to H h,λ , we see that α k cos(hkλ im ω) − 1 . can only be satisfied if the inequalities In particular, upon choosing ε = 1 4 , we see that and hence Since c h → c 0 = 0 as h ↓ 0, the desired inequalities (5.21) follow.
Lemma 5.5. Assume that (HP1),(HP2), (HS), (Hα1) and (Hα2) are satisfied. Then there exists a constant C > 0 such that for all ω ∈ R and 0 < h < h * * and all λ ∈ C with |λ| > λ 0 and | Im λ| ≤ |c h | h π, we have the inequality Proof. We show that H h,λ (iω) is bounded away from 0, uniformly in h, λ and ω. To do so, we show that the real part of H h,λ (iω) can be bounded away from zero, whenever the imaginary part is small, i.e. when (5.21) holds.
Proof. Since Proposition 4.1 implies that L h + λ is Fredholm with index zero, it suffices to prove that L h + λ is injective.
Writev andŵ for the Fourier transforms of v and w respectively. For f ∈ L 2 with Fourier transformf , the identity In particular, we obtainv which using Lemma 5.5 implies that for some constant C > 0 that is independent of h, λ and ω.
Since Ψ is an eigenfunction, (5.17) hence yields Furthermore applying a Fourier Transform to the second line of (5.17), we find Our choice λ 3 ≤ 1 2 ργ implies that −ργ + λ r is bounded away from 0. We may hence writê which yields the bound for some constant C > 0. Therefore we obtain that for some constant C , which is independent of λ, h and v. Clearly this is impossible for v = 0 if Furthermore if v = 0, then clearly also w = 0. Therefore Ψ = 0, allowing us to conclude that L h + λ is invertible.
We conclude our spectral analysis by considering the remaining region R 4 . This region is compact and bounded away from the origin, allowing us to directly apply the theory developed in §3.

Green's functions
In order to establish the nonlinear stability of the pulse solution (u h , w h ), we need to consider two types of Green's functions. In particular, we first study G λ (ξ, ξ 0 ), which can roughly be seen as a solution of the equation where δ is the Dirac delta-distribution. We then use these functions to build a Green's function G for the linearisation of the LDE (2.1) around the travelling pulse solution. An important difficulty in comparison to the PDE setting is caused by the discreteness of the spatial variable j. In particular, we cannot use a frame of reference in which the solution (u h , w h ) is constant without changing the structure of the equation (2.1). The Green's function G hence will be the solution to a non-autonomous problem that satisfies a shift-periodicity condition. Nevertheless, one can follow the technique in [5] to express G in terms of a contour integral involving the functions G λ . A significant part of our effort here is concerned with the construction of these latter functions. Indeed, previous approaches in [3,36] all used exponential dichotomies or variation-of-constants formula's for MFDEs with finite-range interactions. These tools are no longer available for use in the present infinite-range setting. In particular, we construct the functions G λ in a direct fashion using only Fredholm properties of the operators L h + λ. This makes it somewhat involved to recover the desired exponential decay rates and to properly isolate the meromorphic terms of order O(λ −1 ).
From now on, we will no longer explicitly use the h-dependence of our system. To simplify our notation, we fix 0 < h < h * * and write We emphasize that from now on all our constants may (and will) depend on h.
We will loosely follow §2 of [36], borrowing a number of results from [5,34] at appropriate times. In particular, we start by considering the linearisation of the original LDE (2.1) around the travelling pulse solution U (t) given by (2.21). To this end, we introduce the Hilbert space in which Mat 2 (R) is the space of 2 × 2-matrices with real coefficients which we equip with the maximum-norm | · |. For any V ∈ L 2 , we often write V = V (1,1) V (1,2) V (2,1) V (2,2) , when we need to refer to the component sequences V (i,j) ∈ 2 (Z; R). For any t ∈ R we now introduce the linear operator A(t) : L 2 → L 2 that acts as for v ∈ 2 (Z; R) and w ∈ 2 (Z; R). With all this notation in hand, we can write the desired linearisation as the ODE d dt V(t) = A(t) · V(t) (6.6) posed on L 2 . Fix t 0 ∈ R and j 0 ∈ Z. Consider the function that is uniquely determined by the initial value problem Here we have introduced where I ∈ Mat 2 (R) is the identity matrix. We remark that G j0 j (t, t 0 ) is an element of Mat 2 (R) for each j ∈ Z.
This function G is called the Green's function for the linearisation around our travelling pulse. Indeed, the general solution of the inhomogeneous equation where now V (t) ∈ 2 (Z; R 2 ) ∼ = 2 (Z; R 2×1 ) and F (t) ∈ 2 (Z; R 2 ) ∼ = 2 (Z; R 2×1 ), is given by Introducing the standard convolution operator * , this can be written in the abbreviated form The main result of this section is the following proposition, which shows that we can decompose the Green's function G into a part that decays exponentially and a neutral part associated with translation along the family of travelling pulses. Proposition 6.1. Assume that (HP1),(HP2), (HS), (Hα1) and (Hα2) are satisfied. For any pair t ≥ t 0 and any j, j 0 ∈ Z, we have the representation whileG satisfies the found for some K > 0 andδ > 0. The constant Ω > 0 is given by Furthermore for any t ≥ t 0 we have the representation which can be abbreviated as Our first task is to define G λ in a more rigorous fashion. Recall the operator L ∞;λ and the function ∆ L ∞;λ from Lemma 4.6. We will show that L ∞;λ has a Green's function which takes values in the space Mat 2 (R) and that this function has some nice properties. Recall the constantλ from Lemma 4.6. For each λ ∈ C with Re λ ≥ −λ 2 we define G ∞;λ : R → Mat 2 (R) by We also introduce the notation G ∞ = G ∞;0 . (6.20) Here (Hα2) is essential to ensure that that these Green's functions decay exponentially.
Pick λ ∈ C \ σ(−L) with Re λ ≥ −λ 2 . Observe that We know that G ∞;λ (· − ξ 0 ) ∈ L 2 (R, Mat 2 (R)) since it decays exponentially. Therefore also Hence it is possible to define the function G λ by The next result shows that G λ can be interpreted as the Green's function of L + λ. It is based on [36, Lemma 2.6]. with Re λ ≥ −λ 2 we have that G λ (·, y) is continuous on R\{y} and C 1 -smooth on R\{y+kh : k ∈ Z}. Furthermore it satisfies for all ξ ∈ R and all f ∈ H 1 .
The link between our two types of Green's functions is provided by the following key result. It is based on [5,Theorem 4.2], where it was used to study one-sided spatial discretisation schemes for systems with conservation laws. Proposition 6.4. Assume that (HP1),(HP2), (HS), (Hα1) and (Hα2) are satisfied. Let χ > λ unif be given, where λ unif is as in Lemma 6.7. For all t ≥ t 0 the Green's function G j0 j (t, t 0 ) of (6.8) is given by where G λ is the Green's function of λ + L as defined in (6.25).
In the first half of this section we establish several basic facts concerning G λ that can be used to verify the integral representation (6.27). We then establish a meromorphic expansion for G λ that allows us to shift the integration path in (6.27) to the left of the imaginary axis. The decomposition (6.13) for the Green's function G can subsequently be read off from this expression. Lemma 6.5. Assume that (Hα1) and (Hα2) are satisfied. Consider any bounded function f : R → R which is continuous everywhere except at some ξ 0 ∈ R. Then ∆ h f is continuous everywhere except at {ξ 0 + hk : k ∈ Z}. Moreover, if f is differentiable except at ξ 0 and f is bounded, then ∆ h f is differentiable everywhere except at {ξ 0 + hk : Proof. For convenience we set ξ 0 = 0. Pick ξ ∈ R with ξ / ∈ {kh : k ∈ Z}. Then f is continuous in each point ξ + kh for k ∈ Z. Choose ε > 0. Since f is bounded and ∞ j=1 |α j | < ∞, we can pick K > 0 in such a way that For j ∈ {1, ..., K − 1} we can pick δ j > 0 in such a way that for all y ∈ R with |y| < δ j . Let δ = min{δ j : 1 ≤ j < K} > 0. Then for y ∈ R with |y| < δ we obtain for n ∈ Z >0 , we can compute This allows us to estimate In particular the sequence {f n } converges uniformly to ∆ h f from which it follows that Lemma 4.7 implies that we can choose β * > 0 and K * > 0 in such a way that ∆ L ∞;λ (z) −1 ≤ K * 1+| Im z| (6.36) for all λ ∈ R and all z ∈ C with | Re z| ≤ 2β * . In particular, it follows that (y → ∆ L ∞;λ (iy) −1 ) ∈ L 2 (R). By the Plancherel Theorem it follows that G ∞;λ is a well-defined function in L 2 (R). In particular, it is bounded. Shifting the integration path in (6.19) in the standard fashion described in [39,48], we obtain the bound for all ξ, ξ 0 ∈ R and λ ∈ R.
We loosely follow the approach of [34, §5.1], which considers a similar setting for Green's functions for Banach space valued operators with finite range interactions. Pick λ ∈ R. We rewrite the definition of ∆ L ∞;λ given in (4.18) in the more general form For α ∈ R close to 0 we introduce the expression R L ∞;λ ;α by for z ∈ C unequal to α and | Re z| ≤ 2β * . Since we can compute we obtain the estimate for all y ∈ R, possibly after increasing K * . Exploiting the decomposition (6.39), we write where we have introduced and we can easily compute (6.48) Since R L ∞;λ ;α ∈ L 1 (R) we see that R α is continuous. Therefore G ∞;λ is continuous outside of ξ = 0. Similarly to [34,Eq. (5.79)] we observe that (6.50) using [34,Lemma 5.8]. In particular, we see that L ∞;λ G ∞;λ (ξ) = 0 (6.51) for all ξ outside of {hk : k ∈ Z}. Lemma 6.5 subsequently shows that G ∞;λ is C 1 -smooth outside of {hk : k ∈ Z}. Fix f ∈ H 1 . For any δ > 0 we may compute together with for all ξ ∈ R.
Proof. On account of Lemma 6.2 we can lift the results from [48,] to our current infinite range setting. The proof of these results are identical, since the estimate [48,Eq. (5.4)] still holds in our setting on account of (Hα2). A more detailed description for this procedure can be found in [6,.
Proof of Lemma 6.3. Pick λ ∈ C \ σ(−L) and compute The last statement follows immediately from this identity.

Set
Since u ∈ H 1 and hence continuous, we must have that u is continuous. As argued before ∆ h H (1,1) and ∆ h H (1,2) are also continuous. Hence we see that c d dξ H is continuous on R \ {ξ 0 } and thus that d dξ H is continuous on R \ {ξ 0 }. Therefore we obtain that G λ (·, ξ 0 ) is C 1 -smooth on R \ {ξ 0 + kh : k ∈ Z}. Now we show that for λ with sufficiently large real part, G λ (·, ξ 0 ) is bounded uniformly by a constant. This result is based on [5,Lemma 4.1].
We introduce G 0 λ as the Green's function of (λ + c d dξ ) viewed as a map from H 1 to L 2 . Luckily, it is well-known that this Green's function admits the estimate We can look for the Green's function G λ as the solution of the fixed point problem Since λ + L is invertible by Theorem 2.2, G λ must necessarily satisfy the fixed point problem (6.64).
Proof of Proposition 6.4. Fix j 0 ∈ Z and t 0 ∈ R. Since (6.8) is merely a linear ODE in the Banach space L 2 , it follows from the Cauchy-Lipschitz theorem that (6.8) indeed has a unique solution V : [t 0 , ∞) → L 2 . For any Z ∈ C ∞ c (R; L 2 ), an integration by parts yields (6.72) We want to show that the function V j (t) := G j0 j (t, t 0 ) defined by (6.27) coincides with V on [t 0 , ∞). To accomplish this, we define and show that V is a weak solution to (6.8) in the sense that holds for all Z ∈ C ∞ c (R; L 2 ). Indeed, the uniqueness of weak solutions then implies that V = V. Note first that V (t) = 0 for t < t 0 , which can be seen by using (6.61) and taking χ → ∞ in (6.27). We write y = hj 0 + ct 0 , χ − = χ − iπc h and χ + = χ + iπc h . We see that since V (t) = 0 for t < t 0 . Moreover we write G j (t) = G λ (hj + ct, y). (6.76) Using our definition of V (t), we have where The permutation of the summations and integrations is allowed by Legesgue's theorem because Z and dZ dt are compactly supported and G λ is uniformly bounded by (6.61). Fix χ − ≤ λ ≤ χ + and j ∈ Z. Using the change of variable x = hj + ct we obtain Exploiting the fact that Z j and therefore Z j is compactly supported, (6.26) yields Now since Z j is compactly supported, we can exchange sums and integrals in equation (6.77). This allows us to compute as desired.
Proof. Lemma 6.8 implies that |α k |(2e δhk + 2)), (6.89) where the last sum converges by (Hα2), possibly after decreasing δ > 0. Using the fact that we hence see that there exists a constant K 3 > 0 such that The proof for the bound on (Φ − ) is identical. We recall the spaces together with the operators L −1 in the spaces B(X, X) and in B(Y, X) that were defined in Proposition 4.2. We also recall the notation L qinv Θ that was introduced in Corollary 4.4 for the unique solution Ψ of the equation in the space X, which is given explicitly by We now exploit these operators to decompose the Green's function of λ + L into a meromorphic and an analytic part. This result is based on [36, Lemma 2.7].
Combining this estimate with Lemma 6.8 and Lemma 6.12 yields the desired bound for some constants K 13 > 0 and ω > 0.
We write where λ is defined in the proof of Lemma 6.14.
We now use the spectral properties of −L to also decompose the Green's function G(t, t 0 ) into two parts. We obtain the following result, which is based on [36, Corollary 2.8].
Proof of Theorem 2.3. The Green's function bounds from Proposition 6.1 are strong enough to follow the program outlined in the proof of [36, Prop 2.1]. In particular, using standard fixed point arguments one can foliate the p -neighbourhood of the wave U with the one-parameter family formed by the stable manifolds of the translates U (· + ϑ). The full details can be found in [55, §8].