Travelling wave solutions for fully discrete FitzHugh-Nagumo type equations with infinite-range interactions

We investigate the impact of spatial-temporal discretisation schemes on the dynamics of a class of reaction-diffusion equations that includes the FitzHugh-Nagumo system. For the temporal discretisation we consider the family of six backward differential formula (BDF) methods, which includes the well-known backward-Euler scheme. The spatial discretisations can feature infinite-range interactions, allowing us to consider neural field models. We construct travelling wave solutions to these fully discrete systems in the small time-step limit by viewing them as singular perturbations of the corresponding spatially discrete system. In particular, we refine the previous approach by Hupkes and Van Vleck for scalar fully discretised systems, which is based on a spectral convergence technique that was developed by Bates, Chen and Chmaj.


Introduction
In this paper, we consider spatial-temporal discretisations of a class of reaction-diffusion systems that contains the FitzHugh-Nagumo partial differential equation (PDE). This PDE is given by (1.1) Here g is the bistable, cubic nonlinearity g(u; r) = u(1 − u)(u − r) with r ∈ (0, 1), while ρ > 0 and γ > 0 are positive constants. In particular, our goal is to show that travelling waves for the system (1.1) persist under these spatial-temporal discretisations. As such, we contribute to the broad study of numerical schemes and their impact on the solutions under consideration, which has produced an immense quantity of literature. The main distinguishing feature is that we are interested in structures that persist for all time, while almost all of the studies in this area focus on finite time estimates.

Pulse propagation
The system (1.1) was introduced in the 1960s [27,29] as a simplification of the Hodgkin-Huxley equations, which were used to describe the propagation of spike signals through the nerve fibers of giant squids [35]. After observing similar pulse solutions for the system (1.1) numerically [28], a more rigorous, analytical approach to understanding these pulse solutions turned out to be rather delicate. Indeed, many new tools have been developed, some even very recently, to construct these pulses and analyse their stability in various settings. These techniques include geometric singular perturbation theory [7,34,44,45], the variational principle [11], Lin's method [8,9,48], and the Maslov index [15,16]. Pulse solutions for the system (1.1) take the form (u, w)(x, t) = (u 0 , w 0 )(x + c 0 t) (1.2) for some wavespeed c 0 and smooth wave profiles u 0 , w 0 that satisfy the limits Spatially discrete systems It is well-known that electrical pulses can only move through nerve fibres at appropriate speeds if the nerves are insulated with a myelin coating. This coating admits regularly spaced gaps at the so-called nodes of Ranvier [56]. In fact, through a process called saltatory conduction, excitations of these nerves appear to jump from one node to the next [50]. Since the FitzHugh-Nagumo PDE (1.1) does not take this discrete structure into account directly, it has been proposed [47] to, instead, model these phenenomena using a so-called lattice differential equation (LDE). For example, by applying a nearest-neighbour spatial discretisation to (1.1), we arrive aṫ (1.4) where the variable j ranges over the lattice Z. In the system (1.4), the variable u j represents the potential at the jth node of the nerve fibre, while the variable w j describes a recovery component. Finally, we have τ ∼ h −2 , where h > 0 is the distance between subsequent nodes. We emphasize that the time variable remains continuous.
Spatialy discrete travelling pulses for the system (1.4) take the form (u, w) j (t) = (u 0 , w 0 )(j + c 0 t), (1.5) for some wavespeed c 0 , again with the limits (1.3). Plugging the Ansatz (1.5) into the LDE (1.4) yields the functional differential equation of mixed type (MFDE) c 0 u 0 (ξ) = τ [u 0 (ξ + 1) + u 0 (ξ − 1) − 2u 0 (ξ)] + g(u 0 (ξ); r) − w 0 (ξ) (1.6) in which ξ = j + c 0 t. In [37,38], Hupkes and Sandstede developed an infinite dimensional version of the exchange lemma to show that the system (1.4) admits nonlinearly stable travelling pulse solutions. They relied heavily on the existence of exponential dichotomies for MFDEs, which were established in [33,54]. In addition, we established the existence and nonlinear stability of pulse solutions for a spatially periodic version of (1.4) [59] by building on a spectral convergence method developed by Bates, Chen and Chmaj [2]. The spectral convergence method plays an important role in this paper as well and will be treated in more detail later on.
Infinite-range interactions Neural field models aim to describe the dynamic behaviour of large networks of neurons. In neural networks, neurons interact with each other over large distances through their interconnecting nerve axons [3,5,6,55]. It has been proposed [5,Eq. (3.31)] to capture these long distance interactions using an infinite-range version of the system (1.4). To be concrete, we focus our discussion on the prototype systeṁ (1.7) This system can also be obtained directly from the PDE (1.1) by using an infinite-range spatial discretisation.
We emphasize that infinite-range interactions also arise naturally when considering discretisations of fractional Laplacians [13]. Indeed, such operators are intrinsically nonlocal and are used in many physical systems that feature nonstandard diffusion processes, such as amorphous semiconductors [30] and liquid crystals [14].
Spatial-temporal discretisations Our main goal here is to understand the impact of temporal discretisation schemes on the behaviour of travelling wave solutions of the system (1.7). This is a relatively novel area of study, although a handful of results have been established for scalar problems. For example, Bambusi, Faou, Greébert and Jézéquel constructed solutions to fully discrete Schrödinger equations with Dirichlet or periodic spatial boundary conditions in [1,24]. Most other studies have focused on spatial-temporal discretisations of the Nagumo PDE u t = u xx + g(u; r), (1.9) or, equivalently, temporal discretisations of the Nagumo LDĖ The PDE (1.9) and the LDE (1.10) can be seen as scalar versions of the FitzHugh-Nagumo PDE (1.1) and LDE (1.4) respectively.
The early works by Elmer and Van Vleck [17][18][19] provided ad-hoc techniques to understand the impact of spatial-, temporal-and spatial-temporal discretisations of the PDE (1.9) on the dynamics of travelling waves. In addition, Chow, Mallet-Paret and Shen [12] established the existence of travelling wave solutions to temporal discretisations of the LDE (1.10) by considering Poincare return maps for the dynamics of this LDE. These results were later expanded by Hupkes and Van Vleck [40], whose methods allowed them to address issues of uniqueness and parameter-dependence.
Let us also mention the recent series of papers [41][42][43] by Hupkes and Van Vleck, who study spatial discretisation schemes with an adaptive grid. That is, the authors consider a time dependent moving mesh method which aims to equidistribute the arclength of the solution under consideration.
In order to introduce the temporal discretisation schemes that we study in this paper, we briefly discuss the test problemv = λv (1.11) with λ < 0. Applying the forward-Euler discretisation scheme with time-step ∆t > 0 yields where n ∈ Z. Since a nontrivial solution of the test problem (1.11) converges to zero as t → ∞, the convergence v n → 0 should also be enforced. However, this yields the restriction 0 < ∆t < 2|λ| −1 , which cannot be satisfied for all λ < 0 for a fixed time-step ∆t > 0. In contrast, these issues do not occur for the backward-Euler discretisation scheme. For the test problem (1.11), this scheme yields (1.14) In particular, we see that v n → 0 for any value of λ < 0 and time-step ∆t > 0. A numerical scheme is called A(α) stable if this property holds for all λ in the wedge {z ∈ C \ {0} : Arg(−z) < α}. We note that the backward-Euler discretisation is A( π 2 ) stable.
In fact, the backward-Euler discretisation scheme is one of six so-called backwards differentiation formula (BDF) methods. These BDF methods are all A(α) stable for various coefficients 0 < α ≤ π 2 and have several convenient analytical properties. For this reason, we have to chosen to focus on these temporal discretisation schemes in this paper. We do, however, emphasize that there are other stable discretisation schemes which we could have used, see for example [32].
Applied to the Nagumo system, the backward-Euler discretisation scheme yields the evolution A travelling wave solution for the system (1.15) with wavespeed c takes the form with the limits lim As such, the travelling waves need to satisfy the system Hupkes and Van Vleck showed [40] that, for sufficiently large, rational values of M = (c∆t) −1 , the system (1.15) admits travelling wave solutions with wavespeed c. These travelling waves are constructed as perturbations of travelling wave solutions of the LDE (1.10). The corresponding transition from the semi-discrete setting to the fully discrete setting is highly singular, since a derivative is replaced by a difference. The rationality of M plays a key role here, as it ensures that the domain of the variable ξ in the system (1.18) is a discrete subset of the real line. This restriction arises naturally in the analysis, since it ensures we can use finitely many interpolations to go from a fully discrete to a spatially discrete setting.
Spectral convergence In order to analyse this singular perturbation, Hupkes and Van Vleck relied heavily on the previously mentioned spectral convergence method, which also plays an important role in [41-43, 58, 59]. This method was introduced in [2] to construct travelling wave solutions to an infinite-range version of the Nagumo LDE (1.10) in the continuum limit, i.e. when the discretisation distance h ∼ τ − 1 2 is sufficiently small. A key role in [2] is reserved for the family of operators (1.20) In particular, the authors in [2] fixed a constant δ > 0 and used the invertibility of the operator L 0 + δ to establish the invertibility of the operator L h + δ for h > 0 sufficiently small. Indeed, they considered weakly converging sequences {v n } and {w n } with L h v n + δv n = w n and tried to find a uniform (in h and δ) lower bound on the norm of v n in terms of the norm of w n . Such a lower bound prevents the limitless transfer of energy into oscillatory modes, a common concern when dealing with weakly converging sequences. The bistable nature of the nonlinearity g was used to control the behaviour at ±∞, while the local L 2 -norm can be bounded on the remaining compact set. We emphasize that this method requires a detailed understanding of the limiting operator L 0 .
In [40], this method was lifted to the fully discrete Nagumo equation (1.18). Writing M = p q with gcd(p, q) = 1, the corresponding limiting operator resembles a q times coupled version of the operator L h given by (1.19). For q = 2, this limiting operator takes the form 21) where u is the travelling wave solution of the LDE (1.10) with wavespeed c. Here the domain of the variables ζ and ξ is given by ζ ∈ {0, 1 2 } and ξ ∈ R, with the convention that v(ζ + 1, ξ) = v(ζ, ξ). Since the MFDE corresponding to (1.21) admits a comparison principle, the Fredholm properties of the operator K q follow directly from the general results in [39]. Hupkes and Van Vleck generalized the spectral convergence method to lift the Fredholm properties of the operator K q to the operator (1.22) in the regime M 1, again with ζ ∈ {0, 1 2 } and ξ ∈ 1 2 M −1 Z. The operator K M arises as the linearisation of the fully discrete system (1.18) around the travelling wave u, using the additional ζ variable to ensure that all ξ-shifted arguments are multiples of M −1 .
Results In this paper, we consider reaction-diffusion LDEs such as (1.7) and replace the temporal derivative by one of the six BDF discretisation schems. For example, applying the backward-Euler method to (1.7), we arrive at the prototype system (1.23) Our main result states that systems such as (1.23) admit travelling wave solutions. To achieve this, we extend the spectral convergence method that was developed in [40] for scalar LDEs with finite-range spatial interactions to the current setting, which features multi-component systems with infinite-range interactions. This generalisation is far from trivial and requires several technical obstructions to be resolved.
The first main obstacle is that the spectral convergence method hinges on the understanding of the corresponding limiting operator. Indeed, the analog of the operator K q from (1.21) for our system (1.23) does not admit a comparison principle, since this is not available for FitzHugh-Nagumo type systems. As such, very limited a-priori knowledge is available for this limiting operator, which forces us to prove many of its properties from scratch. For this, we mainly employ techniques from harmonic analysis.
The second main obstacle is that the system setting introduces several cross-terms that need to be controlled. Several key techniques from our earlier works [58,59] concerning spatially discrete systems can be adjusted to handle these cross-terms in the present fully-discrete setting. However, several crucial points in the analysis still require these terms to be handled with special care.
The remaining obstacles are directly related to the infinite-range interactions, which introduce several convergence issues that need to be overcome. It also requires us to establish more refined estimates on the decay rates of solutions to our limiting MFDE. We achieve this by employing an explicit representation of the corresponding inverse linear operator that was first introduced in [58].
Loss of uniqueness In [40], Hupkes and Van Vleck extensively studied the uniqueness and parameter-dependence of the travelling wave solutions of (1.15). The key observation is that the rationality of the variable M = (c∆t) −1 breaks the translational symmetry in the travelling wave problem, potentially allowing a family of solutions to exist. For example, one can apply an irrational phase shift to the continuous waveprofiles for (1.10) that underlies the perturbation argument discussed above. In this fashion, one could construct a different fully discrete wave for the same detuning parameter value r in the nonlinearity g(·; r). However, this is a very delicate issue. In particular, M = (c∆t) −1 is fixed in the analysis, so additional work is required to obtain results for fixed time-steps ∆t > 0.
For the backward-Euler discretisation scheme, this nonuniqueness can be made fully rigorous. In particular, Hupkes and Van Vleck showed that, for a fixed time step ∆t > 0 both the r(c) relation and the c(r) relation can be multi-valued. In particular, for a fixed value of c there can be multiple values of r for which a solution to the system (1.15) exists and vice-versa. This can be achieved by embedding the system (1.18) into an MFDE that admits a comparison principle, allowing it to be analysed using the techniques developed by Keener [46] and Mallet-Paret [52].
By contrast, the c(r) relation for travelling wave solutions to the PDE (1.9) and the LDE (1.10) are both single-valued. The same holds for the r(c) relation, with the single exception that it can be multi-valued for (1.10) in the special case c = 0 [21,36]. This reflects the well-known wave-pinning phenomenon caused by the broken translational symmetry of the lattice [4,20,23,36,46,53].
In this paper we study the r(c) and the c(r) relation for a fully-discrete version of the FitzHugh-Nagumo system. For the corresponding PDE (1.1) and LDE (1.4), numerical evidence [10,49] suggests that both these relations are at most 2-valued. In addition, theoretical results [8] for this PDE usually yield a locally unique r(c) relation. For the system (1.23) a comparison principle is not available, rendering a direct analysis similar to the one in [40] infeasible. Instead, we run several numerical simulations to investigate these issues. These computations indicate that both the r(c) and the c(r) relation are typically multi-valued. Indeed, the points (r, c) points at which we were able to find solutions appear to map onto a surface instead of a curve. That is, there exists an entire spectrum of travelling wave solutions with different wavespeeds to the same fully discrete system.
Acknowledgements. Both authors acknowledge support from the Netherlands Organization for Scientific Research (NWO) (grant 639.032.612).

Main result
Our main goal is to study the impact of several important temporal discretisation schemes on travelling wave solutions of reaction-diffusion LDEs of the forṁ This LDE is posed on the one-dimensional lattice j ∈ Z, but may have multiple components in the sense that U j ∈ R d for some integer d ≥ 1. We start by discussing the structural conditions that we impose on the LDE (2.1) and its travelling wave solutions in §2.1 respectively §2.2. In §2.3 we introduce the appropriate temporal discretisation schemes and formulate our main result. Finally, we discuss some numerical results concerning the nonuniqueness of the fully discrete travelling waves in §2.4.

The spatially discrete system
Besides a handful of exceptions [2,25,26,31,57,58], almost all results concerning LDEs of the form (2.1) assume that only finitely many of the coefficients α m in (2.1) are nonzero. However, following [2,58], we will impose the following much weaker conditions.
Assumption (HS1). The coefficients {α m } m∈Z>0 are diagonal d × d matrices and τ > 0 is a positive constant. There for some constant ν > 0, as well as the identity holds for all z ∈ (0, 2π) and all 1 ≤ i ≤ d diff .
In particular, the diffusion matrices {α m } m∈Z>0 only act directly on the first d diff components of U j . For example, for the FitzHugh-Nagumo LDĖ we have d = 2 and d diff = 1, while for the Nagumo LDĖ In particular, (HS1) ensures that (2.5) can be interpreted as the spatial discretisation of the FitzHugh-Nagumo PDE (1.1) on a grid with distance h, where τ = 1 h 2 . Additional remarks concerning this assumption in the scalar case d = 1 can be found in [2, §1].
We now turn to the spatially homogeneous equilibrium solutions to (2.1), which are roots of the nonlinearity G. We will assume that there are two r-independent equilibria P ± , but emphasize that they are allowed to be identical.
The temporal stability of these two equilibria P ± plays an essential and delicate role in our analysis. Indeed, it does not suffice to simply require that the eigenvalues of DG(P ± ) have strictly negative real parts, see the proof of [59,Lem. 4.6] for details. Following [59], we consider two auxiliary assumptions on the triplet (G, P − , P + ) to address this issue. Recalling the constant 1 ≤ d diff ≤ d from (HS1), we first write DG(U ; r) in the block form for any U ∈ R d and r ∈ (0, 1), taking DG [1,1] Assumption (HS3 r ). The triplet (G, P − , P + ) satisfies at least one of the following conditions.

Spatially discrete travelling waves
Our final two assumptions for (2.1) concern the existence and stability of travelling wave solutions that connect the equilibria P − and P + . These solutions take the form for some smooth profile U 0 and nonzero wavespeed c 0 . Substituting the Ansatz (2.10) into (2.1) and writing ξ = j + c 0 t, we see that the pair (c 0 , U 0 ) must satisfy the travelling wave MFDE together with the boundary conditions Assumption (HW1 r ). There exists a waveprofile U 0 and a wavespeed c 0 = 0 that solve the travelling wave MFDE (2.11) for r = r, together with the boundary conditions (2.12).
We now turn to the spectral stability of these travelling wave solutions. To this end, we introduce the operator L 0 : for the linearisation of (2.11) around the travelling wave U 0 , which acts as Here the operator ∆ 0 : where In addition, we introduce the formal adjoint L * 0 : We remark that the spectrum of L 0 is 2πic 0 -periodic on account of the identity see [58,Lem. 5.1]. We impose the following condition on the spectral properties of this operator L 0 .
Assumption (HW2 r ). There exist functions Φ ± 0 ∈ H 1 (R; R d ), together with a constantλ > 0 so that the following properties hold for the LDE (2.1) with r = r.
together with the normalisation ii) The spectrum of the operator −L 0 in the half-plane {z ∈ C : Re z ≥ −λ} consists precisely of the points 2πimc 0 with m ∈ Z, which are all eigenvalues of L 0 . Moreover, we have the identities and 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 nontrivial w ∈ Ker(L * − λ). In particular, the normalisation (2.19) implies that the eigenvalues 2πic 0 Z are all simple eigenvalues of −L 0 .
For the FitzHugh-Nagumo system (2.5), the assumptions (HW1 r ) and (HW2 r ) are both satisfied for all sufficiently small discretisation distances h > 0 and sufficiently small ρ > 0, see [

The fully discrete system
We aim to approximate solutions to (2.1) at discrete time intervals t = n∆t by We need to apply an appropriate discretisation scheme to the temporal derivative in (2.1). Although there are many different approximation schemes available, we mainly focus on the six so-called BDF methods. These methods are based on interpolation polynomials of different degrees. In particular, the BDF method of order k ∈ {1, 2, ..., 6} approximates U in (2.1) at t = n∆t by first constructing an interpolating polynomial of degree k through the k + 1 points {W ((n − n )∆t)} k n =0 and then computing the derivative of this polynomial at W (n∆t). As such, the temporal discretisations of the LDE (2.1) under consideration are of the form (2.23) The coefficients β k and {µ n;k } in (2.23) are given implicitly by the identities which must hold for any scalar function v. Here we have introduced the notation  Table 1: The coefficients µ n;k and β k associated to the BDF discretisation schemes as given by (2.24).
For convenience, the values of the coefficients β k and µ n;k can be found in Table 1. We note that the BDF method of order 1 is the well-known backward-Euler method.
Our main goal is to study travelling wave solutions to the fully discrete system (2.23), utilizing our assumptions for the spatially discrete system (2.1). Such solutions are given by the Ansatz for some wave speed c and profile Φ with the boundary conditions in a sense that we make precise below.
For notational convenience, we introduce the quantity M = (c∆t) −1 . Substituting the Ansatz (2.27) into (2.23) yields the system for all ξ that can be written as ξ = n + jM −1 for (j, n) ∈ Z 2 . Here we have introduced the discrete derivatives for all integers 1 ≤ l ≤ k and all Φ ∈ C l+1 (R; R d ), in which the constant C l ≥ 1 is independent of k, Φ and M . Indeed, this estimate shows that the regular derivative can be approximated by the discrete derivatives as the time step ∆t shrinks to zero. We emphasize that BDF discretisation schemes of order k ≥ 2 do not allow for a comparison principle, even when the original LDE does allow for one. This is a consequence of the existence of coefficients µ n;k > 0 that have n < k.
Most of our results, including our main theorem, require a restriction on the values of M that are allowed. In particular, upon fixing an integer q ≥ 1, we introduce the set M q = { p q : p ∈ N has gcd(p, q) = 1 and p ≥ q}. (2.32) Often, we introduce M = p q ∈ M q , which implicitly defines the integer p = p(M ) = qM . Moreover, we see that the natural domain for the values of ξ in the system (2.29), as well as in the boundary conditions (2.28), is precisely the set p −1 Z. Theorem 2.1. Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Fix a pair of integers 1 ≤ k ≤ 6 and q ≥ 1. Then there exist constants M * 1 and δ r > 0 so that for any that satisfy the following properties.
In addition, there exists δ > 0 such that the following holds true. Any triplet (c, and also enjoys the estimate
The factor p −1 in (2.39) is used to compensate the growing number of terms as p → ∞. In particular, we can view this as a uniqueness result with respect to a scaled L 2 -norm that will be specified later.

Nonuniqueness and numerical examples
Fixing r ∈ [r − δ r , r + δ r ], M = p q ≥ M * and θ ∈ R, the travelling wave (c M (θ, r), U M (θ, r)) is constructed as a perturbation of the travelling wave (c 0 , U 0 (· + θ)) on the domain p −1 Z. Since the wave profiles U 0 (· + θ) and U 0 (· + θ + p −1 ) are simply translates of each other on this domain, the shift-periodicity (2.37) follows easily. However, it is not clear how, specifically, the travelling wave depends on θ. Indeed, in [40, §5], Hupkes and Van Vleck show that it is reasonable to expect that the derivative ∂ θ c M (θ, r) is exponentially small in M . As such, it is unclear how to further analyse this dependence.
We emphasize that in general the travelling wave solution will not necessarily be unique, even up to translation. In particular, fixing θ ∈ (0, p −1 ), we note that the waves U 0 and U 0 (· + θ) are different on the domain p −1 Z. One might be tempted to conclude that if M is sufficiently large, the wave profiles U M (0, r) and U M (θ, r) are different as well. However, a larger value of M means that the grid p −1 Z becomes finer. In particular, since the travelling waves U M (0, r) and U M (θ, r) are perturbations of the waves U 0 and U 0 (· + θ), it could be that these perturbations cancel out the difference between U 0 and U 0 (· + θ).
In addition, since the constant M = (c∆t) −1 is fixed in the statement of Theorem 2.1, fluctuations in c automatically lead to changes in ∆t. This complicates our understanding of the fully discrete system for a fixed timestep ∆t > 0. Our main goal here is to show that the wavespeed c and the detuning parameter r do not depend on each other in a locally unique fashion, which is in major contrast to the corresponding continuous and semi-discrete systems.
However, the lack of a comparison principle for FitzHugh-Nagumo systems heavily complicates a direct analysis. As such, we have chosen to, instead, use numerical simulations to illustrate these phenomena. In particular, we focus on the backward-Euler discretisation of the FitzHugh-Nagumo MFDE, which takes the form These simulations turned out to be rather delicate, since the quality of the initial condition heavily influenced whether a solution could be found. In many cases, the simulation returned the zero solution. Simply augmenting an extra nontriviality condition often produced no solution at all. In addition, the value of c greatly determines the number of points ξ ∈ R for which the values (u, v)(ξ) need to be determined. In particular, upon writing  case, our simulations clearly show that the parameters c and r depend on each other in an intricate fashion. In particular, our results suggest that travelling wave solutions to the system (2.40) are not unique, since we were able to find solutions with a range of different wavespeeds at the same value for r. We refer to [10] and [49] for the corresponding dependence for the FitzHugh-Nagumo PDE and LDE respectively. In both cases, this dependence is given by a curve in the (c, r)-plane that resembles the symbol ∩.

Setup
The fully discrete travelling wave equation (2.29) is a highly singular perturbation of the semi-discrete travelling wave MFDE (2.11), which is the key complication for our analysis. In order to tackle this issue, we start by studying the linear operators that arise when linearizing the fully discrete travelling wave equation (2.29) around the semi-discrete travelling wave (c 0 , U 0 ). In particular, we define the linear expressions Our aim is to establish that the operators L k,M inherit several useful properties from the operator L 0 defined in (2.13) in the small timestep limit ∆t ↓ 0.
In this section we summarize and adept the setup from [40], sticking to the same notation as much as possible. In order to formulate our results, we need to define several function spaces. For any η ∈ R, we write In addition, given a Hilbert space H and any µ > 0, we define the corresponding sequence space   For now, we fix two integers q ≥ 1 and 1 ≤ k ≤ 6, together with a constant M = p q ∈ M q . To streamline our notation, we write Y M to refer to the space 2 p (R d ), i.e., Moreover, we introduce the space Y 1 k,M , which differs from Y M only by its inner product. To be more precise, we write (3.6) In addition, for f ∈ BC −η (R; R d ) with η > 0, we write π Y M for the sequence If moreover f ∈ BC 1 −η (R; R d ) and we wish to be explicit, we often write π Y 1 k,M f to refer to the restriction (3.7). The restriction operators π Y M and π Y 1 k,M are bounded, see Lemma A.1.
We can now consider the operators L k,M appearing in (3.1) as bounded linear maps Our goal is to define new sequence spaces, which allow us to pass to the limit M → ∞ in a controlled fashion. The basic idea is to use L 2 -interpolants for functions in Y M and H 1 -interpolants for functions in Y 1 k,M , so that the sequences in these spaces can be compared regardless of the different values of M . The main difficulty is to control terms of the form q , which is impossible to extract solely from the behaviour of D k,M v.
To tackle this issue, we need to perform q separate interpolations. Each of these interpolations must bridge a gap of size M −1 = q p . In particular, upon fixing an integer q ≥ 1 and writing Z q = {0, 1, 2, ..., q}, we introduce the space equipped with the inner product for M = p q ∈ M q , which both act as Our goal is to interpret L k,M as a map from H 1 k,M into H M . To this end, we pick n ∈ Z and 0 < ϑ ≤ 1 in such a way that 1 = (n + ϑ)M −1 . (3.20) Since M = p q ∈ M q , we see that ϑ = p−nq q , which yields In fact, because gcd(p, q) = 1, it follows that gcd(p, ϑq) = 1.
With these preparations in hand, we now write K k,M : H 1 k,M → H M for the linear operator that acts as for ζ ∈ q −1 Z q and ξ ∈ M −1 Z. Here the operator ∆ M is given by where we have introduced the twist operator T M : H M → H M that acts as taking into account the convention In particular, we see that the shift ϑ acts as a rotation number, connecting the different components of Φ in the ζ-direction. The inequality for Φ ∈ H M is almost trivial to verify in the finite-range setting, but turns out to be much harder to establish when dealing with infinite-range interactions; see Lemma A.5.
Finally, we introduce the notation to refer to the multiplication operator In fact, it is easy to see that which shows that K k,M and L k,M are equivalent.
Since the operator K k,M is not self-adjoint, we need to introduce the formal adjoint K * k,M : together with the map which constructs a function π ⊥ f ∈ L 2 (R, 2 q,⊥;∞ ) from a function f ∈ L 2 (R; R d ).

The limiting system
Our goal here is to exploit our understanding of the operator L 0 in order to determine the Fredholm properties of the limiting operator K q,ϑ . Due to the lack of a comparison principle we cannot immediately appeal to a general Frobenius-Peron-type result as was possible in [40]. The theory in this section aims to fill these gaps and can be considered the key technical contribution of this paper. We collect the main results in the following Proposition, which plays an essential role in Lemma 5.3 below.
Proposition 4.1 (cf. [40,Lem. 3.6]). Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Fix an integer q ≥ 1, together with a constant ϑ ∈ q −1 Z q that has gcd(ϑq, q) = 1. Then the operators K q,ϑ and K * q,ϑ are both Fredholm operators with index 0 and we have the identities Moreover, recalling the constantλ appearing in (HW2 r ), the operator K q,ϑ + λ is invertible for each λ ∈ C that has Re λ ≥ −λ and λ / ∈ 2πic 0 q −1 Z. Finally, there exists constants C > 0 and δ 0 > 0 so that for each 0 < δ < δ 0 and each Θ ∈ L 2 (R, 2 q,⊥;∞ ) we have the bound The first step towards proving Proposition 4.1 is to find the eigenvalues of the operator K q,ϑ . After that, we will focus on the essential spectrum of this operator. The idea behind the proof of Lemma 4.2 below can best be illustrated by considering the case q = 2. In this case, we have ϑ = 1 2 , together with In particular, if Θ is in the kernel of K 2, 1 2 + λ, the functions X 0 (ξ) = [Π 0 Θ](ξ), are eigenfunctions of the operator L 0 with eigenvalues −λ and −λ − c 0 πi respectively. Since −λ and −λ − c 0 πi cannot both be eigenvalues of L 0 at the same time in view of (HW2 r ), this means that at least one of the functions X 0 or X 1 is identically 0.
Without loss, we assume that X 0 = 0. In this case, the function Θ can explicitly be identified as As such, the eigenfunctions of K q,ϑ can be expressed in terms of those of L 0 , thus providing an upper bound on the dimension of the corresponding eigenspace. In addition, we have the identity Proof. Fix λ ∈ C with Re λ ≥ −λ. Suppose that Θ is in the kernel of the operator K q,ϑ + λ. For n ∈ {0, ..., q − 1} we set (4.11) with ζ q = exp[2πi/q] the q-th root of unity. Recalling that gcd(ϑq, q) = 1, it follows that this sum contains each of the functions Θ 0, ξ , ..., Θ (q − 1)q −1 , ξ exactly once. Recalling the definitions of the operators T 0 and T q,ϑ from (2.15) and (3.37), we can compute This allows us to obtain the identity (4.14) Suppose first that λ / ∈ 2c 0 πiq −1 Z. Then it follows from (HW2 r ) that −2c 0 πinq −1 − λ is no eigenvalue of L 0 for all 0 ≤ n ≤ q − 1. In particular, we must have X n = 0 for all 0 ≤ n ≤ q − 1. This means that the functions Π n Θ for 0 ≤ n ≤ q − 1 are also identically 0. Since the q × q Vandermonde matrix Z given by Z n,n = ζ n·n q is invertible, we obtain Θ(nϑ, ·) = 0 for all 0 ≤ n ≤ q − 1 from which (4.8) follows.
Proof. For fixed y ∈ R we introduce the matrix By decreasingλ if necessary, we can assume that −DG ρ − DG T ρ + Re λ is positive definite. It follows that X is the sum of a positive semi-definite matrix and a positive definite matrix and as such, it is positive definite itself. As a consequence, ∆ ρ;λ is positive definite as well and hence we obtain det ∆ ρ;λ (iy) = 0. Lemma 4.4. Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HW1 r ) and (HW2 r ) are satisfied. Assume that the triplet (G, P − , P + ) satisfies (HS3 r (b)). Pick λ ∈ C with Re λ > −λ and 0 ≤ ρ ≤ 1. Then we have det ∆ ρ;λ (iy) = 0 for all y ∈ R.

Linear theory for ∆t → 0
In this section, we apply the spectral convergence method to lift the Fredholm properties of the semi-discrete system to the fully discrete system in the small timestep limit ∆t → 0. In particular, we establish the main result below, which gives a quasi-inverse for the operators L k,M . This turns out to be the key ingredient in the construction of the discrete waves, which can subsequently be proved by means of a standard fixed point argument.
so that for all Ψ ∈ Y M the pair is the unique solution to the problem that satisfies the normalisation condition In addition, for all Ψ ∈ Y M we have the bound Proof of Theorem 2.1. On account of Proposition 5.1, the procedure outlined in [40, §4.1] can be followed to arrive at the desired result.
In order to facilitate the reading, we first outline our strategy and formulate two intermediate results in §5.1. This strategy heavily follows the program in [40], allowing us to simply refer to these results in many cases. However, due to the lack of a comparison principle and the many cross-terms we need to control, there are several key points in the analysis that need a fully new approach, which we develop in §5.2. In addition, the infinite-range setting forces us to obtain an extra order of regularity on the operator (L 0 + δ) −1 , which we achieve in §5.3.
The key step towards proving Proposition 5.1 is the establishment of lower bounds for these quantities. This procedure is based on [2,Lem. 6]. Our strategy to prove it is essentially the same, but some major modifications are needed to incorporate the difficulties arising from the discrete derivatives.
Proposition 5.2 (cf. [40,Prop. 3.7]). Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Fix a pair of integers 1 ≤ k ≤ 6 and q ≥ 1. Then there exists κ > 0 such that for all 0 < δ < δ 0 we have We are now ready to start our interpolation procedure. For any ξ ∈ R, we pick two quantities ξ ± M (ξ) ∈ M −1 Z in such a way that Using these quantities, we can define two interpolation operators that act as for all ζ ∈ q −1 Z q and all ξ ∈ R. These operators can be seen as interpolations of order zero and one respectively, both acting only on the second coordinate of φ. We refer to [40, for some useful estimates involving these interpolations.
With these preparations in hand, we start the proof of Proposition 5.2 using the methods described in the proof of [2,Lem. 6]. We focus on the quantity κ(δ) defined in (5.7), noting that κ * (δ) can be treated in a similar fashion. In particular, we find a lower bound for κ(δ) by constructing sequences that minimize this quantity. At this point it becomes clear why we work on the spaces H 1 (R, 2 q,⊥ ) and L 2 (R, 2 q,⊥ ), as we exploit the fact that bounded closed subsets of these spaces are weakly compact. Lemma 5.3 (cf. [40,). Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Fix a pair of integers 1 ≤ k ≤ 6 and q ≥ 1, as well as 0 < δ < δ 0 . Then there exist two functions Φ * ∈ H 1 (R, 2 q,⊥;∞ ), Ψ * ∈ L 2 (R, 2 q,⊥;∞ ), (5.12) together with three sequences and two constants ϑ ∈ q −1 Z q \ {0} and K 1 > 0 that satisfy the following properties.
(ii) The identity holds for all j ∈ N.
(iii) Recalling the constant κ(δ) defined in (5.7), we have the limit (iv) As j → ∞, we have the weak convergences (v) For any compact interval I ⊂ R, we have the strong convergences as j → ∞.
(vi) The function Φ * is a weak solution to (K q,ϑ + δ)Φ * = Ψ * and we have the bound In order to prove Proposition 5.2, we need to establish a lower bound on the norm Φ * H 1 (R, 2 q,⊥;∞ ) on account of (5.18). In Proposition 5.4 we follow the approach of [40,Lem. 3.18] in order to obtain this lower bound. Here we have to deal with both the cross-terms arising from the system setting as well as the infinite-range interactions.
Using this decomposition, we can expand the inner product as where we have introduced the cross-terms Recalling DG [1,2] = −Γ(DG [2,1] ) T and exploiting the identity we can rewrite the cross-terms to obtain which yields the desired bound.
Lemma 5.9. Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Then for each 0 < δ < δ 0 and each G ∈ H 1 (R; R d ) we have Then we can rewrite the equation (L 0 + δ)F = G in the form From this representation it immediately follows that F ∈ L ∞ (R; R d ). Differentiating both sides yields for n ∈ Z >0 , we can compute This allows us to estimate In particular, the sequence {F n } converges uniformly to ∆ 0 F , from which it follows that Since F, G ∈ H 1 (R; R d ), these considerations yield that F ∈ L 2 (R; R d ), from which the desired result follows.
We now turn to the desired exponential decay. The assumptions (HW1 r ) and (HW2 r ) yield the following useful properties of the operator L 0 . Lemma 5.10. Assume that (HS1) and (HS2) are satisfied and pick r in such a way that (HS3 r ), (HW1 r ) and (HW2 r ) are satisfied. Then the following properties hold for the LDE (2.1) with r = r.
(ii) Upon introducing the spaces In addition, there exists a constantη > 0 in such a way that for each 0 < η <η the map L 0 maps Proof. The proof of the statements (i)-(ii) follows the procedure described [58,Lem. 4.15, 6.8, 6.9] and will hence be omitted. It hence suffices to prove that ∆ 0 maps BC −η (R; R d diff ) into itself for η small enough. Upon picking F ∈ BC −η (R; R d diff ) and K ∈ R >0 in such a way that the bound in the space X 0 , which is given explicitly by The proof of [58,Prop. 5.2] provides the representation for each 0 < δ < δ 0 and each G ∈ L 2 (R; R d ). In addition, we can use Lemma 5.10 to pick constants K > 0 andα > 0 in such a way that Φ + 0 (y) dy (5.72) holds for each G ∈ L 2 (R; R d ). The following three results use (5.70) and (5.72) to establish the desired pointwise bound for (L 0 + δ) −1 .
Lemma 5.12. Consider the setting of Lemma 5.11. Then there exist constants 0 < δ * ≤ δ 0 and K > 0 so that for each 0 < δ < δ * and each G ∈ BC 1 −η (R; R d ) we have the bound For n ∈ Z >0 a calculation similar to (5.74) yields Using [58, Lem. 6.6] and Lemma 5.11, we obtain (5.78) Continuing in this fashion we see that the estimate holds for all n ∈ Z >0 and for some constant K 2 > 0. If we set then for each n ∈ Z >0 and each 0 < δ < δ * we have In particular, it follows that Since H 1 (R; R d )-convergence implies pointwise convergence we see that (5.83) Corollary 5.13. Consider the setting of Lemma 5.12. There exists a constant K > 0 so that for each 0 < δ < δ * and each G ∈ BC 1 −η (R; R d ) we have the bound (5.84) Proof. Fix 0 < δ < δ * and G ∈ BC 1 −η (R; R d ). Write F = (L 0 + δ) −1 G. The representation (5.70) together with Lemma 5.12 immediately yields the bound

A Auxiliary results
In this section we collect several useful results that we use throughout this paper. The first three results concern the sequence spaces Y M and Y 1 k,M and their associated inner products (3.5)-(3.6). Lemma A.1 ([40, Lem. 3.1]). Fix a pair of integers 1 ≤ k ≤ 6 and q ≥ 1, together with a constant η > 0. Then there exists a constant C ≥ 1 for which the bounds hold for all M ∈ M q and all functions f ∈ BC −η (R; R) and g ∈ BC 1 −η (R; R). Lemma A.2 ([40, Lem. 3.4]). Fix an integer q ≥ 1. Then there exists C > 1 so that the bound holds for all M ∈ M q and all functions f, g ∈ BC 1 −η (R; R d ). that hold for smooth, localized functions u. When dealing solely with nearest-neighbour interactions as in [40] the inequality ∆ M Φ, Φ H M ≤ 0 follows immediately from the Cauchy-Schwarz inequality. However, in our setting, some of the coefficients α k may not be positive definite, preventing us from taking them out of the inner product. This motivates the indirect approach that is taken in the proof of Lemma A.5. and observing that 1 = p q p q −1 = (p − nq) + nq M −1 q −1 , (A.8) we may computẽ for arbitrary ξ ∈ M −1 Z, ζ ∈ q −1 Z • q ∪ {0} and m ∈ Z. In particular, for m ∈ Z we obtain the identity (A.11) The desired result now follows from [2,Lem. 3].
We now show that K k,M approaches K q,ϑ in a more rigorous fashion. The infinite-range interactions cause complications here, because we need to interchange a limit and an infinite sum. For M ∈ M q we introduce the notation ϑ(M ) to refer to the value of ϑ in (3.20) with M =M .