Travelling Corners for Spatially Discrete Reaction-Diffusion System

We consider reaction-diffusion equations on the planar square lattice that admit spectrally stable planar travelling wave solutions. We show that these solutions can be continued into a branch of travelling corners. As an example, we consider the monochromatic and bichromatic Nagumo lattice differential equation and show that both systems exhibit interior and exterior corners. Our result is valid in the setting where the group velocity is zero. In this case, the equations for the corner can be written as a difference equation posed on an appropriate Hilbert space. Using a non-standard global center manifold reduction, we recover a two-component difference equation that describes the behaviour of solutions that bifurcate off the planar travelling wave. The main technical complication is the lack of regularity caused by the spatial discreteness, which prevents the symmetry group from being factored out in a standard fashion.


Introduction
In this paper we construct travelling corner solutions to a class of planar lattice differential equations (LDEs) that includes the Nagumo LDĖ u i,j = u i+1,j + u i−1,j + u i,j+1 + u i,j−1 − 4u i,j + g cub (u; ρ) (1.1) posed on the two-dimensional square lattice (i, j) ∈ Z 2 , in which the nonlinearity is given by the bistable cubic g cub (u; ρ) = (u 2 − 1)(ρ − u), −1 < ρ < 1. (1.2) Such corners can be seen as interfaces that connect planar waves travelling in slightly different directions. In particular, our analysis does not require the use of the comparison principle, but merely requires a number of spectral and geometric conditions to hold for the underlying planar travelling waves. This allows our results to be applied to a wide range of LDEs, highlighting the important role that anisotropy and topology play in spatially discrete settings.
Reaction-diffusion systems The LDE (1.1) can be seen as a nearest-neighbour spatial discretization of the Nagumo PDE u t = u xx + u yy + g cub (u; ρ).
(1. 3) In modelling contexts one often uses the two stable equilibria of the nonlinearity g to represent material phases or biological species that compete for dominance in a spatial domain. Indeed, the diffusion term tends to attenuate high frequency oscillations, while the bistable nonlinearity promotes these. The balance between these two dynamical features leads to interesting pattern forming behaviour. As a consequence, the PDE (1.3) has served as a prototype system for the understanding of many basic concepts at the heart of dynamical systems theory, including the existence and stability of planar travelling waves, the expansion of localized structures and the study of obstacles. Multicomponent versions of (1.3) such as the Gray-Scott model [19] play an important role in the formation of patterns, generating spatially periodic structures from equilibria that destabilize through Turing bifurcations. Memory devices have been designed using FitzHugh-Nagumo-type systems with two components [31], which support stable stationary radially symmetric spot patterns. Similarly, one can find stable travelling spots [46] for three-component FitzHugh-Nagumo systems, which have been used to describe gas discharges [38,42].
At present, a major effort is underway to understand the impact that non-local effects can have on reaction-diffusion systems. For example, many neural field models include infinite-range convolution terms to describe the dynamics of large networks of neurons [10,11,39,43], which interact with each other over long distances. The description of phase transitions in Ising models [3,4] features non-local interactions that can be both attractive and repulsive depending on the length scale involved.
It is well-known by now that the topology of the underlying spatial domain can have a major impact on the dynamical behaviour exhibited by such non-local systems. For example, nerve fibers have a myeline coating that admits gaps at regular intervals [40], which can block signals from propagating through the fiber [15,30,34]. In order to study the growth of plants, one must take into account that cells divide and grow in a fashion that is influenced heavily by the spatial configuration of their neighbours [20]. Finally, the periodic structure inherent in many meta-materials strongly influences the phase transitions that can occur [13,14,44] as a consequence of the visco-elastic interactions between their building blocks.
We view the planar LDE (1.1) as a prototype model that allows the impact of such non-local spatially-discrete effects to be explored. Indeed, the spatial R 2 → Z 2 transition breaks the locality but also the translational and rotational symmetry of (1.3), leading to several interesting phenomena and mathematical challenges.
Existence of planar waves It is well-known that the balance between the diffusion and reaction terms in the PDE (1.3) is resolved through the formation of planar travelling wave solutions u(x, y, t) = Φ(x cos ζ + y sin ζ + ct); Φ(−∞) = −1, Φ(∞) = 1, (1.4) which connect the two stable equilibria u = ±1. When c = 0, these waves can be thought of as a mechanism by which the fitter species or more energetically favourable phase invades the spatial domain. The existence of these waves can be obtained by applying a phase-plane analysis [18] to the travelling wave ODE cΦ = Φ + g cub (Φ; ρ), Φ(±∞) = ±1, (1.5) which results after substituting (1.4) into (1.3).
(1.7) The broken rotational invariance in the transition from (1.3) to (1.1) is manifested by the explicit presence of the propagation direction in (1.7). The broken translational invariance causes the wavespeed c to appear in (1.7) as a singular parameter.
A comprehensive existence theory for solutions to (1.7) was obtained in [36]. In particular, for every ζ ∈ [0, 2π] and ρ ∈ [−1, 1] there exists a unique wavespeed c = c ρ,ζ for which (1.7) admits a solution. However, it is a delicate question to decide whether c = 0 or c = 0. Indeed, a sufficient energy difference between the two stable equilibrium states is needed for the propagation of waves [4,6,17,33], a phenomenon referred to as propagation failure. In fact, due to the angular dependence in (1.7), planar waves can fail to propagate in certain directions that resonate with the lattice, whilst travelling freely in others [12,25,37].
Linearization It is well-known that planar travelling waves can be used as a skeleton to describe the global dynamics of the PDE (1.4) [1]. In particular, they have been used as building blocks to construct other more complicated types of solutions. A key ingredient in such constructions is to understand the dynamics of the system that arises after linearizing (1.3) around the planar waves (1.4).
Performing this linearization for ζ = 0, we obtain the system ∂ t v(x, y, t) = ∂ xx v(x, y, t) + ∂ yy v(x, y, t) + g cub Φ(x + ct); ρ v(x, y, t), (1.8) which can be transformed to the temporally autonomous system ∂ t v(x, y, t) = ∂ xx v(x, y, t) + ∂ yy v(x, y, t) − c∂ x v(x, y, t) + g cub Φ(x); ρ v(x, y, t) (1.9) by the variable transformation x = x + ct. Since this system is also autonomous with respect to the y-coordinate, which is transverse to the motion of the wave, it is convenient to apply a Fourier transform in this direction. Upon introducing the symbol L z p (x) = ∂ xx p(x) + z 2 p(x) − c∂ x p(x) + g cub Φ(x); ρ p(x), (1. 10) we readily find ∂ tvω (x, t) = L iωvω (·, t)](x). (1.11) Inspecting (1.10), we readily see that the spectrum of L z can be obtained by rigidly shifting the spectrum of L 0 by z 2 . In particular, writing λ z = z 2 , we find that (1.12) Noting that λ iω = −ω 2 , we hence see that perturbations of the form v(x, y, 0) = θ(y, 0)Φ (x) evolve under (1.9) according to the heat semiflow θ t = θ yy . These perturbations are important because they correspond at the linear level with transverse deformations of the planar wave interface.
(1.17) The theory developed in [5,7] essentially justifies this formal calculation and confirms that the spectral properties of L z can be used to understand the dynamics of the time-dependent problem (1.14). Upon writing λ z = 2(cosh(z) − 1), we again have To find the evolution of perturbations of the form v ij (0) = θ j Φ under (1.14), we must now solve the discrete heat equationθ The situation is hence similar to that encountered for the PDE (1.3). More material changes arise however when considering the diagonal direction ζ = π 4 . Following a similar procedure as above, one arrives at the linear operator [L diag z p](ξ) = −cp (ξ) + 2 cosh(z) p(ξ + 1) + p(ξ − 1) − 4p(ξ) + g cub Φ * (ξ); ρ p(ξ), (1.20) which has a spectrum that can no longer be directly related to that of L diag 0 . It is hence no longer clear how to formulate an analogue for (1.19) to describe the linear evolution of interface deformations. However, it is still the case that λ z = O(z 2 ) as z → 0 for the curve of eigenvalues that bifurcates from the zero eigenvalue L diag 0 Φ = 0. For general rational angles ζ this quadratic behaviour need no longer be true. In fact, we obtain the relation for the quantity that is often referred to as the group velocity. A similar relation was found in [21] for planar PDEs with direction-dependent diffusion coefficients. However, in this case it is always possible to change the coordinate system in such a way that λ z = O(z 2 ) holds again. Such a transformation is not possible in the spatially discrete setting (1.1), since this would require the transverse spatial coordinate to become continuous. However, we do remark here that the function φ → c ρ,ζ can behave rather wildly in the critical regime where ρ is small, allowing the group velocity to vanish at specific values for ρ even if ζ / ∈ π 4 Z.
Stability of planar waves The realization that transverse interfacial deformations are governed by a heat equation led to the development of two main approaches to establish the nonlinear stability of the planar waves (1.4). Both approaches exploit the coordinate system u(x, y, t) = Φ x + ct + θ(y, t) + v(x, y, t) (1.22) in the neighbourhood of the planar travelling wave and require the initial perturbations θ(y, 0) and v(x, y, 0) to be localized in a suitable sense.
The first approach was pioneered by Kapitula in [32], where he used semigroup methods and fixed-point arguments to show that θ tends algebraically to zero, while v decays exponentially fast. The advantage of this approach is that only weak spectral assumptions need to be imposed on the underlying system. However, the crude estimates on the nonlinear terms lead to rather weak estimates for the basin of attraction.
The second approach leverages the comparison principle to obtain stability for a much larger class of initial perturbations. By slowing down the natural decay-rate of the fundamental solution of the heat equation, the authors of the landmark paper [8] were able to construct explicit super and sub-solutions to (1.3) that trap perturbations that can be arbitrarily large (but localized). In fact, the authors use their construction to show that these planar waves can pass around large compact obstacles and still eventually recover their shape.
In [23,24] these approaches were generalized to the discrete setting of (1.1), thereby continuing the early work by Bates and Chen [2] featuring a related four-dimensional non-local problem. In both cases the key technical challenge was the analysis of troublesome non-selfadjoint terms spawned by the anisotropy of the lattice, especially in situation where the group velocity does not vanish. These terms have slower decay rates than their PDE counterparts and hence require special care to close the nonlinear bootstrapping procedure. For example, the sub-solutions in [8] consist of only two terms, while 33 terms were required in [23] to correct for the slower decay.
Spreading phenomena The classic result [1,Thm. 5.3] obtained by Weinberger for the PDE (1.3) states that large compact blobs with u ≈ 1 inside and u ≈ −1 outside can expand throughout the plane. The proof of this result relies on the construction of radially expanding sub-and supersolutions by glueing together planar travelling waves.
In [23] a weak version of this expansion result was established for the LDE (1.1) in the special case that no direction is pinned. However, the underlying sub-and super-solutions expand at the speeds min 0≤ζ≤2π c ρ,ζ and max 0≤ζ≤2π c ρ,ζ respectively, which still leaves a considerable hole in our knowledge of the expansion process. Indeed, the numerical results in [45] provide strong evidence that the limiting shape of the expanding blob can be found by applying the Wulff construction [41] to the polar plot of the ζ → c ρ,ζ relation. For a large subset of parameters ρ this limiting shape resembles a polygon.
The main motivation behind the current paper is to take a step towards understanding this expansion process by looking at the evolution of a single corner. Indeed, when the expanding blob is sufficiently large, it would seem to be very reasonable to assume that the corners of the polygon behave in an almost independent fashion.
Corners for PDEs Assuming for concreteness that ρ < 0, the horizontal planar wave (c, Φ) given by (1.4) with φ = 0 satisfies c > 0, which means that it travels towards the left. In [22] Haragus and Scheel construct travelling corner solutions to (1.3) by 'bending' this planar wave to the left in the spatial limits y → ±∞, so that the interface resembles a > sign.
In particular, for any small opening angle ϕ > 0, the authors establish the existence and stability of solutions of the form Here v(·, y) H 2 = O(ϕ 2 ) uniformly in y, while the phase θ satisfies the limits lim y→±∞ θ (y) = ± tan ϕ. (1.24) Notice that the horizontal speed c cos ϕ of these corners is faster than the original speed of the planar wave.
The result is obtained by using the change of variable x = x +ct to recast (1.3) as and subsequently demanding u t = 0. The resulting system can be written in the first-order form which admits the family of y-independent equilibriã The linearization of (1.26) around (c, u) = (c, Φ) can be written as (1.28) This system admits the y-independent solutions (Φ , 0) caused by the translational invariance, together with the linearly growing solution (yΦ , Φ ). In particular, the desired corner (1.23) lives on the two-dimensional global center manifold associated to the family (1.27). The solutions on this manifold can be represented in the form One can subsequently obtain two skew-coupled ODEs to describe the dynamics of the scalar functions κ and θ. A relatively straightforward analysis shows that these ODEs have solutions for which θ satisfies the limits (1.24), while κ remains small. This suffices to establish the existence of the corners (1.23).
In addition, in [21] anisotropic effects were introduced into the problem by allowing the nonlinearity g to depend on the gradient of u and considering non-diagonal diffusion coefficients. In such cases the group velocity c g defined by the quantities (1.21) need not vanish, but it can be removed by applying a coordinate transformation y = y − c g t in the transverse direction.
By restricting their attention to small opening angles ϕ and using center manifold arguments, Scheel and Haragus were also able to apply their techniques to multi-component reaction-diffusion PDEs such as the FitzHugh-Nagumo and Gray-Scott equations [21,22] However, it is also possible to consider large opening angles when considering equations that admit a comparison principle. Indeed, in [9] explicit sub-and super-solutions are used to construct corners for the Nagumo PDE (1.3) that can be arbitrarily sharp.

Corners for LDEs
The crucial point in the analysis outlined above for the corners (1.23) is that the phase shift θ(y) can be completely factored out from the system. This implies that the ODE for κ does not depend on θ. In addition, it allows the center manifold to be constructed by a standard fixed point argument analogous to the local case. This is possible because the right-hand side of (1.28) maps H 2 × H 1 into H 1 × L 2 , which roughly means that its inverse gains an order of regularity in both components. This precisely compensates for the loss of regularity that arises by factoring out the phase-shift.
However, when attempting to mimic this procedure for the LDE (1.1) one runs into a fundamental difficulty. Indeed, the analog of (1.28) has a right-hand side that now maps H 1 ×H 1 into H 1 ×L 2 due to the lack of second derivatives in the equation. This forces us to construct a full two-dimensional global center manifold that takes into the account the dynamics of θ and v simultaneously.
A similar situation was encountered by one of the authors in [27], where modulated travelling wave solutions were constructed to a class of non-local systems. However, the analog of (1.28) is a difference equation rather than a differential equation. The order of this difference equation can become arbitrarily large depending on the height of the fraction tan ζ, which we require to be rational. Nevertheless, the final step in our analysis requires us to uncover a first-order difference equation for the center variables.
The main technical contribution in this paper is that we adapt the spirit of the approach in [27] to construct global center manifolds for the differential-difference systems that we encounter here. This approach uses two intertwined fixed point procedures to separate the flow problem for the two center variables from the task of capturing the shape of the remainder function h * . The underlying linear problems have non-autonomous slowly-varying coefficients, for which we develop appropriate solution operators.
In this paper we do require the group velocity (1.21) to vanish. Unlike in the spatially continuous setting, this cannot always be arranged by a simple variable transformation. Indeed, such a transformation would force the spatial variable transverse to the propagation direction to become continuous, destroying the difference structure of the system. This prevents us from exploiting the ω → ω + 2π periodicity in the Fourier variable. As a result, resonances start to appear in the spectrum that are very hard to control. A similar situation was encountered in [27], which forced the authors to add a smoothening term to the underlying system.
We emphasize that the group velocity for (1.1) vanishes automatically in the directions ζ = 0 and ζ = π 4 . In addition, directions where the wavespeed is minimal (and hence the group velocity is zero) play an important role in the Wulff construction, which is the primary motivation for our analysis here. In any case, the delicate behaviour of the c ρ,ζ map for the Nagumo LDE (1.1) leads to a much richer class of behaviour than that displayed by its continuous counterpart (1.3). For example, the latter only features interior corners, while the former can also admit exterior corners. The former also allows for so-called bichromatic corners, which connect spatially homogeneous equilibria to checkerboard patterns.
While we are confident that our center manifold construction will also allow us to establish the (linearized) stability of the corners constructed here along the lines of the approach in [21], we do not pursue this in the present paper. The main reason is that there is no coordinate transformation that can freeze our corners and also leave the discrete structure of the equation intact. One would need to generalize the approach developed in [5,7,24] to accommodate solutions that vary in two directions instead of just one, which we expect to be a tedious task.
Organization Our main results are formulated in §2 and applied to the Nagumo LDE (1.1) in §2.1-2.2. In §3 we derive the differential-difference system that the pair (θ, v) must satisfy and formulate the global center manifold result. We proceed in §4 by deriving a representation formula for solutions to the linearized problem with constant phase. This requires us to compute a convoluted spectral projection operator that arises from the second order pole that the operator L −1 z has in z = 0. In §5- §6 we combine this representation formula with Fourier analysis to construct a solution operator for the linearized problem where the phase is allowed to vary slowly. Finally, we setup the fixed point problems required to build the global center manifold in §7, appealing at times to the results in [27] for overlapping parts of the program. posed on the planar lattice (i, j) ∈ Z 2 , in which u takes values in R d . For convenience, we introduce the operator π + ij : ∞ (Z 2 ; R d ) → (R d ) 5 that acts as for any (i, j) ∈ Z 2 , which allows us to rewrite (2.1) in the condensed forṁ 3) The plus sign corresponds with the fact that a "+"-shaped stencil is used to sample u. The conditions we impose on the nonlinearity f are summarized in the following assumption.
(Hf) The nonlinearity f : (R d ) 5 → R d is C r -smooth for some r ≥ 2 and there exist two points We emphasize that the two points u ± are allowed to be equal. These two equilibria are required to be connected by a planar travelling wave solution to (2.1). In particular, we pick an arbitrary rational direction (σ A , σ B ) ∈ Z 2 with gcd(σ A , σ B ) = 1 and impose the following condition.
(HΦ) There exists a wave speed c * = 0 and a wave profile Φ * ∈ C r+1 (R, R d ) so that the function In addition, we have the limits Upon introducing the operator τ : C(R; R d ) → C(R; (R d ) 5 ) that acts as we note that the pair (c * , Φ * ) must satisfy the functional differential equation of mixed type (MFDE) (2.8) In particular, the C r+1 -continuity mentioned in (HΦ) is automatic upon assuming that Φ * is merely continuous. For convenience, we now introduce the new coordinates which are parallel respectively orthogonal to the direction of motion of the wave (2.5). Upon introducing the notation which admits the travelling wave solution u nl (t) = Φ * (n + c * t). (2.12) A standard approach towards establishing the stability of the wave (2.12) under the nonlinear dynamics of the LDE (2.11) is to consider the linear variational probleṁ Looking for a solution of the form v nl (t) = e λt e zl p(n + c * t), (2.14) we readily find that p must satisfy the eigenvalue problem Here the linear operator L z : with shifts r j and functions A z,j that are given by (2.17) Since Φ * (ξ) approaches u ± as ξ → ±∞, it is possible to define the characteristic C d×d -valued functions Our first spectral assumption states that these characteristic functions cannot have roots on the imaginary axis whenever z is purely imaginary.
(HS2) For any ω = 0 the operator L iω is invertible as a map from Since the Fredholm index varies continuously, (HS1) and (HS2) together imply that the Fredholm index of L 0 is zero. The translational invariance of the problem implies that L 0 Φ * = 0, which means that zero is an eigenvalue for L 0 . Our next assumption states that this eigenvalue is in fact algebraically simple. For any z ∈ C we now introduce the linear operator that acts as An easy computation shows that holds for all pairs p, q ∈ W 1,∞ (R, C d ). For these reason, we refer to this operator as the formal adjoint of L z . Using [35,Thm. A] together with (HS3), one sees that the kernel of L adj 0 must also be onedimensional. In particular, it is spanned by a function ψ * ∈ W 1,∞ (R, R d ) that can be uniquely fixed by the identity on account of (2.21). We note that (HS1) implies that both Φ * (ξ) and ψ * (ξ) decay exponentially as ξ → ±∞. We now explore two important consequences of the algebraic simplicity condition (HS3). The first of these states that the zero eigenvalue can be extended to a branch of eigenvalues λ z for L z when |z| is small.  (2.26) defined for each z ∈ C with |z| < δ z , such that the following hold true.
(i) The characterization together with the algebraic simplicity condition hold for each z ∈ C with |z| < δ z .
The second consequence is that wave (c * , Φ * ) travelling in the rational direction (σ A , σ B ) can be perturbed to yield waves travelling in nearby directions. In particular, we introduce the constants (ζ * , σ * ) by writing Looking for solutions to the LDE (2.11) of the form a short computation shows that the pair (c ϕ , Φ ϕ ) must satisfy the MFDE in which we have introduced the notation (2.33) In order to translate these waves back to the original coordinates, we remark that any solution to (2.32) yields a solution to the original LDE (2.1) by writing with the rescaled quantities . Assume that (Hf ), (HΦ) and (HS1)-(HS3) are all satisfied. Then there exists a constant δ ϕ > 0 together with pairs defined for each ϕ ∈ (−δ ϕ , δ ϕ ), such that the following hold true.
satisfies the LDE (2.11) for all t ∈ R.
(iv) We have the identities We remark that the first quantities in (2.39) can be interpreted as a so-called group velocity, which represents the speed at which long-amplitude perturbations travel in the transverse direction. Indeed, expanding (2.14) with p = φ z and λ = λ z we find (2.40) Our final condition requires λ z to depend quadratically on z, which means that this group velocity has to vanish. We emphasize that the inequality [∂ 2 z λ z ] z=0 > 0 was required in [24] to obtain the nonlinear stability of the planar wave (c * , Φ * ).
(HM) We have the identities As a final preparation, we introduce the directional dispersion Assuming that the original wave travels in the horizontal direction ζ * = 0, the quantity d ϕ represents the speed at which level-sets of the wave (c ϕ , Φ ϕ ) travel along the horizontal axis; see Figure 1. An easy calculation using (2.41) shows that Our main result establishes the existence of travelling corners in the setting where [∂ 2 ϕ d ϕ ] ϕ=0 = 0. Assuming again that φ * = 0 and that [∂ 2 z λ z ] z=0 and [∂ 2 ϕ d ϕ ] ϕ=0 are both strictly positive, the levelsets resemble a > sign. In particular, when c * > 0 this resembles an interior corner travelling to the left.
together with two angles ϕ − < 0 < ϕ + that satisfy the following properties.
satisfies the LDE (2.11) for all t ∈ R.
(iii) We have the identities have the same sign, then we have the limits On the other hand, if these quantities have opposing signs, then we have the limits

The Nagumo LDE
As an example, we return to the Nagumo LDĖ in which the nonlinearity is given by the scaled cubic for some detuning parameter ρ ∈ (−1, 1). In the terminology of (2.1), we hence have which shows that (Hf) is satisfied upon picking u ± = ±1.
Turning to (HΦ), we note that the results in [36] show that for each ζ ∈ [0, 2π] and ρ ∈ (−1, 1) there is a unique wavespeed c = c ρ,ζ for which the system admits a monotonic solution Φ = Φ ρ,ζ that also satisfies the limits (2.6). Figure 2 contains polar plots of the ζ → c ρ,ζ relation, which can be very delicate whenever |ρ| is small. By symmetry, we have c ρ,ζ = −c −ρ,ζ and hence c 0,ζ = 0 for all angles ζ ∈ [0, 2π]. Upon writing ρ * (ζ) = sup{ρ : c ρ,ζ = 0}, (2.56) the results in [36] show that 0 ≤ ρ * (ζ) < 1 for all ζ ∈ [0, 2π]. In particular, this means that (HΦ) is satisfied whenever tan ζ is rational (or infinite) and ρ * (ζ) < |ρ| < 1. Under the same conditions, the discussion in [24, §6] uses arguments based on the comparison principle to show that also (HS1)-(HS3) are valid. The verification of the conditions in (HM) is much more subtle. In order to make the angular dependence fully explicit, we first pick and consider the operatorsL that act as Writingλ z for the branch of eigenvalues forL z bifurcating fromλ 0 = 0 and comparing this to the branch λ z defined in Lemma 2.1, we haveλ with σ * > 0 the smallest number so that In particular, we have In view of the similar rescalings (2.35) and the fact that the statements in Theorem 2.3 merely depend on the signs of the quantities [∂ 2 and focus on the eigenvalues λ z;ρ,ζ and eigenfunctions φ z;ρ,ζ bifurcating from (0, Φ ρ,ζ ) for ρ * (ζ) < |ρ| < 1. We write ψ ρ,ζ for the solution to the adjoint equation In our context, the operators defined in (4.13) and (4.13) act as (2.65) In particular, Lemma 4.2 allows us to compute which in turn allows us to find [∂ z φ z;ρ,ζ ] z=0 by solving the MFDE In addition, item (iv) of Lemma 2.2 shows that Turning to the second derivatives, we again use Lemma 4.2 to compute (2.69) We remark that the last line vanishes in principle if the normalization (2.29) is imposed. However, numerically it is convenient to be free to utilize a different normalization, in which case this term should be included.
Notice the extra minima that start to form in the directions tan ζ = 1 and subsequently tan ζ = 2 3 as ρ is decreased.
Finally, in view of the fact that that D 2 f acts only on its fifth argument, we can use Lemma 4.5 to obtain (2.70) The last line can be ignored if indeed [∂ ζ c ρ,ζ ] = 0. In the special cases where ζ = k π 4 for some k ∈ Z, we have A 1 = 0 and hence For ζ = 0 this allows us to write On the other hand, for ζ = π 4 we have Since ψ ρ,ζ and Φ ρ,ζ are strictly positive we hence see that . The sharp spikes occur at the critical value ρ * (ζ) where pinning sets in. We note that sign changes appear for ζ = π 2 but not for ζ = 0. In particular, the identity c g ≡ 0 for these directions implies that interior and exterior corners can both occur for ζ = π 2 , while the horizontal direction ζ = 0 features interior corners only. The right panel contains numerically computed values for c g (ρ). Notice the zero-crossings for tan ζ = 3 4 and tan ζ = 4 5 , which indicates the presence of interior corners at these two critical values for ρ.
for all k ∈ Z. The numerical results in [24, §6] suggest that this inequality extends to a wide range of (ρ, ζ) and we take this for granted for the remainder of our discussion. However, even for the straightforward expressions (2.72)-(2.73), it is not clear whether the quantity c ρ,ζ + ∂ 2 ζ c ρ,ζ has a sign.
For any fixed ζ * , we now introduce the notation for the group velocity and second derivative of the directional dispersion that play a role in Theorem 2.3. In particular, to apply this result we need c g (ρ) = 0 and κ d (ρ) = 0. Since c ρ,ζ ≤ 0 whenever ρ ≥ 0, we have an interior corner for κ d (ρ) < 0 and an exterior corner for κ d (ρ) > 0. In both cases the corner travels in the rightward direction (provided |ζ| < π 2 ). In Figure 3 we have numerically computed the quantities (2.75) for a range of rational directions. In all cases, we also confirmed numerically that [∂ 2 z λ z,ρ,ζ ] z=0 > 0. In particular, the results predict interior corners travelling in the horizontal direction ζ * = 0, while both types of corners can travel in the diagonal direction ζ * = π 4 . In addition, for two critical values of ρ > 0 there are interior corners that travel in the direction ζ * = arctan(3/4) respectively ζ * = arctan(4/5). To obtain these results, we simultaneously solved the systems (2.55), (2.64) and (2.67). For well-posedness reasons, we added the extra terms γΦ , γψ respectively γ[∂ z φ z;ρ,ζ ] to the right-hand side of each equation, taking γ = 10 −6 . For the precise procedure, we refer to [24, §6].

Bichromatic Nagumo LDE
We here reconsider the Nagumo LDĖ  cos(ζ−ζ * ) , with ζ * = 0 for the left column and ζ * = π 4 for the right column, again with ρ = 0. These results strongly suggest that [∂ 2 ζ d(ζ)] ζ=ζ * can take both signs as the diffusion coefficient α is varied. In particular, both the horizontal and diagonal directions can have interior and exterior corners.
but are now interested in so-called bichromatic planar travelling wave solutions. Such solutions can be written in the form for some wavespeed c ∈ R and R 2 -valued waveprofile These waves fit into the framework of this paper, since they can be seen as travelling wave solutions for the 'doubled' LDĖ We now introduce the notation The results in [26] show that there exists an open set of values (ρ, α) with ρ ∈ (−1, 1) and α > 0 for which (2.84) admits solutions (u bc , v bc ) ∈ (0, 1) 2 that are stable spatially homogeneous equilibria for (2.79). By applying the theory in [28], we hence obtain the existence of solutions to (2.83) that satisfy the limits together with similar solutions that connect (u bc , v bc ) with (1, 1). We remark that the existence theory in [28] does not prescribe whether c = 0 or c = 0. However, in the special case ζ = π 4 the travelling wave system can be written as there exists an open set of values (ρ, α) for which the wavespeed does not vanish, allowing us to verify (HΦ). By continuity in ζ, this hence also holds for nearby angles. In addition, our numerical results suggest that the diagonal direction is the first to become pinned as α is decreased; see Figure 4. This contrasts the situation encountered in the monochromatic case, where the diagonal waves satisfy the same travelling wave MFDE as the horizontal waves, but with a doubled diffusion coefficient. As a result, the monochromatic horizontal waves pin earlier than their diagonal counterparts.
Since the spectral conditions (HS1)-(HS3) can be verified with techniques similar to those used for the monochromatic case, we now turn our attention to (HM). In particular, writing (c ρ,α,ζ , Φ ρ,α,ζ ) for the solution to (2.83) that satisfies the limits (2.85), we introduce the operator . We now introduce the notation A mc 1 , A mc 2 and B mc 1 for the operators (2.65) defined for the monochromatic equation. In addition, we write A bc 1 , A bc 2 and B bc 1 for the operators defined in (4.13) and (4.13) associated to the bichromatic problem (2.83). It is not hard to verify the relations We write (λ z;ρ,α,ζ , φ z;ρ,α,ζ ) and ψ ρ,α,ζ for the analogs of the similar named expressions defined in §2.1. Since A bc 1 = 0 for ζ = 0 and ζ = π 4 , we again have For ζ = 0 this allows us to write On the other hand, for ζ = π 4 we have (2.91) Since both components of ψ ρ,α,ζ and Φ ρ,α,ζ are strictly positive, we again see that for all k ∈ Z.
As before, it is unclear if the derivatives of c have a sign. As a consequence of the increasing number of components in the MFDEs, our numerical method is at present unable to compute the desired derivatives in the same fashion as above. Instead, we simply compute the directional dispersion relation directly and determine by inspection whether it is concave or convex; see Figure  4. Interestingly enough, we find that this characterization flips at least twice as α is decreased, both for the horizontal direction ζ = 0 and the diagonal direction ζ = π 4 . In contrast to the monochromatic case, we hence see that interior and exterior bichromatic corners can both travel in the horizontal direction.

Problem setup
In this section we setup a differential-algebraic equation to describe solutions to the LDĖ that can be written in the form u nl (t) = Ξ l (n + ct) (3.2) for some sequence Ξ : The elements Ξ l will be required to lie in the orbital vicinity of the waveprofile Φ * . In particular, we formulate a global center manifold reduction that allows us to find an equivalent two component difference equation of skew-product form.
For any sequence Ξ of the form (3.3), we introduce the notation to refer to the expanded sequence In addition, for any p = (p 1 , . . . , we introduce the function τ p that is given by This allows us to write π × nl u(t) = τ Ξ l (n + ct) + τ Ξ l (n + ct) holds for all l ∈ Z and ξ ∈ R.
From now on, we often drop the explicit ξ-dependence. For instance, we simply write instead of the longer form (3.9). For technical reasons it is advantageous to recast (3.10) as a 2d-component system. To this end, we introduce the first differences In addition, we introduce the notation for the summed sequence in which we make the convention that sums where the lower index is strictly larger than the upper index are set to zero. For example, in the special case (σ A , σ B ) = (1, 0) we have These definitions allow us to observe that . (3.15) In particular, the system (3.10) can now be rewritten in the equivalent form (3. 16) In the special case c = c * , the travelling wave solution (2.12) gives rise to l-independent solutions to (3.16) in which ϑ ∈ R can be chosen arbitrarily. Here we have introduced the left-shift operator for any p ∈ C(R; R d ).
We now look for a branch of solutions to (3.16) that bifurcates from the travelling waves (3.17) for c = c * . In particular, we consider the Ansatz for three sequences (θ, v, w) : (3.21) In order to close the system, we augment (3.21) by demanding that for all l ∈ Z.
We now set out to isolate the linear and nonlinear parts of the system (3.21). For anyṽ ∈ C(R; (R d ) 5 ) and ϑ ∈ R we therefore introduce the function N f ;ϑ (ṽ) ∈ C(R; R d ) that is given by (3.23) In addition, for any phase ϑ ∈ R and difference θ d ∈ R we introduce the function Using these functions, the system (3.21) can be recast as For any v ∈ L 2 (R; R d ), we now introduce the notation Applying a difference to (3.22), we obtain Substituting the equation for v, we arrive at In order to formulate this in a more compact fashion, we write together with ev l p = p l−σ * +1 , . . . , p l+σ * −1 , p l+σ * (3.31) for any sequence p. In addition, we introduce the shorthand notation With a slight abuse of notation, we introduce the function S : R 2σ * × H 1 × H 1 → R that acts as This allows us to rewrite (3.29) in the form For any triplet (θ, v, w) : Z → R × H 1 × H 1 , we now introduce the sequences In addition, we introduce the nonlinear function that acts as (3.39) Finally, we introduce the operators and the associated projections ϑ (v, w) = (0, P ϑ w). (3.41) This allows us to represent the full problem as  and consider the pair (Ξ, Υ) defined by (3.19). Then the differential-algebraic system (3.16) and the identity are both satisfied for all l ∈ Z, if and only if (3.42) holds for all l ∈ Z.
Proof. The computations above shows that indeed (3.42) holds whenever (3.16) and (3.45) are satisfied. The converse implication can be checked by using (3.35) to compute We now proceed to obtain estimates on the nonlinear terms. These are mostly standard bounds for quadratic nonlinearities that will be used in §7 for the center manifold construction. Notice however that any dependencies on the phase θ always involve differences in θ. (3.47) In addition, for any pair (ϑ A , ϑ B ) ∈ R 2 and any pair we readily see that Since f is at least C 2 -smooth, there exists C 1 > 0 so that the pointwise bound This yields for some C 2 > 0, from which (3.47) follows. Upon writing a short computation shows that Under the assumption that (3.53) holds for both (Φ A , v A ) and (Φ B , v B ), we hence find for some C 3 > 0. The bound (3.49) now follows from the fact that for some C 4 > 0.
Lemma 3.3. Assume that (Hf ), (HΦ) and (HS1) hold and pick a sufficiently large constant K > 0. Then for any pair ϑ ∈ R and θ d ∈ R we have the bound In addition, for any pair (ϑ A , ϑ B ) ∈ R 2 and any pair (θ A d , θ B d ) ∈ R 2 we have the bound Proof. Assumptions (Hf), (HΦ) and (HS1) imply that Φ * is a continuous function that decays exponentially, from which (3.59) follows. Writing (3.62) The inequality (3.60) follows directly from this representation.
for all l ∈ Z, we have the bounds 65) for all l ∈ Z. In addition, for any pair of triplets for all l ∈ Z, the quantities satisfy the bounds Proof. This follows from Lemma's 3.2-3.3 upon inspecting the definitions of S and R.
We are now in a position to state our main center manifold result. The fact that this manifold is two dimensional is related to the observation that the linear problem has the constant solutions (Φ * , 0) and ([∂ z φ z ] z=0 , Φ * ), as we will see in §5.
such that the following properties are satisfied. for κ → 0 and c → c * , uniformly for ϑ ∈ R.
(ii) The functions f θ and f κ are C r−1 -smooth. In addition, we have the behaviour as κ → 0 and c → c * .
(iv) Pick a c ∈ (c * − δ, c * + δ) and consider a triplet of sequences that satisfies (3.42) and admits the bound for all l ∈ Z. Then upon writing together with the difference equation are both satisfied for all l ∈ Z.
Proof of Theorem 2.3. For explicitness, we assume that [∂ 2 ϕ d ϕ ] ϕ=0 > 0 and [∂ 2 z λ z ] z=0 > 0. Whenever c − c * > 0 is sufficiently small, the identity [∂ ϕ d ϕ ] ϕ=0 = 0 allows us to use a Taylor expansion to show that there exist ϕ − < 0 < ϕ + for which d ϕ− = d ϕ+ = c. Noting that we use (iii) of Proposition 3.5 to define two quantities κ ± = κ ϕ± for which we have the identity In addition, possibly after further restricting the size of c − c * , we can ensure that the inequalities hold for all κ ∈ (κ − , κ + ). As a consequence, for any such κ we have the bounds (3.85) In particular, for any κ ∈ (κ − , κ + ) one can apply the contraction mapping principle to the fixed point problemκ and obtain a unique solution inκ ∈ (κ − , κ). For any choice ofκ 0 ∈ (κ − , κ + ), the problem can therefore be iterated backwards and forwards with respect to l to yield a solution κ : Z → (κ − , κ + ). This solution is strictly increasing and satisfies the limits lim l→±∞ κ l = κ ± .

Preliminaries
In this section we obtain a number of preliminary results related to the constant-coefficient linear system In particular, we study the Fourier symbol ∆(z) associated to this system and obtain a representation formula for solutions that are allowed to grow at a small exponential rate. As a preparation, we introduce the notation for the linearization of the travelling wave MFDE (2.8) around Φ * and the corresponding projection onto the kernel element Φ * . In addition, for any z ∈ C we define the vector recalling the convention that sums where the lower index is strictly larger than the upper index are set to zero. By construction, this allows us to write s [e z· w] = e z· s z w (4.4) for any w ∈ H 1 . Upon introducing the linear operators ∆(z) : is invertible and satisfies the bound In order to gain insight regarding solutions to the homogeneous linear system we briefly discuss the maximal Jordan chain associated to ∆(z) at z = 0. In particular, we set out to construct an analytic function z → J (z) ∈ H 1 × H 1 with J (0) = 0 so that ∆(z)J (z) = O(z m ) (4.10) for the largest possible value of m.
As a preparation, we note that whenever z = 0. Recalling the definition (2.16), we hence see that (4.12) Upon introducing the notation A k = ∂ k z L z z=0 (4.13) for any integer k ≥ 1, we may differentiate (4.12) to find (4.14) In particular, we may write Using (HS3) we immediately see which allows us to pick J (0) = (Φ * , 0). This chain can be extended by exploiting the following preliminary identities.

(4.18)
Proof. Differentiating the definition L z φ z = λ z φ z , we obtain the identities (4.19) Evaluating these expressions at z = 0 we find (4.17). The deriatives (4.18) can then be obtained by recalling the normalization ψ * , φ z = 1 and using the fact that ψ * , L * y = 0 for all y ∈ H 1 .

Indeed, we now write
for a pair of analytic functions z → v(z), w(z) ∈ H 1 × H 1 . Using the identity we may exploit (4.17) to compute (4.22) In particular, we have achieved m = 2 in (4.10). This corresponds with the presence of two solutions to the linear homogeneous problem (4.9). However, it is not possible to achieve m = 3. Indeed, setting the O(z 2 ) term in (4.22) to zero, we obtain and hence Taking the inner product with ψ * , we may use (HM) to obtain the contradiction Our second main result confirms that there are no other linearly independent solutions to (4.9) that are bounded by e ηmax|l| . In addition, it provides a representation formula for solutions to the inhomogeneous system (4.1) that share such an exponential bound.
Whenever L z : H 1 → L 2 is invertible, a short calculation shows that the same is true for ∆(z) with It is hence crucial to understand the behaviour of L −1 z for small |z| > 0, which we set out to do by exploiting the Fredholm properties of L * .
As a preparation, we introduce the notation for the unique v ∈ H 1 that has ψ * , v = 0 and satisfies the problem This allows us to rephrase the identities (4.17) in a more explicit form.  Proof. These expressions follow from (4.17), noting that L qinv Φ * = 0.
At this point, it is natural to briefly turn our attention to the angular dependence of the waves (c ϕ , Φ ϕ ), which can be analyzed using techniques that are similar to those used above. Unfortunately, the expressions for the second derivatives are somewhat more involved. In order to accomodate this, we introduce the notation    For any p ∈ H 1 , we can compute Based on (4.14), we may hence write (4.47) Using these identities to evaluate the expressions (4.43)-(4.44) at ϕ = 0 readily leads to the identities (4.41). The deriatives (4.42) can then be obtained by using the fact that ψ * , L * y = 0 for all y ∈ H 1 . Item (iv) of Lemma 2.2 follows directly by comparing the first identities in (4.17) and (4.18) with those in (4.41) and (4.42).
We now construct a preliminary inverse for L z that behaves as z −2 as z → 0. As a preparation, we implicitly define the remainder expressions R L;i by writing whenever h ∈ L 2 and 0 < |z| < δ z .
Proof. We set out to seek a solution a solution to (4.50) of the form for some κ ∈ R and w ∈ H 1 that satisfies ψ * , w = 0. Writing we may use (4.17) to compute The identity (4.50) is hence equivalent to the system

(4.54)
Whenever |z| is sufficiently small, we may use the quantity to rewrite the first line of (4.54) in the form Substituting this into the second line of (4.54), we find The desired properties now follow from the fact that I +zM (z) is invertible whenever |z| is sufficiently small.
In the following result we explicitly identify the singular O(z −2 ) and O(z −1 ) terms in the expansion of L −1 z . In order to express these in a convenient fashion, we introduce the operator Γ * : L 2 → R that acts as so that we have Proof. Fix z = 0 together with h ∈ L 2 and consider the functions Performing the expansion L z v # = E #;0 + zE #;1 + z 2 R #;2 (z) (4.63) and demanding that R #;2 (z) is analytic in zero for # ∈ {A, B, C}, we can use Lemma 4.2 to compute together with and finally Summing these expressions we see that where we used (4.59) to simplify the second expression.
In particular, Lemma 4.6 allows us to write from which the desired expansion follows.
Proof of Proposition 4.1. For any δ z > 0 and sufficiently small η max > 0, (HS1) implies that L z is invertible for | z| ≤ η max and δ z ≤ | z| ≤ π. The result now follows from the expansion obtained for L −1 z in Lemma 4.7. We now proceed towards establishing the representation formula (4.32). To this end, we introduce the left-shift operator S that acts on sequences as (4.69) Lemma 4.8. Consider any sequence y ∈ BX µ,ν (H) and pick an integer k ≥ 1. Then we have the identities Proof. Upon computing j≥0 e −zj y j+k = j ≥k e zk e −zj y j = e zk j ≥0 e −zj y j − e zk k−1 j =0 e −zj y j (4.72) together with j≥0 e −zj y j−k = j ≥−k e −zk e −zj y j = e −zk j ≥0 e −zj y j + e −zk k j=1 e zj y −j , the identities (4.70) readily follow. In addition, we compute j≥1 e zj y −j+k = j ≥1−k e zk e zj y −j = e zk j ≥1 e zj y −j + e zk k−1 j=0 e −zj y j from which (4.71) follows.
Proof. For the first compoment, we compute

(4.79)
For the second component, we note that The desired expression follows directly from these computations, noting that the third and fourth component can be obtained by flipping σ A and σ B in the expressions for the second respectively first component.  Proof. The expressions follow from the direct computation we introduce the notation This allows us to take the two discrete Laplace transforms of the main linear system (4.1).
Proof. Writing V = (v, w) and H = (g, h), we may use Lemma 4.8 and Corollary 4.9 to compute which is equivalent to (4.87). The remaining identity (4.88) follows in a similar fashion.
For convenience, we introduce the linear operators together with the projections Proof. Pick any z ∈ Z for which L z is invertible. Upon introducing the representation Pick a sufficiently small η max > 0 together with a sufficiently large K > 0. Then for every 0 < η < η max there exists a C r−1 -smooth map that satisfies the following properties.
(i) Pick any 0 < η < η max and ϑ ∈ R. For any H ∈ BS η (H 1 × L 2 ), the function V = K cc η (ϑ)H satisfies (5.1) and admits the orthogonality conditions (ii) Pick any 0 < η < η max and assume that V ∈ BS η (H 1 × H 1 ) satisfies (5.1) with H = 0 for some ϑ ∈ R. Then there exists a pair (a 1 , a 2 ) ∈ R 2 for which we have (iii) For any 0 < η < η max and ϑ ∈ R we have the bound (iv) Pick any 0 < η < η max . Then for any pair (ϑ 1 , ϑ 2 ) ∈ R 2 , we have the identity (v) Consider a pair (η 1 , η 2 ) ∈ (0, η max ) 2 together with a function Then for any ϑ ∈ R we have K cc η1 (ϑ)H = K cc η2 (ϑ)H. (5.9) Our strategy is to exploit the representation formula derived in §4 for the unprojected problem In particular, we first use the Fourier symbols ∆(z) defined in (4.5) to construct an inverse in the sequence spaces where again H is a Hilbert space. This can subsequently be used to obtain an inverse in the spaces by exploiting the fact that interactions between lattice sites decay exponentially with respect to the separation distance.
(ii) We have the bound Λ inv (5.14) (iii) We have the explicit expression Proof. This follows directly from Proposition 4.1 and standard properties of the Fourier transform; see for example [29, §3]. Pick a sufficiently small η max > 0 together with a sufficiently large K > 0. Then for every pair 0 < η 1 < η 2 < η max and any In addition, for any Proof. Following the approach in [27,Lem. 5.8], we introduce the sequences In view of the convergence k∈Z H (k) = H ∈ 2 η2 (H 1 × L 2 ), (5.22) the boundedness of Λ inv η2 implies that also for all l ∈ Z. We now pick two constants η ± in such a way that By construction, we have Recalling the left-shift operator S defined in (4.69), we note that In particular, we have Here the last identity follows from Proposition 4.3, since the sequence together with the homogeneous problem We are now able to use item (ii) of Lemma 5.2 to compute for some C 1 > 0. Summing over k, we hence find 33) The result follows directly from this bound, possibly after decreasing the size of η max > 0.
For any H ∈ BS η (H 1 × L 2 ) we now introduce the splitting for some small > 0, which by construction implies that V = K up η;I H satisfies the unprojected problem (5.10) with ϑ = 0. In addition, Lemma 5.3 implies that V ∈ BS η (H 1 × H 1 ). In order to allow for any ϑ ∈ R, we introduce the operator In view of the orthogonality conditions (5.4), we finally write Pick a sufficiently small η max > 0 together with a sufficiently large K > 0. Then for any 0 < η < η max , any ϑ ∈ R and any H ∈ BS η (H 1 × L 2 ), the function V = K up η (ϑ)H satisfies the unprojected problem (5.10) and admits the orthogonality conditions In addition, properties (iii) -(v) from Proposition 5.1 are satisfied after replacing K cc η by K up η .
Proof. In view of the discussion above, the statements follow directly from the fact that the set of solutions to the homogeneous problem (4.9) in BS η (H 1 × H 1 ) is two-dimensional as a consequence of Proposition 4.3. A detailed discussion can be found in the proof of [27,Prop. 5.1].
We now set out to lift the results above from the unprojected system (5.10) to the full system (5.1). A key role is reserved for the summation operator J that acts on a sequence W as (5.40) with the usual remark that sums are set to zero when the lower bound is strictly larger than the upper bound. for the set of solutions to the homogenous version of (5.1). By relating this set to its counterpart for (4.9) we show that N cc η is also two-dimensional for small η > 0.
Lemma 5.6. Assume that (Hf ), (HΦ), (HS1)-(HS3) and (HM) are all satisfied. Then for all sufficiently small η > 0, we have the identification Proof. A direct computation shows that In particular, we find (I − P 46) which verifies that the right-hand-side of (5.43) is indeed contained in N cc η .
Conversely, let us consider a sequence W ∈ N cc η , which implies that

Slowly varying coefficients
In this section we study the properties of the bounded linear operator that for any sequence θ : Z → R acts as We are specially interested in cases where the sequence θ varies slowly with respect to l ∈ Z, which means (S − I)θ ∞ < δ θ (6.3) for some small δ θ > 0. Our first main result states that the kernel of Λ(θ) is again two-dimensional, provided that (6.3) holds. For technical reasons, we also extend the two basis functions for the kernel to situations where (6.3) fails to hold. Proposition 6.1. Assume that (Hf ), (HΦ), (HS1)-(HS3) and (HM) are all satisfied. Pick a sufficiently small constant δ θ > 0 together with a sufficiently small η > 0. Then for every θ : Z → R there exist two functions that satisfy the following properties.
(iii) Pick any 0 < η < η max and suppose that Λ(θ)V = 0 for some V ∈ BS η (H 1 × H 1 ) and θ : Z → R for which (S − I)θ ∞ < δ θ . Then we have the identity Our second main result constructs operators K η (θ) that can be seen as an inverse for Λ(θ) whenever (6.3) holds. Naturally, the kernel elements above obtained in the result above need to be projected out, which is performed in (6.9). Special care needs to be taken when considering the smoothness with respect to θ. Indeed, the smoothness criteria below are based on the arguments involving nested Banach spaces argument that are traditionally used to establish the smoothness of center manifolds; see for example [16,§IX.7]. We remark that the notation L (p) stands for bounded p-linear maps. Proposition 6.2. Assume that (Hf ), (HΦ), (HS1)-(HS3) and (HM) are all satisfied. Recall the integer r appearing in (Hf ) and pick two sufficiently small constants 0 < η min < 2rη min < η max . Then for every η min < η < η max , there exists a map that satisfies the following properties.
Our strategy is to exploit the inverses K cc η (ϑ) for the constant-coefficient problem (5.1) to introduce an approximate inverse for Λ(θ) by writing [K apx η (θ)H] j = pev j K cc η (θ j )H. (6.19) In order to turn this into an actual inverse, we need to establish bounds for the remainder term To this end, we introduce the coordinate projection π 2 [v, w] = w together with the sequence A short computation shows that In a similar spirit, we introduce the sequences and compute These computations allows us to obtain the identity (6.25) In order to formulate appropriate bounds for this expression, we introduce the notation cev l θ = ev l θ − 1θ l = θ l−σ * +1 − θ l , . . . , θ l+σ * − θ l . Pick a sufficently small constant η max > 0 together with a sufficiently large K > 0. Then for any 0 < η < η max and any H ∈ BS η (H 1 × L 2 ), the following estimates hold.
(i) For any sequence θ : Z → R we have the bound Proof. As a consequence of the bound (5.6) and the smoothness of the map ϑ → K cc η (ϑ), there exists C 1 > 0 for which In particular, we see that for all |j| ≤ σ * we have The bound (6.27) hence follows immediately from the representations (6.22), (6.24) and (6.25).
In addition, (6.29) also implies that The second bound (6.28) readily follows from this.
(i) For any sequence θ : Z → R we have the bound (ii) For any η 1 > 0 and any pair of sequences θ A , θ B ∈ BS η1 (R), we have the bound After pickingδ θ > 0 to be sufficiently small, the bound (6.35) allows us to define the full inverse We remark that the normalization conditions (6.9) hold as a direct consequence of (5.4) and the construction of K apx η .
Turning to identify the kernel of Λ(θ), we introduce the operator for any sequence θ : Z → R. We now write Proof of Proposition 6.1. Item (i) follows from the fact that together with a similar identity for V B hom (θ). Item (ii) follows directly from the normalization (6.9). Finally, (iii) can be established by following the proof of [27, Lem. 6.4].

The center manifold
Our goal here is to construct and analyze a global center manifold for the system (3.42) that captures all the solutions where the pair (v, w) remains small. In particular, we set out to establish Proposition 3.5. While the main spirit of the ideas in [27, §7] can be used to establish the existence of the manifold, we need to take special care to identify the reduced equation that is satisfied on the center space.
The key issue is that we wish to recover a first order difference equation from a differential-difference system of order 2σ * .
Proposition 7.1. Assume that (Hf ), (HΦ), (HS1)-(HS3) and (HM) are all satisfied and pick two sufficiently small constants 0 < η min < η max . In addition, pick a sufficiently large constant K > 0 together with a sufficiently small constant δ v > 0 and write Then there exist C r−1 -smooth maps so that the following properties hold true.