Travelling Waves for Adaptive Grid Discretizations of Reaction Diffusion Systems III: Nonlinear Theory

In this paper we consider a spatial discretization scheme with an adaptive grid for the Nagumo PDE and establish the existence of travelling waves. In particular, we consider the time dependent spatial mesh adaptation method that aims to equidistribute the arclength of the solution under consideration. We assume that this equidistribution is strictly enforced, which leads to the non-local problem with infinite range interactions that we derived in Hupkes and Van Vleck (J Dyn Differ Eqn, 2021). Using the Fredholm theory developed in Hupkes and Van Vleck (J Dyn Differ Eqn, 2021) we setup a fixed point procedure that enables the travelling PDE waves to be lifted to our spatially discrete setting.


Introduction
Our goal in this paper is to complete the program initiated in [26,27] to analyze the impact of adaptive discretization schemes on scalar bistable reaction-diffusion PDEs of the form u t = u xx + g(u). (1.1) In particular, for any discretization distance h > 0 and any j ∈ Z, we write x jh (t) for the time-dependent location of the relevant gridpoint and U jh (t) for the associated approximation for u x jh (t), t . We then study the systeṁ x (j+1)h (t)−x jh (t) + g U jh (t) , (1.2) in which x(t) is defined implicitly by demanding that and imposing the boundary constraint lim j→−∞ [x jh (t) − jh] = 0. (1.4) This means that mesh locations are determined dynamically in time so that there is equidistribution of a finite difference approximation of the arclength of the solution of (1.2). We show that this system has solutions of the form U jh (t) = Φ(x jh (t) + ct), (1.5) which can be interpreted as travelling waves. For concreteness, we will use the cubic nonlinearity g(u) = g cub (u; a) = u(1 − u)(u − a), 0 < a < 1 (1.6) throughout this introduction to explain the main ideas.
Travelling waves The pair (Φ, c) that we construct will be close to the travelling wave (Φ * , c * ) for the PDE (1.1). Using (1.6), this pair must satisfy the travelling wave ODE cΦ * = Φ * + g cub (Φ * ; a), Φ * (−∞) = 0, Φ * (+∞) = 1. (1.7) Such solutions provide a mechanism through which the fitter biological species (corresponding to the deepest well of the potential − g cub ) can become dominant throughout a spatial domain. For this reason they are sometimes referred to as invasion waves.
It is well-known that these waves play an important role in the global dynamics of (1.1). For example, using the comparison principle one can show that these waves are nonlinearly stable under a large class of perturbations [14] and that they determine the spreading speed of localized structures [39]. In addition, they have been used extensively as building blocks to construct general time dependent solutions of reaction-diffusion systems. For example, planar versions of (1.1) support (sharp) travelling corners [4,16], expansion waves [33], scattering waves [3] and modulated waves [9] that connect periodic travelling waves of nearby frequencies.
Uniform spatial discretizations In order to set the stage, let us return to the lattice differential equation (LDE)U j (t) = 1 h 2 [U j−1 (t) + U j+1 (t) − 2U j (t)] + g cub U j (t); a , (1.8) which can be used to describe the uniformly discretized approximants U j (t) ∼ u(jh, t). Mathematically speaking, the transition from (1.1) to (1.8) breaks the continuous translational symmetry of the underlying space. Indeed, (1.8) merely admits the discrete group of symmetries j → j + k with k ∈ Z. As a consequence, travelling wave solutions U j (t) = Φ(jh + ct) (1.9) can no longer be seen as equilibria in an appropriate comoving frame. Instead, they must be treated as periodic solutions modulo the discrete shift symmetry discussed above. The resulting challenges occur frequently in similar discrete settings and general techniques have been developed to overcome them [2,7,15]. Direct substitution of (1.9) into (1.8) yields the travelling wave equation to which we again append the boundary conditions Φ(−∞) = 0, Φ(+∞) = 1. (1.11) Due to the presence of the shifted arguments such equations are known as functional differential equations of mixed type (MFDEs). Note that the unbounded second derivative operator in (1.7) has been replaced by a bounded second-difference operator. In addition, the transition c → 0 is now singular since it changes the structure of the equation. As a consequence, there is a fundamental difference between standing and moving wave solutions to (1.8).
In the anti-continuum regime h 1, the second-difference operator can be treated as a small perturbation to the remaining ODE. An elegant construction pioneered by Keener [28] allows one to construct standing waves for a = 1 2 that satisfy the boundary conditions (1.11) and block the two stable background states Φ ≡ 0 and Φ ≡ 1 from invading the domain. In particular, the simple geometric condition [27,Eq. (1.5)] is violated in this setting. This phenomenon is often referred to as pinning or propagation failure and has attracted a considerable amount of attention [1,8,10,11,18,22].
In the intermediate h ∼ 1 regime the shifted terms cannot be handled so easily and one needs to understand the full MFDE. Such equations are ill-posed as initial value problems and hence must be handled delicately. Several important tools have been developed to accomplish this, such as Fredholm theory [29] and exponential dichotomies [17,31,34,35]; see [21] for a detailed overview.
Using a global homotopy argument together with the comparison principle, Mallet-Paret [30] constructed a branch of solutions Φ(a), c(a) to (1.10) with (1.11), in which c(a) is unique and Φ(a) is unique up to translation when c(a) = 0. For the uniform spatial discretization of the FitzHugh-Nagumo PDE [26,Eq. (1.3)], a generalization of Lin's method can be used to establish a version of the exchange Lemma for MFDEs and construct stable travelling pulses [23,24]. Further results in this area can be found in [5,6,28,30,40] and the survey [21].
Uniform spatial-temporal discretizations The so-called backward differentiation formula (BDF) are a family of six schemes that can be used for discretizing the temporal derivative in (1.8). These are well-known multistep methods that are appropriate for parabolic PDEs due to their numerical stability properties. As an illustration, we note that the two lowest order schemes prescribe the substitutionsU j (t) → 1 ∆t U j n∆t − U j (n − 1)∆t , U j (t) → 1 2∆t 3U j n∆t − 4U j (n − 1)∆t + U j (n − 2)∆t , (1.12) in which ∆t > 0 denotes the timestep and t = n∆t. The first scheme is also known as the backward Euler method and has the advantage that is preserves the comparison principle, unlike the five other members of the family.
In [25] we constructed fully discretized travelling wave solutions U j (n∆t) = Φ(j + nc∆t), Φ(−∞) = 0, Φ(+∞) = 1 (1.13) for the coupled map lattices arising from these discretization schemes. The relevant travelling wave equations for the two lowest order schemes can be obtained by making the replacements (1.14) in the MFDE (1.10). In the first case we leveraged the comparison principle to obtain global results. We established that the c(a) relation can become multi-valued, which clearly distinguishes the fullydiscrete regime from its spatially-discrete counterpart. The same behaviour occurs for the five other BDF methods, but here we only have results for small ∆t > 0. We remark that related phenomena have been observed in monostable KPP systems [32] in the presence of inhomogeneities.
These non-uniqueness results should be seen as part of the program that was initiated in [11][12][13] to study the impact of temporal and full discretization schemes on various reaction-diffusion systems. Indeed, these papers studied versions of (1.1) with various smooth and piecewise linear bistable nonlinearities. The authors used adhoc techniques to obtain rigorous, formal and first order information concerning the change in the dynamics of traveling wave solutions. In addition, in [7] the authors considered the forward-Euler scheme and used Poincaré return-maps and topological arguments to obtain the existence of fully-discretized waves.
To appreciate the advantage of this indirect approach, we note that any attempt to use ξ will lead to an equation for the waveprofile Φ with shifts that depend on the waveprofile Φ itself. In particular, the resulting wave equation is a state-dependent MFDE with infinite range interactions. At the moment, even state-dependent delay equations with a finite number of shifts are technically very challenging to analyze, requiring special care in the linearization procedure [38]. Indeed, linearizations typically involve higher order (continuous) derivatives, making it very hard to close fixed-point arguments.
Physical frame It turns out that there is a close relation between the two wave Ansatzes (1.5) and (1.19). In order to see this, let us assume for the moment that we have found a triplet (Φ, c, x) for which x and the function U defined in (1.5) satisfy (1.2) together with (1.3)- (1.4). Let us also assume that for each ϑ ∈ R there is a unique increasing sequence y jh;ϑ with y 0;ϑ = ϑ for which holds for all j ∈ Z. This can be arranged by imposing a-priori Lipschitz bounds on Φ and Φ and picking h > 0 to be sufficiently small. Finally, let us assume for definiteness that c > 0 and that the wave outruns the grid in the sense thatẋ 0 (t) + c > > 0. A direct consequence of this inequality is that The uniqueness property discussed above hence implies that (1. 25) we see that in fact for all j ∈ Z. Taking the limit j → −∞, the boundary conditions (1.4) imply that cT = h. Exploiting the well-posedness of our dynamics in forward and backward time, we conclude that holds for all j ∈ Z and t ∈ R. Writing Ψ x (ϑ) = x 0 (ϑ/c), we hence find for all j ∈ Z and t ∈ R. Upon introducing the function this allows us to obtain the representation In fact, we show that for arbitrary solutions U to (1.15) for which U (t 0 ) is close to Ψ * (hZ+ϑ), we indeed have the pointwise inequalities |ẋ(t 0 )| < |c| whenever c is sufficiently close to c * . This can be used to show that the coordinate transformation (1.30) can be inverted, allowing us to reconstruct the profile Φ(ξ) from Ψ U (τ ).
Fixed-point setup In order to construct our travelling waves, we write Ψ = Ψ * +v and decompose (1.20) into the form (1.32) using the linear operators L h that were introduced in [27]. In the limit h ↓ 0, these operators reduce formally to the operator L * associated to the linearization (1.17) around the wave (1.18). The singular nature of the transition between (1.15) and (1.17) is fully encoded in the transition between L * and L h , which was studied at length in [27]. As a result, our analysis in this paper can be seen as the construction of a regular fixed point problem. However, there are two main challenges that need to be overcome. The first complication is that L h is not the 'exact' linearization of G, which is far too complicated to handle. Instead, we recover our operator L h after several simplification steps, which each introduce h-dependent errors that need to kept under control. In order to achieve this, we reapply the approximation framework developed in [26] in order to systematically bound the global errors that arise by modifying the individual factors of the products that appear in the definition of G.
The second obstacle is that the nonlinearity G acts on sequences U : hZ → R, while the fixed-point problem (1.32) is formulated in terms of functions v ∈ L 2 . Since the relevant transitions between supremum and L 2 -based norms cost a factor of h −1/2 , special care must be taken to construct appropriate function spaces that allow uniform bounds for h ↓ 0. This is particularly dangerous for the terms that are quadratic in the second differences of U , which correspond roughly to the − γ −4 w 2 θθ term in (1.17). In fact, we need to exploit the special structure of G and take a discrete derivative of (1.32) in order to close our problem. We hence need to obtain estimates on the discrete derivative of the nonlinear residual G nl , which requires an elaborate bookkeeping system.
Outlook We view our work here as a first step towards understanding the impact of adaptive discretization schemes on travelling waves and other patterns that exist for all time. In particular, we believe that the waves constructed here can be seen as a slow manifold for the dynamics of the full system (1.2) with the non-instantaneous gridpoint behaviour prescribed by the MMPDE5 scheme [20]. Here σ > 0 is a tunable speed parameter, which we effectively set to zero by passing to (1.3).
Using the Fredholm theory developed in [27] for the operators L h one should be able to leverage the ideas in [37] to effectively track the fast grid-dynamics in the 0 < σ 1 regime. A further step in the program would be to also handle temporal discretizations, inspired by the approach developed in [25] that we described above. Finally, we feel that it is important to understand the stability of the discretized waves under the full dynamics of the numerical scheme. To achieve this, one could follow the approach in [36] to transfer information from the operators L h to the linearization around the actual adaptive travelling waves constructed in this paper.
We are specially interested here in the pinning phenomenon. Indeed, numerical observations indicate that the set of detuning parameters a for which c(a) = 0 shrinks dramatically when using adaptive discretizations. Understanding this in a rigorous fashion would give considerable insight into the theoretical benefits of adaptive grids compared to the practical benefits of increased performance. Preliminary results in this direction can be found in [19].
Let us emphasize that the application range of our techniques does not appear to be restricted to the scalar problem (1.1) or the specific grid-update scheme (1.33) that we use. Indeed, using the framework developed in [36], it should be possible to perform a similar analysis for the FitzHugh-Nagumo equation PDE and other multi-component reaction-diffusion problems. In addition, any numerical scheme based on the arclength monitor function will share (1.3) as the instantaneous equidistribution limit.
Overview This paper is organized as follows. After formulating our main results in §2, we introduce our notational framework and recap the key contributions from [26,27] in §3. In §4 we simplify the nonlinear functions that appear as factors in the product structure of G and obtain estimates on all the resulting errors. These estimates are used in §5 to compute tractable expressions for the linearization of G and its discrete derivative G + around Ψ * and obtain bounds on the residuals. We conclude in §6 by combining all these ingredients with the theory developed in [27]. In particular, we develop an appropriate fixed-point argument to construct our desired travelling waves.
In order to develop the main story in a reasonably streamlined fashion that focuses on the key ideas, we have chosen to transfer many of the tedious underlying estimates and algebraic manipulations to the appendices. In order to keep this paper as self-contained as possible, these appendices also summarize some of the fundamental auxiliary bounds that were obtained in [26,27].
Acknowledgements. Hupkes acknowledges support from the Netherlands Organization for Scientific Research (NWO) (grant 639.032.612). Van Vleck acknowledges support from the NSF (DMS-1419047 and DMS-1714195). Both authors wish to thank W. Huang for helpful discussions during the conception and writing of this paper and an anonymous referee for providing valuable feedback.

Main results
The main results of this paper concern adaptive-grid discretizations of the scalar PDE u t = u xx + g(u). (2.1) Throughout the paper, we assume that the nonlinearity g satisfies the following standard bistability condition.
(Hg) The nonlinearity g : R → R is C 3 -smooth and has a bistable structure, in the sense that there exists a constant 0 < a < 1 such that we have It is well-known that the PDE (2.1) admits a travelling wave solution that connects the two stable equilibria of g [14]. The key requirement in our next assumption is that this wave is not stationary, which can be arranged by demanding 1 0 g(u) du = 0. (HΦ * ) There exists a wave speed c * = 0 and a profile Φ * ∈ C 5 (R, R) that satisfies the limits and yields a solution to the PDE (2.1) upon writing In [26] we derived an effective equation for the dynamics of the sequence U (t) : hZ → R featuring in the adaptive scheme (1.2)-(1.4) for (2.1) that no longer explictly depends on the location of the gridpoints. In order to formulate this reduced equation, we recall the discrete derivatives together with the first-order differences and the second order analogues This allows us to introduce the auxiliary functions and subsequently write These ingredients allow us to formulate the effective reduced system [26,Eq. (2.25)] for the dynamics of (1.2)-(1.4) asU (t) = G U (t) , (2.13) which will be the main system that we analyze in this paper. We recall the arclength parametrization ξ * (τ ) defined by the identity together with the stretched waveprofile Ψ * : R → R given by The main result of this paper states that for sufficiently small h > 0, the reduced problem (2.13) admits a travelling wave solution in an appropriate sense. These waves are locally unique up to translation. We note that items (iv) and (v) use the notation for functions v. In addition, we use the shorthands L 2 = L 2 (R; R) and H 1 = H 1 (R; R), together with the Heaviside sequence H jh = 1 j≥0 .
In addition, the identity (2.13) and the strict inequality ∂ + U (t) ∞ < 1 both hold for all t ∈ R.
(v) Pick any 0 < h ≤ δ h and consider a pair

25)
together with the strict inequality ∂ + U ∞ < 1 for all t ∈ R. In addition, if U is a solution to the system (2.13) for all t ∈ R, then we must have for some ϑ ∈ R.
We emphasize that the location of the gridpoints for the waves (2.16) can be determined by using that satisfy the following properties.
(i) Upon writing (ii) For every t ∈ R and j ∈ Z, the functions defined in (2.29) satisfy the relation We remark that if (2.16) and (2.30) both hold, simple substitutions yield the identity (2.31) In particular, the main assertion in Corollary 2.2 is that the perturbed coordinate transformation is invertible for sufficiently small h > 0, allowing us to transfer the waves back to the original physical framework.

Setup and notation
In this section we recall several crucial results and notational conventions introduced in the prequel papers [26,27]. This will ensure that the current paper can be read reasonably independently. As a preparation, we recall the sequence spaces that were introduced in [26, §3.3], together with the higher order norms and their counterparts For a single fixed h > 0 all these norms are naturally equivalent to the 2 h -norm or the ∞ h -norm. The point here is that the h −1 factor in the definition of ∂ + introduces a natural scaling that will allow us to formulate h-independent bounds.
In addition, we pick a reference function U ref; * ∈ C 2 (R, [0, 1]) that satisfies the properties For any κ > 0, this allows us write and introduce an open subset We can now recall the affine subset [26, §3.4] that captures the admissable states of the waves that we are interested in and provides adequate control on the necessary difference operators. Indeed, for each U ∈ Ω h;κ we have the crucial bound ∂ + U ∞ h ≤ 1 − κ, which ensures that our grid points are well-defined. In addition, the norms ∂ + U 2,1 h , U ∞;2 h and g(U ) 2 h are all bounded uniformly in h > 0. Finally, it is possible to pick 0 > 0 and κ > 0 in such a way that for any 0 < h < 1 and any v ∈ H 1 that has

Linear operators
In [27] we analyzed several important linear operators that will turn out to be closely related to our travelling wave system (1.20). To set the stage, we recall the sequences  (3.11) together with their twisted counterparts L U : 2 h → 2 h defined by always taking U ∈ Ω h;κ . A special role is reserved for the discrete derivative ∂ + M U , which we approximate by the linear operator Based on the computations in [27,Prop. 5.5], this decomposition uses the expressions that feature at most second differences of V , together with which contains third differences that needs to be treated carefully. Several crucial bounds for these operators are collected in Proposition C.1.
We are now ready to recall the linear operator L h : H 1 → L 2 that acts as where we are slightly abusing notation. Indeed, recalling the discrete evaluation operator that 'samples' a function f on the grid ϑ + hZ, the identity L h v = f should be interpreted as the statement that holds for each ϑ ∈ [0, h]. We remark that the right-hand side above is continuous in 2 h as a function of ϑ as a consequence of (A.6) and the continuity of the translation operator on H 1 .
The key purpose of [27] was to construct a quasi-inverse for the operator L h . Indeed, [27, Thm.

2.3] establishes the existence of two linear maps
defined for small h > 0, so that for each f ∈ L 2 the pair is the unique solution to the problem up to a normalization condition that can be used to fix the phase of our constructed wave. The crucial point is that we obtain h-uniform bounds which will provide the required control on the second and third differences of our travelling wave. We remark that these difference operators cannot be replaced by the corresponding derivatives, which forces us to develop a rather delicate fixed point argument in §6. In addition, we note that the spectral convergence framework used to obtain (3.22) relies strongly on the inner product structure of L 2 , which explains why we do not have L ∞ -based estimates (yet). This is the reason that we go to great lengths throughout the paper to work with 2 h -bounds as much as possible. Indeed, the results in §A show that these mix well with L 2 -functions, unlike supremum bounds.

Error functions
The errors that need to be controlled during our reduction steps arise from various sources that we briefly discuss here. As preparation, we recall the translation operators and the sums and products S ± a = 1 2 (a + T ± a), P ± a = aT ± a. These allow us to recall the function from [26,Eq. (7.28)], which measures the smoothness of U in some sense. Indeed, it becomes small whenever third differences of U can be controlled, which is the case when taking U = Ψ * . In a similar vein, we introduce the error functions which can be used to 'shift' function evaluations back and forth between neighbouring lattice sites. We note in general that overlined symbols will be used for expressions related to G + , which naturally involve higher order differences than those related to G. Indeed, the nonlinearities in our problem will be controlled by the product (3.28) Observe here that the supremum norms are always at least one order smaller than the highest 2 -based norms. In addition, there are no squares of third-differences or products involving only supremum bounds. These facts will turn out to be crucial when passing between sequences and functions in order to apply the estimates (3.22) in §6.
Our final error functions are given by together with its approximate first difference (3.30) The relation between these two functions is explored in Proposition C.2. We view both expressions as a measure for the difference between U and the stretched travelling wave Ψ * . Indeed, upon introducing the notation which resembles the continuum limit of (3.29). This can be differentiated to yield the natural limit of (3.30). Together with the smoothness term E sm , the functions (3.29) and (3.30) can be used to define our final remainder terms These are small when taking U = Ψ * and describe the additional error contributions generated in this paper that cannot be absorbed by the terms in [26].

Initial approximants for G and G +
The expression (2.12) for G(U ) is too convoluted for practical use, featuring third differences and double sums. It hence needs to be simplified, at the cost of introducing error terms. An initial step in this direction was performed in [26, Eq. (6.10)], where we decomposed G(U ) into a number of products featuring nonlinearities from the set which were all defined in [26, §6] and contain at most second differences. A similar decomposition was obtained for G + (U ) = ∂ + G(U ) in [26,Eq. (6.16)], but now with nonlinearities from the set together with an explicit third-difference term. In addition, for each of the nonlinearities f ∈ S nl;short we (implicitly) defined an approximation f apx (U ) and an approximate linearization f lin;U [V ] in [26, §8].
In fact, the full definitions of the nonlinearities f ∈ S nl;short turn out to be irrelevant for our purposes here, so there is no need to repeat them from [26, §6]. However, we do need to manipulate their approximations, which we therefore evaluate in full here by substituting the relevant expressions from [26, §7] into the definitions [26, Eq. (8.1)- (8.4)]. This yields the approximants together with the approximate linearizations The corresponding expressions for the two remaining second-difference operators can be copied from [26,Eq. (7.22)] and read for G(U ), featuring the four components    One of our main aims in [26] was to develop a framework to control the errors that arise by these types of approximations. In particular, we needed to track the propagation of errors on the individual factors of (2.12) through the full sums and exponents. The bounds in [26, Lems. 8.1-8.3] provide a constant K = K(κ) > 0 so that these errors satisfy for all U ∈ Ω h;κ . In addition, the nonlinear residuals can be estimated as for any U ∈ Ω h;κ and any V ∈ 2 h for which U + V ∈ Ω h;κ . These initial approximants for G and G + are already much easier to work with than (2.12) and enabled us to establish the well-posedness of our reduced system (2.13) in [26]. However, they are still unwieldy on account of the shifts and the sums. In addition, several simplifications can be made that only become apparent when looking at the full combinations (3.41), (3.43) and (3.45). This will be the main focus of §4- §5.
Convention Throughout the remainder of this paper, we use the convention that primed constants (such as C 1 , C 2 etc) that appear in proofs are positive and depend only on κ and the nonlinearity g, unless explicitly stated otherwise.

Component estimates
The first important task in this paper is to build a bridge between the linear theory described in §3.1 and the approximation framework outlined in §3. 3. This requires us to refine the approximants introduced in the latter section. We carry out the first step of this procedure here, focusing our attention on the nonlinearities introduced in (3.35)-(3.36).
In particular, for any f ∈ S nl;short we introduce the further decompositions The expressions with the label 'expl' are the actual explicit simplifications that will play a key role in our further computations. The label 'sh' is used for terms which are always small, which we will be able to absorb into the error terms E sh and E sh defined in (3.26). Finally, the label 'rem' is used for remainder terms that are small when using U = Ψ * . The explicit decompositions (4.1) are provided in §D. Our main contribution here is to summarize the errors that arise in a structured fashion that resembles the main spirit of the framework developed in [26, §7]. This will allow us to replace all the occurrences of f apx and f lin;U in §3.3 by their refinements f apx;expl and f lin;U ;expl , leading to a second round of approximations G apx;II , G + apx;II , G lin;U ;II and G + lin;U ;II . In order to achieve this, we define a preferred exponent set for each f ∈ S nl;short , together with its counterpart for each f ∈ S nl;short . This is done in such a way that we can write for a set of bounded multi-linear maps for all 1 ≤ j ≤ k i . Stated more informally, the 2 h norm of G apx;I (U ) can be bounded in terms of products of q 2 -norms of nonlinearities f ∈ S nl;short , where each q is taken from the preferred set of exponents. This is the direct analogue of [26,Cor. 6.4].
A short inspection of the products (3.42), (3.46) and (3.47) readily shows that there is some freedom as to which factors should be measured in 2 h . In fact, it is possible to put an 2 h norm on any chosen factor, at the price of possibly having to flip the exponent of a companion factor that has 2 ∈ Q f ;pref from two to infinity.
This freedom is essential to obtain sharp estimates and hence requires us to deviate from the preferred exponents from time to time. The main focus of [26, §5, §7-8] was to develop a bookkeeping framework to keep track of this procedure. We build on this investment here and follow the spirit of [26, §5.2] to define further exponent sets for each f ∈ S nl;short ∪ S nl;short . The first of these contains all values of q for which f apx maps into q h . On the other hand, the set Q f ;lin contains all q for which we need to evaluate the q h -norm of f lin;U ;expl and f lin;U ;sh , while Q f ;lin;rem contains these exponents for f lin;U ;rem .
In order to illustrate these points, let us consider the example which appears (after dropping a shift for notational clarity) as a factor in the component G B;lin;U ;I that needs to be evaluated in the supremum norm; see (E.17). Our goal is to simplify this expression by writing where we have introduced the terms We note here that the preferred exponent sets are defined in (D.12), (D. 19) and (D.30) and given by Recalling the remainder function introduced in (3.34), we may readily use these preferred exponents to compute Here we use property (4.33) below with q = ∞ (see (D.36)) together with the a-priori bounds (4.21). For the second term we can use the same properties, but now with q = 2 (see (D.36)). In particular, we obtain (4.14) Note that this required us to swap the first two exponents, which is made possible by the demand ∞ ∈ Q f for each f ∈ S nl;short ; see Proposition 4.1. This swap is also required for the final term, which can be controlled by Indeed, simply using ∞ h on the middle factor would lead to a contribution proportional to V ∞ h ; see the third line of (D.32). Such a term would lead to problems in §6 and hence is not contained in E rem;U . This scenario is covered in our structural results below by using options (b) from both Proposition 4.2 and 4.3. We feel that this relatively small example already clearly illustrates the benefits of utilizing an abstract bookkeeping scheme instead of direct estimates.

Summary of estimates
In order to state our results, we introduce the expressions and finally These expressions are all related to the size of the f apx functions and play a very similar role as the quantities S full and S 2;fix that were defined in [26, §7]. In particular, the 'full' terms correspond to all the exponents that we need to use, while the 'fix' expressions reflect the contributions that are only allowed to be evaluated in 2 h ; see (4.26). In addition, we recall the quantities that are associated to the approximate linearizations f lin . Here T ∞;opt represents the contributions where the use of the supremum norm is optional, in the sense that they could also be measured in 2 h . The remaining contributions are all reflected in T safe . We emphasize that the main point of our bookkeeping scheme is to ensure that products of the form S #;full T ∞;opt are never needed, where # ∈ {sh, rem, diff}. Our main results summarize the structure that the decompositions described in §D will adhere to. Propositions 4.1 and 4.2 state that the approximants f apx;# are all uniformly bounded and that the full linear approximants f lin;U share the structure and estimates of the nonlinearities in the sets S nl ∪ S nl analyzed in [26]. These can be interpreted as the counterparts of [  . For every f ∈ S nl;short we have ∞ ∈ Q f together with (4.20) In addition, there exists K > 0 so that for each q ∈ Q f , the bound holds for all h > 0 and U ∈ Ω h;κ . The same properties hold upon replacing (S nl;short , Q f ;pref ) by (S nl;short , Q f ;pref ).
. Assume that (Hg) is satisfied and fix 0 < κ < 1 12 . For any f ∈ S nl;short , any # ∈ {expl, sh, rem} and any q ∈ Q f ;pref , at least one of the following two properties hold true.
(a) There exists K > 0 so that holds for every h > 0, U ∈ Ω h;κ and V ∈ 2 h .
(b) We have q = ∞ and there exists K > 0 so that the bounds hold for every h > 0, U ∈ Ω h;κ and V ∈ 2 h . The same properties hold upon making the replacement (4.24) . Assume that (Hg) is satisfied and fix 0 < κ < 1 12 . Then there exists K > 0 so that for every f ∈ S nl;short , q ∈ Q f ;pref and # ∈ {sh, rem}, we have for any h > 0 and U ∈ Ω h;κ . In addition, if 2 ∈ Q f ;pref then for every # ∈ {sh, rem} at least one of the following two properties hold true.
(a) There exists K > 0 so that holds for every h > 0 and U ∈ Ω h;κ .
(b) There exists K > 0 so that holds for every h > 0 and U ∈ Ω h;κ .

Refined approximants for G and G +
We now introduce the expressions for each f ∈ S nl;short ∪ S nl;short . The full explicit forms can be found in §E-F, but are not important for our purposes here. We leave the expressions for A a intact and simply write  In particular, the results below describe the residuals that arise when replacing the initial approximants defined in (3.41), (3.43) and (3.45) by these refined versions. Since our bookkeeping framework has the same overall structure as in [26], we can reuse the analysis developed there in a streamlined fashion.
Lemma 4.7. Assume that (Hg) is satisfied and fix 0 < κ < 1 12 . There exists a constant K > 0 together with sequences defined for every h > 0 and U ∈ Ω h;κ , so that the following properties hold true.
(i) For every h > 0 and U ∈ Ω h;κ we have the identities   Recalling the general identity and its extensions, we write      Arguing as above, Proposition 4.6 yields

Estimates for G and G +
In this section we exploit the component estimates from §4 to analyze the function G defined in (2.12), together with its first difference G + . In particular, recalling the operator L U defined in (3.12), we introduce our final approximants (5.1) and write Using the discrete calculus outlined in §A, one may readily verify the identities Our main result quantifies the approximation errors in terms of the functions E sh;U , E rem;U and E prod and their counterparts E sh;U , E rem;U and E prod;U defined in (3.26), (3.27), (3.28) and (3.34). For convenience, we also reference the quantities (4.16)-(4.17).
Proposition 5.1. Suppose that (Hg) is satisfied and fix 0 < κ < 1 12 . Then there exists K > 0 so that the following properties hold.

Refinement strategy
We recall the refined approximants G apx;II (U ) and G lin;U ;II [V ] that we defined in §4.2. The main task in this section is to track the errors that accumulate as we reduce these expressions even further to our relatively simple approximants (5.1). In contrast to the abstract approach in §4, we achieve this in a direct fashion through several explicit computations. Indeed, in §E-F we obtain the following representations.
Proposition 5.2 (see §E-F). Assume that (Hg) is satisfied and fix 0 < κ < 1 12 . There exists a constant K > 0 together with sequences defined for every h > 0 and U ∈ Ω h;κ , so that the following properties hold true.
(i) For every h > 0 and U ∈ Ω h;κ we have the identities (5.14) (ii) For every h > 0, U ∈ Ω h;κ and V ∈ 2 h we have the bounds   We now turn towards the Lipschitz bounds for G nl .
Lemma 5.5. Assume that (Hg) is satisfied and pick 0 < κ < 1 12 . There exists a constant K > 0 so that the estimate Proof. By definition, we have In particular, we get In view of (5.18), we find

Travelling waves
Formally substituting the travelling wave Ansatz (2.20) into the reduced system (2.13) leads to the nonlocal differential equation cΨ = G(Ψ). (6.1) In this section we set out to construct solutions to this equation for small h > 0 that can be written as Ψ = Ψ * + v, c = c * +c (6.2) for pairs (c, v) that tend to zero as h ↓ 0. Care must be taken to ensure that the expression G(Ψ) is well-defined, but based on our preparations we are able to provide a relatively streamlined fixed-point argument here, which allows us to prove the results stated in §2.
In order to control the size of the perturbation (c, v) ∈ R × H 1 , we introduce the norms for h > 0 and write Z h for the set R × H 1 equipped with this new norm. Observe that for fixed h this norm is equivalent to the usual one on R × H 1 .
Recalling the discussion at the start of §3, we pick 0 < κ < 1 12 and 0 > 0 in such a way that the inclusion holds for all 0 < h < 1, all ϑ ∈ [0, h] and all v ∈ H 1 that satisfy (3.8). In order to accommodate this, we pick two parameters δ > 0 and δ + v > 0 and introduce the set Since ∂ + h is bounded on H 1 and L 2 for each fixed h, we note that this is a closed subset of Z h . Substituting (6.2) into (6.1), we obtain which should be interpreted in a sense similar to that of (3.18).

Upon introducing the nonlinearity
and inspecting the definitions (3.16) and (5.1), we can rewrite (6.6) as Recalling the two solution operators (3.19), we now introduce the map W h : which allows us to recast (6.8) as the fixed point problem In order to show that W h is a contraction mapping on Z h;δ,δ + v , we study the two expressions H h and c * Ψ * − G(Ψ * ) separately in our first results. The sampling bounds from §A play a key role here, as they enable us to extract L 2 -based bounds on G nl;Ψ * from the sequence estimates obtained in §5.
Notice that the control (6.15) on v 2;2 h would not have been possible using only bounds on (6.3), since L 2 -norms cannot directly be turned into 2 h -norms. This would prevent us from bounding the terms that are quadratic in these second differences. In fact, this is the reason that we needed to obtain such detailed bounds on G + in this series of papers. Indeed, the additional third-differences only appear in a linear fashion, which does allow us to easily pass between sequences and functions; see (A.18). Lemma 6.1. Suppose that (Hg) and (HΦ * ) are satisfied. There exists K > 0 so that for any pair (δ, δ + v ) ∈ (0, 1) 2 and any 0 < h < 1 the estimates (1) , v (2) − v (1) ) Z h (6.12) holds for each set of pairs (c (1) , v (1) ) ∈ Z h;δ,δ + v and (c (2) , v (2) ) ∈ Z h;δ,δ + v . Proof. The first term in H h can be handled by the elementary estimates together with (6.14) Using Corollary A.1 we see that for all (c, v) ∈ Z h;δ,δ + . For any ϑ ∈ R, we may hence exploit Propositions C.3 and 5.1 to obtain the estimate A second application of Corollary A.1 yields the bound For any ϑ ∈ R, we may hence use Propositions C.3 and 5.1 to find We now apply Lemma A.2 to obtain (6.20) Using (A.6) we note that

(6.21)
Applying Lemma A.2 once more, we obtain 22) The desired bounds follow readily from these estimates. Lemma 6.2. Suppose that (Hg) and (HΦ * ) are satisfied. There exists K > 0 so that for each 0 < h < 1 we have the bounds Proof. Applying Lemma A.2 together with Propositions C.3 and 5.1, we find We now compute together with (6.26) Applying (A.8) we see that as desired.
Utilizing the linear theory from [27] that we outlined in §3.1, we are now in a position to study the full nonlinear term W h . Our main result subsequently follows in a relatively standard fashion from the uniqueness properties of the contraction mapping theorem. Proof. Using the estimates (3.22) together with the a-priori bounds (h, δ, δ + v ) ∈ (0, 1) 3 , we obtain the inequalities we see that δ = δ + v ≤ h 1/2 0 for all sufficiently small h > 0. In addition, we find The result hence follows from the contraction mapping theorem.
Proof of Theorem 2.1. We write (c h , v h ) for the unique solution to the fixed point problem (6.10) that is provided by Lemma 6.3. This allows us to define For fixed h > 0, we claim that the map is continous. Indeed, this follows from the smoothness of Ψ * together with (A.6) and the fact that the translation operator is continuous on H 1 . Since the map is continuous on a subset of 2 h that contains ev ϑ v h for all ϑ ∈ [0, h], we conclude that is continuous. The travelling wave equation (6.1) now implies the inclusion (2.21).
In a similar fashion, the inclusion (2.25) follows from (A.6) and the continuity of the translation operator on H 1 . The remaining statements are a direct consequence of Lemma 6.3.
We now turn to the proof of Corollary 2.2, which asserts the existence of a waveprofile Φ h in the original physical coordinates. The key tool for our purpose here is [26,Prop. 4.2], which states that the gridpoints associated to a solution U of (2.13) satisfẏ Here the sequence M can be written as see [26, Eq. (6.31) and (6.33)] where this function was referred to as Y. Notice the strong resemblance with the structure of (3.42). Indeed, we see that see also [26,Eq. (6.9)] for comparison. In view of the identities from [26,Eq. (7.29)], (D.1), (5.1) and (D.11), it makes sense to formally factor out ∂ 0 U and introduce the approximant This allows us to extract a crucial lower bound for the speed of the gridpoints.
Proof of Corollary 2.2. Upon defining Notice that the terms appearing in (6.39) also all appear in (3.42) after making the replacement Z − → Y 1 . The reduction Z − → Z − apx leads to error terms that are covered by the theory in [26]. Since Z − apx does not need to be reduced further, we can follow all the computations in the present paper to obtain the error bound which is the natural analogue of (5.5). Substituting U = Ψ h and applying the Lipschitz bounds (C.6), we find In a similar fashion, we may exploit (B.4) to conclude Together, these observations yield the pointwise bound Assuming for clarity that c * > 0, this implies the pointwise inequality Since γ Ψ * is strictly bounded away from zero, uniformly in h, we conclude that for all sufficiently small h > 0 and all τ ∈ R. This shows that the coordinate transformation (2.32) is invertible, as desired.

A Discrete calculus
In this appendix we collect several useful identities and bounds from [26,27] related to the interplay between discrete and continuous calculus. In particular, we state a discrete version of the product rule, provide two summation-by-parts identities and show how taking discrete samples of functions in L 2 and H 1 affects the various norms.
Recalling the notation introduced at the start of §3.1 and §3.2, a short computation yields the basic identities together with the product rules which hold for a, b ∈ ∞ h . As in [26, §3.1], these can subsequently be used to derive the second-order product rule Recalling the discrete summation operators (2.10), one can read-off the identities for a ∈ 1 (hZ; R). In addition, the discrete summation-by-parts identities hold whenever a, b ∈ 2 h ; see [26,Eq. (3.13) and (3.15)]. Turning to sampling issues, we repeat the useful estimates [26, Eq. (A.6), (A.4)] which state that for any u ∈ H 1 . On the other hand, for any q ∈ {2, ∞} and u ∈ W 1,q , we have for any h > 0; see [26,Eq. (A.3), (A. 13)]. For any q ∈ {2, ∞} and h > 0 we also obtain the error estimate whenever u ∈ W 2;q (see [27,Lem. 4.1]), together with max{ ∂ whenever u ∈ W 3;q (see [27,Cor. 4.2]) and finally whenever u ∈ W 4;q (see [27,Cor. 4.3]). We now recall the sampling operator ev ϑ defined in (3.17). Our final two results here show how to pass back and forth between discrete and continuous estimates.
Corollary A.1 ([26, Cor. A.3]). There exists K > 0 so that for any ϑ ∈ R, any v ∈ H 1 and any 0 < h < 1, we have the bounds . Consider any f ∈ C(R; R) and any g ∈ H 1 . Then the following properties hold for all h > 0.
The gridspace function γ U The gridpoint spacing function plays an important role throughout this paper and was analyzed at length in the prequels [26,27]. We recall some of these results here and also obtain several novel bounds related to the sums that are evaluated in §E and §F. Recalling the definitions (3.24) for sums and products, we first state some useful identities pertaining to powers of γ U .
Turning to estimates, we first note that holds for any U (a) , U (b) ∈ Ω h;κ ; see [26,Eq. (C.4)]. This can be used [26,Cor. D.2] to obtain the Lipschitz bound for q ∈ {2, ∞}, where K depends on κ but not on h. In addition, it can be exploited to establish the following approximation errors for various expressions involving γ U .
We now continue the discussion from [26, §D] and consider discrete versions of the integral identities (B.10) Instead of computing the corresponding sums exactly, we obtain useful approximation that are O(h)-accurate.
Lemma B.4. Fix 0 < κ < 1 12 . Then there exists K > 0 so that for any h > 0 and any U ∈ Ω h;κ , the two linear expressions satisfy the pointwise estimate for all V ∈ 2 h . Proof. Using (A.4) we first observe that The summation-by-parts identity (A.5) allows us to compute (B.14) Upon writing we use the identity together with (B.6) to obtain We now write In particular, (B.6) yields We now transfer the S + using the summation-by-parts identity (A.5) to obtain We hence see that Using the fact that the desired estimate follows.
Lemma B.5. Fix 0 < κ < 1 12 . Then there exists K > 0 so that for any h > 0 and any U ∈ Ω h;κ , we have the pointwise estimate Proof. Since [γ U ] jh → 1 as j → −∞, we have In particular, writing we may use the estimates (B.5) to obtain Using the second summation-by-parts identity in (A.5), we can transfer the S + to obtain In particular, writing we see that from which the desired estimate follows.

C Operator bounds
Our goal here is to establish several crucial bounds on the linear operators and error functions introduced in §3.1-3.2. The errors that arise when approximating ∂ + M U and ∂ + E tw by M + U ;apx and E + tw;apx are of special importance.
Proposition C.1 ([27, Prop. 5.1]). Assume that (Hg) is satisfied and fix κ > 0. There exists K > 0 so that for any h > 0, U ∈ Ω h;κ and V ∈ 2 h we have the a-priori bounds together with the estimate In addition, for any h > 0, any pair (U (1) , U (2) ) ∈ Ω 2 h;κ and any V ∈ 2 h , we have the Lipschitz bound Proposition C.2. Assume that (Hg) and (HΦ * ) are satisfied and fix 0 < κ < 1 12 . There exists K > 0 so that for any h > 0 and U ∈ Ω h;κ we have the a-priori bounds together with the estimate while for any U (1) ∈ Ω h;κ and U (2) ∈ Ω h;κ we have the Lipschitz bounds In order to establish (C.5), we compute and notice that Upon estimating we can use (B.5) together with (B.6) to obtain the desired bound.
Proposition C.3. Assume that (Hg) and (HΦ * ) are satisfied. Then there exists K > 0 so that for any h > 0 we have the estimates Proof. We have Ψ * ∈ W 3;q for q ∈ {2, ∞}, which allows us to apply (A.7) and (A.6) to obtain This yields the first bound.
Using the fact that Ψ * ∈ W 4,2 ∩ W 4,∞ , which allows us to apply (A.10), we may argue in a similar fashion as above to conclude for q ∈ {2, ∞}. The differentiated travelling wave equation (3.33) allows us to write (C.20) Using (C.19) together with (C.5) we may hence conclude which yields the third bound.
D Decompositions for f ∈ S nl;short Our goal here is to provide the explicit decompositions (4.1) for the nonlinearities (3.35) and (3.36).
In addition, we validate the bookkeeping claims made in Propositions 4.1-4.6, providing the underpinning for the estimates in §4.2. For efficiency purposes, we combine our treatment of nonlinearities that admit similar bounds.

D.1 Decompositions for Y 1 and X A
Recalling the definitions from (3.37) and (3.39), we realize the splittings (4.1) by writing Proof. The bounds follow from inspection.
Lemma D.2. Fix 0 < κ < 1 12 and pick f ∈ {Y 1 , X A }. There exists a constant K > 0 so that the following properties are true.
Lemma D.3. Fix 0 < κ < 1 12 and pick f ∈ {Y 1 , X A }. There exists a constant K > 0 so that the following properties are true.
(i) For any h > 0, any pair U ∈ Ω h;κ and any V ∈ 2 h , we have the bound (ii) For any h > 0, any pair (U (1) , U (2) ) ∈ Ω 2 h;κ and any V ∈ 2 h , we have the bounds Proof. Recalling the Lipschitz bound (C.6), the estimates follow by inspection.
Proof. This follows from Proposition C.1.
Lemma D.5. Fix 0 < κ < 1 12 . There exists a constant K > 0 so that the following properties are true.

D.4
Decompositions for X B , X C and X D 23) from (3.38), we realize the splittings (4.1) by writing

Recalling the definitions
together with for the explicit terms. The shift terms are given by while the remainder terms are given by X B;apx;rem (U ) = X C;apx;rem (U ) = X D;apx;rem (U ) = 0, (D. 28) together with hold for all h > 0, U ∈ Ω h;κ and V ∈ 2 h . Proof. These bounds follow by inspection.
Lemma D.9. Fix 0 < κ < 1 12 and pick f ∈ {X B , X C , X D }. There exists a constant K > 0 so that the following properties are true.

(D.35)
Proof. These bounds follow from the discrete derivative expressions in Lemma B.1 and the Lipschitz bounds for γ U in (B.4).
Lemma D.10. Fix 0 < κ < 1 12 and pick f ∈ {X B , X C , X D }. There exists a constant K > 0 so that the following properties are true.

D.5 Decomposition for Y + 1
Recalling the definitions (D. 38) from (3.37) and (3.39), we realize the splittings (4.1) by writing Notice that we have eliminated the T + [∂ (2) U ] term in the explicit expressions, while keeping the T + [∂ (2) V ] dependency. This inconsistency is deliberate as it will help us to make a useful substitution in the sequel. In addition, we introduce the sets h . Proof. These bounds follow by inspection.
Lemma D.12. Fix 0 < κ < 1 12 . There exists a constant K > 0 so that the following properties are true.
(i) For any h > 0 and U ∈ Ω h;κ we have the bound for any h > 0, any pair U ∈ Ω h;κ and any V ∈ 2 h . Proof. These estimates follow by inspection.

D.6 Decomposition for Y + 2b
Recalling the definitions from (3.37) and (3.39), we realize the first splitting in (4.1) by writing The second splitting (4.1) is obtained implicitly by recalling the definition (3.15) and writing for any h > 0, any pair U ∈ Ω h;κ and any V ∈ 2 h . Proof. Recalling (3.14), we make the decomposition (D.56) Introducing the function

E Reductions for G
Our goal here is to construct the functions G apx;sh;b , G apx;rem;b , G lin;U ;sh;b and G lin;U ;rem;b and demonstrate that they satisfy the corresponding bounds in Propositions 5.2-5.3. We proceed in a relatively direct fashion, treating each of the components in (3.42) separately and subsequently combining the results.

E.1 Simplifications for G A
We recall the definition Substituting the relevant expressions from §D we find We now make the decomposition We also recall the definition Substituting the relevant expressions from §D, we find We now make the decomposition We summarize our results by writing and obtaining the following bound.
Lemma E.1. Assume that (Hg) is satisfied and pick 0 < κ < 1 12 . Then there exists a constant K > 0 so that the following properties hold true. Proof. In view of Proposition C.1, the bounds follow by inspection.

E.2 Simplifications for G B
We recall the definition Substituting the relevant expressions from §D we find

E.4 Final decomposition
Recalling the definitions (5.1), we observe that (E.30) Proof of Propositions 5. F Reductions for G + Our goal here is to construct the functions G + apx;sh;b , G + apx;rem;b , G + lin;U ;sh;b and G + lin;U ;rem;b and demonstrate that they satisfy the corresponding bounds in Propositions 5.2-5.3. As in the previous section, we treat each of the components in (3.46) and (3.47) separately and subsequently combine the results.

A b
We recall the definition G + A b;apx;II (U ) = 1 − Y 1;apx;expl (U )X A;apx;expl (U ) Y + 2b;apx;expl (U ). (F.1) Substituting the relevant expressions from §D, we find We now make the decomposition We conclude by writing and obtaining the following bound.
Lemma F.1. Assume that (Hg) is satisfied and pick 0 < κ < 1 12 . Then there exists a constant K > 0 so that the following properties hold true.  and obtaining the following bounds.
Lemma F.2. Assume that (Hg) is satisfied and pick 0 < κ < 1 12 . Then there exists a constant K > 0 so that the following properties hold true.   and obtaining the following bounds.
Lemma F.3. Assume that (Hg) is satisfied and pick 0 < κ < 1 12 . Then there exists a constant K > 0 so that the following properties hold true.

F.4 Final decomposition
Arguing as in §E.3 we see that