Approximate localised dihedral patterns near a turing instability

Fully localised patterns involving cellular hexagons or squares have been found experimentally and numerically in various continuum models. However, there is currently no mathematical theory for the emergence of these localised cellular patterns from a quiescent state. A key issue is that standard techniques for one-dimensional patterns have proven insufficient for understanding localisation in higher dimensions. In this work, we present a comprehensive approach to this problem by using techniques developed in the study of radially-symmetric patterns. Our analysis covers localised planar patterns equipped with a wide range of dihedral symmetries, thereby avoiding a restriction to solutions on a predetermined lattice. The context in this paper is a theory for the emergence of such patterns near a Turing instability for a general class of planar reaction-diffusion equations. Posing the reaction-diffusion system in polar coordinates, we carry out a finite-mode Fourier decomposition in the angular variable to yield a large system of coupled radial ordinary differential equations. We then utilise various radial spatial dynamics methods, such as invariant manifolds, rescaling charts, and normal form analysis, leading to an algebraic matching condition for localised patterns to exist in the finite-mode reduction. This algebraic matching condition is nontrivial, which we solve via a combination of by-hand calculations and Gröbner bases from polynomial algebra to reveal the existence of a plethora of localised dihedral patterns. These results capture the essence of the emergent localised hexagonal patterns witnessed in experiments. Moreover, we combine computer-assisted analysis and a Newton–Kantorovich procedure to prove the existence of localised patches with 6 m-fold symmetry for arbitrarily large Fourier decompositions. This includes the localised hexagon patches that have been elusive to analytical treatment.


Introduction
In this work, we are interested in stationary localised patches of planar patterns embedded in a quiescent state bifurcating from a Turing instability. These patterns have the fascinating property that outside of some compact region in the plane they resemble a homogeneous, or background, state, while inside the compact region they can take on intricate and striking spatial arrangements. Particular examples of such spatial arrangements that continually arise in applications are those of cellular hexagons and squares, as illustrated in Figure 1(a). These localised structures are known to occur in the quadratic-cubic Swift-Hohenberg equation (SHE) [44,57,58] given by where u = u(r, θ) is a function of the polar coordinates (r, θ) in the plane, ∆ := ∂ rr + 1 r ∂ r + 1 r 2 ∂ θθ is the polar Laplacian operator, 0 < µ 1 acts as the bifurcation parameter, and γ ∈ R \ {0} is fixed, as well as in two-component reaction-diffusion (RD) systems near a Turing instability of the form  where D is the diffusion matrix, f is a nonlinear function and µ is the bifurcation parameter. For example, in the von Hardenberg RD model for dryland vegetation [69] one can find localised hexagon patches, as depicted in Figure 1(b), which represent the density of vegetation in a water-scarce environment. Such localised cellular patterns are similarly found for other models of vegetation in arid climates [3,37], nonlinear optics [48,49,68], phase-field crystals [52,62], water waves [15,16], neural field equations [55], granular dynamics [4,65], binary fluid convection [45], and peaks on the surface of a ferrofluid [42,58]. Despite the prevalence and importance of localised planar patterns, little is known about them from a mathematical perspective.
To understand why this problem remains elusive, we first explore the theory of localised patterns in one spatial dimension. In one dimension, the existence of localised structures can be explained using spatial dynamics [22] or symmetry arguments [24], and even more complicated behaviour can be understood from bifurcation theory and energy arguments [9,17,19,71]. Naturally, one might question whether the aforementioned techniques could be extended to higher dimensions. Such an approach has yielded success in the study of localised planar fronts, where we now have a two-dimensional pattern (such as stripes or hexagons) but the localisation remains restricted to a single direction [28,41]. For example, Doelman et al. [28] proved the existence of modulated hexagon fronts connected to an unstable flat, or patterned, state. Proving this result required a centre-manifold reduction as well as finding connecting orbits to the resultant amplitude equations, utilising techniques from spatial dynamics, bifurcation theory, and geometric singular perturbation. However, this approach still requires that the localisation is in a single direction, such that the interface between states is a straight line. In fact, in [28] the authors note the difficulty in rigorously explaining the existence of fully-localised patches of hexagons, acknowledging that 'the approach used in this paper certainly fails'. Alternatively, one could formally impose an ansatz on a fixed hexagon spatial lattice and derive 2D amplitude equations. However, trying to find connecting orbits in these equations is even more difficult than the planar hexagon equations. Furthermore, this formal reduction still results in trying to find fully-localised solutions to a two-dimensional equation, and so does not serve to simplify the problem.
Despite this pessimistic outlook, recent attempts at understanding fully-localised planar patterns have yielded significant progress by focusing on patterns that are radially-symmetric [34,43,46,47,60]. In particular, the works [43,46] rigorously establish the existence of radially-symmetric solutions in the SHE (1.1). We note that, since we are exclusively interested in time-independent solutions, the SHE (1.1) takes the form of an RD equation as in (1.2) by simply setting u = (u, (1 + ∆)u). In the radially-symmetric case, i.e. u = u(r), the SHE (1.1) reduces to a fourth order ordinary differential equation (ODE) in r, meaning that one may interpret localised radially-symmetric patterns as ODE solutions which decay to zero as r → ∞. Using radial centre manifold theory, three types of stationary radially localised patterns have been shown to exist in (1.1) [43,46]: spot A which has a maximum at the core, spot B which has a minimum at the core, and rings which have their maximum/minimum away from the core. These patterns are shown to bifurcate from µ = 0 and for 0 < µ 1 each pattern is constructed via asymptotic matching of solutions in a core manifold, containing all small-amplitude solutions that remain bounded as r → 0, and a far-field manifold, which contains all small-amplitude solutions with exponential decay to zero as r → ∞. The construction of the spot A solution is the simplest to understand as a quadratic order expansion of the core manifold is matched with solutions to the linear flow in the far-field. We refer the reader to Figure 1(c) for a visualisation of the nondegenerate quadratic tangency between the core and far-field manifolds at µ = 0.
In this manuscript we attempt to move beyond the radially-symmetric case of u = u(r) by searching for steady-state solutions to two-component RD systems, including the SHE, that have a nontrivial dependence on the phase variable θ. To this end, we focus on localised dihedral patterns which depend on θ through a rotational m-fold symmetry. Precisely, such a pattern is comprised of peaks arranged in a compact region of the plane such that it is symmetric with respect to reflection in the x-axis and rotations of angle 2π/m about its centre, r = 0. We term these solutions D m patches, in reference to the dihedral symmetry group generated by the above rotations and reflections, and note that a localised square pattern would be a D 4 patch whereas a localised hexagon pattern would be D 6 patch. In [44], Lloyd et al. numerically observed localised D 6 solutions to (1.1) bifurcating from the flat state, which then connects to the curve of a spot A solution, as shown in Figure 1(d). Such numerical schemes involve solving a Galerkin system coming from an N -term truncated Fourier expansion of the solution in the phase variable θ. In Figure 1(e) we observe that different choices of N possess distinct solution curves of localised D 6 patches bifurcating from the flat state at µ = 0, which is a key motivation of this work. We emphasise that our work herein goes beyond just localised square and hexagonal patterns and accounts for both even and odd m.
Here we leverage the radially-symmetric analysis in [43,46,60] by following a similar approach to find D m patches bifurcating from the flat state at a Turing bifurcation point. We study approximations of small amplitude localised dihedral patterns by expanding a D m patch solution as the truncated Fourier series u(r, θ) = u 0 (r) + 2 N n=1 u n (r) cos (mnθ) , (1.3) where N ∈ N is the truncation order. It should of course be pointed out that by definition a D m -lattice pattern is invariant under rotations of 2π/m about its centre, meaning u(r, θ + 2π/m) = u(r, θ) for all (r, θ), as well as under reflection in the x-axis, meaning u(r, −θ) = u(r, θ) for all (r, θ), and so any such solution of (1.2) is captured by the Fourier cosine-series (1.3) upon letting N → ∞. However, the truncation to order N in the above Fourier series is necessary for our analysis, thus leading to the stipulation that we study approximate localised dihedral patterns here. Putting (1.3) into the RD system (1.2) results in a nonlinearly coupled system of non-autonomous ODEs in the radial variable r in terms of the Fourier coefficients u n (r). Our goal in this work is to use the radial centre-manifold theory developed in [43,46,47] to demonstrate the existence of exponentially decaying solutions u n (r) to our coupled ODEs, which through (1.3) lead to approximate localised planar patterns in the planar RD equations of the form (1.2). In particular, our analysis will be restricted to parameter values µ in the neighbourhood of a Turing instability and we will show that our planar patterns bifurcate from the homogeneous state undergoing such an instability. We will exclusively look for localised D m patch solutions that are analogous to the radially-symmetric spot A solutions from [43] and do not attempt to prove the existence of all types of small amplitude localised radial solutions. The major difference between our work and the work on radially-symmetric patterns is that our asymptotic matching between the core and far-field manifolds requires solving (N + 1) ≥ 2 nonlinearly coupled algebraic equations, while the latter has only a single nonlinear equation (N = 0). Our main results in the following section make this connection precise, and we show that these matching equations can be solved explicitly for N ≤ 4, while in the case that m is a multiple of 6 we employ a computer-assisted proof to help demonstrate the existence of solutions for all finite N 1.
Our approach possesses a number of advantages in studying the emergence of localised planar patterns. Firstly, our choice of Fourier decomposition (1.3) allows for patterns with a wide choice of possible symmetries. A common approach is to carry out a weakly nonlinear analysis where one derives slowly varying amplitude equations over a pre-determined periodic lattice in Cartesian coordinates by expanding u as where c.c. denotes complex conjugate, and k i are dual lattice vectors (e.g. a hexagon lattice would have 3 vectors k 1 = (−1, 0), k 2 = 1 2 (1, 3)). One can then formally derive a set of 3 complex Ginzburg-Landau type equations for the amplitudes A i ; see for instance [35,Chapter 9] for the hexagonal case. However, this approach has a number of significant hurdles to overcome, including rigorous justification for the reduction and then the proof of a localised solution to the amplitude equations. Rather than considering the localised version of a predetermined domain covering solution, our approach provides a local bifurcation theory for a multitude of dihedral patterns that extend beyond the well-studied stripes, squares and hexagons. Moreover, through this approach we are able to reduce a planar PDE problem to an algebraic matching condition, which can be solved numerically with relative ease. Even at small truncation orders, which we demonstrate can be solved without numerical assistance, our approximate solutions provide excellent initial conditions for numerical continuation and exhibit a strong likeness to examples of localised patterns observed in experiments. Finally, this approach has value not just in its results, but also in its limitations. We derive an upper bound for the bifurcation parameter in terms of the truncation order N , providing useful intuition for the difficulties encountered as N → ∞. We observe infinitely many localised patterns bifurcating from the trivial state, which highlights an issue in understanding fully-localised planar patterns and helps to motivate further study in this area.
The key limitation of this approach is that we cannot hope to explain the emergence of D m localised patterns in full planar RD systems. However, our solutions closely resemble those that have been documented in the literature, thus leading to the belief that they closely resemble true solutions of RD systems. Furthermore, we are able to say when localised patches of cellular patterns exist and bifurcate from the trivial state in a numerical scheme based on the finite mode Fourier decomposition (1.3). Through our analysis we are able to uncover new approximate localised solutions which can be used as initial conditions for numerical path-following routines, allowing one to continue localised solutions further into µ > 0 where the finitemode decomposition becomes a good approximation for the fully localised numerical solutions, and we observe exponential decay in the maximum of the amplitudes of u n . Numerically, we find these patches then undergo a process similar to homoclinic snaking, as discussed in [6,44]. See also the special issue on Homoclinic snaking [21] for a proper introduction to the topic and review of the literature.
In this work, we demonstrate that a local change of variables brings distinct RD equations into a normal form in the neighbourhood of a Turing bifurcation. Notably, the SHE (1.1) is a specific instance of the general RD equations (1.2) which can be written in this canonical form without any change of variables, where the only difference is that now the nonlinearity in (1.1) should be considered as a truncated Taylor expansion about u = 0, representing the deviation from the homogeneous state that undergoes a Turing instability. This instance in which the SHE arises allows it to be considered as a truncated normal form for systems undergoing a Turing bifurcation, and is exactly why so many pattern formation investigations have employed the SHE in the first place. Turing instabilities have been documented in spatially-extended models from chemistry [20,53], biology [31,40,63], ecology [39,50], and fluid dynamics [25,51], to name a few, and so our results may be applicable to more complicated systems if a suitable normal form reduction can be found. There has been some progress in this direction for radially-symmetric localised solutions, such as in neural-field equations [29] and on the surface of a ferrofluid [34], which suggests the approach presented here may also be extended in these cases. Hence, our work in this manuscript provides evidence for the emergence of localised dihedral patterns resulting from Turing bifurcations on planar domains, while also giving explicit forms for finding such patterns numerically. This paper is organised as follows. In Section 2 we present our main results. We begin with the necessary hypotheses to assume a non-degenerate Turing bifurcation is taking place in system (1.2) at some parameter value and then proceed to state our results for the Galerkin truncated system arising from assuming the form (1.3). In Section 3 we present our numerical findings, beginning with numerical verification of our analysis in §3.1 and then in §3.2 we provide numerical continuations of the localised dihedral patterns far into µ > 0 where they develop larger amplitudes and greater localisation. In Section 4 we define and quantify the core and far-field manifolds of our coupled radial ODE resulting from introducing (1.3) into (1.2). Furthermore, we reduce the problem of asymptotically matching these manifolds when 0 < µ 1 to solving a system of nonlinear matching equations in N + 1 variables. Section 5 is entirely dedicated to solving these matching equations. We begin by explicitly solving them for small values of N and then we demonstrate the existence of solutions to these matching equations for N 1. The latter is achieved by first employing a computerassisted proof (whose details are left to the appendix) to solve a limiting nonlocal integral equation and then using this solution to demonstrate the existence of solutions for large, but finite, N . The results of Sections 4 and 5 together are the proofs of our main results for the RD systems near a Turing instability, Theorems 2.2 and 2.4 below. Finally, we conclude in Section 6 with a discussion of our findings and some future areas of work.

Main Results
We begin with steady planar two-component RD equations, which we express as where u = u(r, θ) ∈ R 2 for the standard planar polar coordinates r ∈ R + and θ ∈ (−π, π]. Throughout we will assume that D ∈ R 2×2 is an invertible matrix and f ∈ C k R 2 × R, R 2 , k ≥ 3, describes the reaction kinetics of the system. The parameter µ plays the role of the bifurcation parameter. We assume that f (0, µ) = 0 for all µ ∈ R such that u(r, θ) ≡ 0 is a homogeneous equilibrium for all µ ∈ R. Since D is invertible, we can apply D −1 to (2.1) to normalise the diffusive term ∆u.
With the following hypothesis we make the assumption that (2.1) undergoes a Turing instability from the homogeneous equilibrium u = 0 with non-zero wave number. In what follows 1 n will denote the n × n identity matrix and for simplicity we will assume that this bifurcation takes place at µ = 0.
We now proceed by expanding D −1 f (u, µ) as a Taylor expansion about (u, µ) = (0, 0) to write (2.1) as In the above we have where Q and C are symmetric bilinear and trilinear maps, respectively. We note that any remainder terms do not affect the subsequent analysis since we are in the region where |u| and |µ| are small, and so we have neglected them in h.o.t., representing the higher order terms. Notice that Hypothesis 1 gives that the linearization of the right-hand-side of (2.2) about (u, µ) = (0, 0) has a zero eigenvalue, giving way to the Turing bifurcation. Let us now make the following non-degeneracy assumptions regarding this bifurcation.
The non-degeneracy conditions in Hypothesis 2 are important for our proofs in Section 4 and so we will formally describe how they arise. We suppose u = AÛ 0 + BÛ 1 and project onto each eigenvector; then, (2.2) can be transformed into the following normal form where A := (A, B) T , and Assuming |B| |A| 1, (2.3) can be reduced to the leading-order PDE which is analogous to the SHE (1.1). Here we have defined κ := Û * 1 , −C(Û 0 ,Û 0 ,Û 0 ) 2 with c 0 and γ as in Hypothesis 2. In [43] spot A solutions were found to exist in the SHE if and only if c 0 > 0 and γ = 0. We note that (2.2) is invariant under the transformationT : (u, Q) → (−u, −Q), and so the assumption γ = 0 is equivalent to assuming γ > 0, up to an application ofT . Hence, Hypothesis 2 provides nondegeneracy assumptions for the existence of localised spots in the SHE (2.4), which we relate to (2.2) in the neighbourhood of a Turing instability.
Remark 2.1. The SHE (1.1) fits into our RD framework and satisfies Hypothesises 1 and 2 by setting We see that M 1 is identical toM 1 in the RD normal form (2.3) for patterns close to a Turing instability, and so we will focus exclusively on the SHE (1.1) for our numerical demonstrations.
Under the above hypotheses, our goal in this work is to obtain approximate dihedral solutions to (2.1) for µ ≈ 0 that are bounded as r → 0, have lim r→∞ u(r, θ) = 0, and satisfy u(r, θ + 2π/m) = u(r, θ) for all (r, θ) and some m ∈ N. Previous works on the SHE have proven the existence of such solutions which are independent of the azimuthal component [43,46] and having seen (2.4) one should be convinced that these results can be extended to the more general system (2.1) with little issue. We aim to initiate the study of solutions of (2.1) that exhibit a nontrivial dependence on the azimuthal component, resulting in the property u(r, θ + 2π/m) = u(r, θ). Such patterns are invariant with respect to the symmetry group D m , generated by rotations of 2π/m in the plane about the origin and reflections over the horizontal axis. To this end, we introduce a D m , N -truncated Fourier approximation of solutions to (2.2) by Upon projecting onto each Fourier mode cos(mnθ), (2.2) is reduced to the (finite-dimensional) Galerkin system for all n ∈ [0, N ] and Since the case when N = 0 has already been thoroughly studied, we therefore restrict our analysis to the case when N ≥ 1 for the remainder of this manuscript and consider only nontrivial solutions of (2.6). We present our first main result of the manuscript detailing the existence of solutions to (2.6) for N = 1, 2, 3, 4, which represent approximate small patch solutions of the steady planar RD system (2.1).
Theorem 2.2. Assume Hypotheses 1 and 2. Fix m, N ∈ N and assume the constants {a n } N n=0 are nondegenerate solutions of the nonlinear matching condition a n = 2 for each n = 0, 1, . . . , N . Then, there exist constants µ 0 , r 0 , r 1 > 0 such that the Galerkin system (2.6) has a radially localised solution of the form for each µ ∈ (0, µ 0 ), n ∈ [0, N ], where ψ n (r) := k c r − mnπ 2 − π 4 and J mn is the (mn) th order Bessel function of the first kind. In particular, such localised solutions exist for N = 1, 2, 3 with the a n given in Proposition 5.4. Furthermore, in the case that N = 4 and 6 | m, these results again hold with the a n given in Proposition 5.6. Remark 2.3. We note that, much like in the axisymmetric case [43], the bifurcation of these localised patterns is entirely determined by the quadratic and linear coefficients γ, c 0 . Hence, the cubic nonlinearity of (1.2), which is essential in determining the bifurcation of localised one-dimensional patterns, has no effect on the emergence of these solutions.
The proof of Theorem 2.2 is broken down over Sections 4 and 5. In particular, Section 4 decomposes the Galerkin system on r > 0 into three disjoint intervals, making up the three components of the solution presented in (2.9). In each of these regions the dynamics can be captured by geometric blow-up methods. Section 4 presents the analysis for general N ≥ 1 and concludes by showing that solutions in the three distinct regions of the Galerkin system can be patched together by finding nondegenerate solutions of (2.8).
The resulting values of a n that satisfy these matching equations are exactly the a n presented in Theorem 2.2. Section 5 is entirely dedicated to further understanding and solving these matching equations. We present a number of symmetry results that help to eliminate redundant solutions of (2.8) and then proceed to solve these equations for small N by hand. The result of this work is the proof of Theorem 2.2.
As one can see, the matching condition (2.8) depends on the value of cos( mπ 3 ) and so there are four distinct cases depending on whether m is divisible by 2 and 3. In the case of N = 1, Proposition 5.4 below shows that there is a unique localised D m solution to (2.6) so long as 3 | m or 2 | m, up to a half-period rotation, with a peak at its centre if 6 | m and a depression otherwise. We emphasise that, although the a n are equivalent for any two choices of m leading to the same value of cos( mπ 3 ), the full solution u n (r) defined in (2.9) has a distinct m-dependent radial profile for each n = 0, . . . , N . This is evident in Figure 2 where we present the N = 1 solutions for rhombic (D 2 ), triangular (D 3 ), square (D 4 ) and hexagonal (D 6 ) patterns; although the rhombic and square patterns possess the same a n since cos( 2π 3 ) = cos( 4π 3 ) = − 1 2 , they produce fundamentally different solutions as a result of the Bessel functions J mn (r) in (2.9).
For N > 1, we note that matching condition (2.8) possesses an extra symmetry when 6 | m, which is analogous to associating 'bright' solutions u * (r, θ) with their 'dark' counterparts v * (r, θ) := U 0 (r) − u * (r, θ), where U 0 (r) is the localised axisymmetric solution found in [43,Theorem 2]. This means that the full set of solutions to (2.8) for 6 | m are described by a smaller subset of solutions than in the case when 6 m because more solutions can be recovered from their respective symmetries. We emphasise that our work in Section 5 shows that the matching condition (2.8) does not have a unique solution for 2 ≤ N ≤ 4, even after we quotient by those related by symmetry.
Solutions given by Theorem 2.2 for N = 2 are presented in  seen as a natural extension of the N = 1 solution, whereas the latter is fundamentally different from its lower-dimension counterpart. This suggests that the second solution is not related to the 'standard' localised D m pattern. Such a solution may be a transitory state that connects to another solution, such as the radial spot, or an artefact of our Galerkin approximation which only appears for a given finite truncation N . We leave such questions for a follow-up investigation.
For N = 3, Theorem 2.2 gives five distinct solutions to (2.6) when m is even or equal to 3. These distinct solutions can be seen in Figure 4 for m = 2, 3, 4, 6, as well as m = 5. Notice that we now find the emergence of a single D 5 solution. We emphasise that, much like when 2 | m or 3 | m, we obtain more solutions for 2 m, 3 m as N increases. However, there does not appear to be any connection between solutions at different truncation orders, and so there does not seem to be a 'standard' localised D m pattern when 2 m, 3 m. We further note that super-lattice structures, caused by a superposition of lattice symmetries, begin to emerge as solutions of (2.8). This is most apparent in the D 6 patterns of Figure 4 where, as well as the standard cellular lattice solution (H 1 ), we obtain solutions comprised of a hexagon patch surrounded by other hexagon patches (H 2 ) or faint triangular patches (H 3 ), respectively.
As N increases, it becomes more difficult to solve (2.8) explicitly, and so our analysis of N = 4 is focused on the simpler case when 6 | m. Theorem 2.2 (via Proposition 5.6 below) gives five distinct solutions to (2.6) which are presented in Figure 5 for the SHE. The super-lattice structures seen in N = 3 are more apparent, including a very striking pattern of localised triangular patches surrounding a hexagon in (H 3 ).
Solving the matching equations (2.8) for general N is a difficult task, especially if one seeks to identify all such solutions for a given N . To provide a more general result, we will restrict ourselves to the case 6 | m, which includes the case of hexagonal symmetries coming from the group D 6 . In this case the matching equations are slightly simplified since cos( mπ(n−j) 3 ) = cos( mπ(n−2j) 3 ) = 1 for each integer n and j. For low truncation orders N = 1, . . . , 4 our results in Section 5 below prove that there is a unique solution for which a n > 0 for all n ∈ [0, N ]. In particular, for m = 6 this solution corresponds to the 'standard' cellular hexagon patch (H 1 ) in Figures 2-5, where peaks are arranged in a uniform hexagonal tiling. As N increases, we numerically observe that the strictly positive solution possesses a scaling of the form a n = O(N −1 ) for all n ∈ [0, N ], and, as can be observed in Figure 6(a), we find that the rescaled positive solutions appear to converge to a continuous function as N becomes very large. Interestingly, for N 1 the matching equations (after rescaling the a n ) resemble a Riemann sum which in the limit as N → ∞ formally yields the nonlocal for each t ∈ [0, 1]. Hence, one expects that if α * (t) satisfies (2.10) then there exists a solution of (2.8) such that a n ≈ α * (n/(N + 1))/(N + 1), for each n = 0, . . . , N . We make this precise with the following theorem, for which the details of moving between the continuum equation (2.10) and the matching equation (2.8) are left to Section 5.    2 ) for each µ ∈ (0, µ 0 ). Precisely, there exists r 0 , r 1 > 0 such that 24c 0 k c a n µ where J 6m0n is the (6m 0 n) th order Bessel function of the first kind. The a n are such that for all ε > 0, there exists N ε > 0 such that for all N ≥ N ε we have where α * (t) is a positive and continuous solution of (2.10) for all t ∈ [0, 1].
Unlike Theorem 2.2, we are only able to identify a single solution of the Galerkin system (2.6) when N 1. This is limited by the fact that a n 's are obtained from a continuous solution to the continuum matching problem (2.10). To obtain the continuum solution we employ a computer-assisted proof, as detailed in Section 5.3, whose details are primarily left to the appendix. It is possible that similar computer-assisted proofs could produce other solutions to (2.10), in which case one can follow the work in Section 5.3 with relative ease to arrive at further solutions of the Galerkin system with large N . We emphasise that although the positive solution corresponds to standard localised cellular hexagons (H 1 ) when m = 6, Theorem 2.4 holds for any m = 6m 0 with m 0 ∈ N. For any m 0 > 1, the pattern given in Theorem 2.4 has 6m 0 -fold N=3 N=4 (a 2 ) symmetry and so corresponds to a localised quasicrystalline structure, such as those seen in [32,62]. We refer the reader to Figure 6

Numerical Investigation of Localised Patterns
In this section we present our numerical results for the localised patterns from the previous section. Throughout we will exclusively focus on the SHE (1.1) and restrict our larger-N system to be N = 10. This allows us to investigate patches for N = 1, . . . , 4 embedded into a higher-dimensional system, while also maintaining computation efficiency. The choice of N can be made considerably higher, however for our chosen radial domain, we have found that the choice of N = 10 is sufficient.
Our numerical procedure is described in Appendix A, which takes the radial domain to be 0 ≤ r ≤ r * , discretised into T mesh points {r i } T i=1 allowing us to numerically solve (1.1) using finite difference methods. We take Neumann boundary conditions at the outer radial boundary r = r * , and standard radial boundary conditions at r = 0 such that u(r, θ) is smooth at the origin. Following this, we employ a secant continuation code similar to [5] in order to continue solutions beyond the limited parameter regions from the results of the previous section. We fix γ = 1.6 in (1.1) as this is the same choice of parameters as seen in [44, Figure  24] where localised hexagon patches are observed undergoing snaking behaviour.
In what follows we first verify our analytical results from the previous section by investigating numerical solutions to (1.1) as µ → 0 + for N = 1, . . . , 4 and N = 10. We then conclude this section by continuing localised dihedral patterns beyond the small µ > 0 parameter regimes of our main results to observe snaking bifurcation curves.

Verification of Analysis
We recall that the theoretical results presented in Section 2 are for the parameter region 0 < µ 1, and so here we aim to support our theoretical results by numerically investigating solutions of (1.1) in the limit as µ → 0. As µ decreases towards 0, any localisation effects become weaker, meaning that solutions begin to grow in width. As a result, solutions will inevitably be affected by the width of the radial boundary r = r * as µ → 0. In order to minimise these boundary effects, we choose r * = 2000, with T = 6000 mesh points, for the rest of this subsection. We investigated various values for the mesh step-size δr with no observable differences in the results; as such, we choose T = 3 r * rather than T = 10 r * in order to improve computation speeds. We begin by introducing the 'numerical amplitude' u * n of each v n , defined by We investigate the ratios of u * n /u * 0 as µ → 0 and compare with |a * are solutions of the nonlinear matching condition (2.8). We complete this investigation for (a) D 2 , (b) D 3 , and (c) D 6 solutions to (1.1) for N = 3, as well as a D 6 solution for N = 4. In each case we solve for the 'standard' pattern, indicated by (R 1 ), (T 1 ), and (H 1 ) in Figures 2-5, respectively, and we present our findings in Figure 7. Notably absent are similar results for the standard D 4 pattern (S 1 ). The reason for this is that the predicted values of a * n for this D 4 pattern from Proposition 5.4 below are the same as for the D 2 pattern, and so their inclusion would clutter Figure 7. We do however comment that we have similarly confirmed our theoretical results for the D 4 case, despite them not being presented here.
In Figure 7 we present results for N = 3 (left) and N = 4 (right), for (a) D 2 , (b) D 3 and (c) D 6 solutions to (1.1). In each case, we have a straight horizontal line valued at 1 on the y-axis, representing the ratio of u * 0 /u * 0 , and so is naturally fixed for each choice of N and m. The horizontal dotted and dashed lines represent the ratios given by solutions of the matching equation (2.8); the (a) dotted, (b) sparsely-dashed, and (c) densely-dashed lines represent the cases when (a) 2 | m and 3 m, (b) 2 m and 3 | m, and (c) 6 | m, respectively. We observe that the ratio u * n /u * 0 of a D 2 solution, indicated by (a n ), tends to the associated dotted line from our theoretical results as µ → 0 + , for each n ∈ [0, N ]. The ratio u * n /u * 0 of a D 3 solution, indicated by (b n ), begins further away from the sparsely dotted line than (a n ) for moderate µ values. However, these solutions also appear to tend to their respective theoretical predictions as µ → 0 + , for each n ∈ [0, N ]. Finally, the ratio u * n /u * 0 of a D 6 solution, indicated by (c n ) in each figure, is plotted with respect to µ. We note that all of the values of (c n ) are less than 1, in contrast to the values of (a n ) and (b n ), and tend to their associated dashed lines from our theoretical results as µ → 0 + . We briefly comment on the reason why some ratios u * j /u * 0 in Figure 7 begin further from their theoretical predictions than others. While this cannot be fully understood without a more thorough numerical investigation, it can be partially explained by our definition of u * j . In (3.1), we do not account for any localisation in our solutions, which will naturally create some error for moderate values of µ. Furthermore, since the maximum value of the Bessel function J mn (r) is further from the origin for larger values of mn, the localisation causes a greater error than for smaller values of mn. Hence, in Figure 7 we see a larger discrepancy between numerical and theoretical  results for moderate values of µ as either m or n is increased.
In Figure 8 we present results for the 'standard' D 6 solution to (1.1) when N = 10. In Figure 8 (a) we plot N u * n against n/N , where u * n is defined in (3.1), and compare with N a n and α * (n/N ), where a n are numerical solutions of the matching condition (2.8) for N = 10, m = 6, and α * is a numerical solution of (2.10). The amplitudes are in quite good agreement, which could also be improved by further increasing r * or N . In Figure 8 (b) we plot log( v n T ) against n/N and observe exponential decay as n → N . This suggests that the coefficients of each Fourier mode cos(mnθ) decay exponentially as n increases, thus providing numerical evidence that the Fourier series (1.3) might converge to a continuous solution as N → ∞.

Continuation of Solutions
One of the major benefits of this work is its application to the numerical study of localised planar patterns. To illustrate, we can begin with an initial guess of the form Here β is a scaling term that accounts for our choices of (µ, γ), the values of a n are given by our theoretical solutions in Section 5.2, and the exponential term gives us an approximation for localisation in the far-field. In this parameter regime, solutions are more strongly localised and hence are smaller in width. Hence, for this subsection we choose r * = 100 with T = 1000 mesh points for computational speed. Then, by first solving the nonlinear matching condition (2.8) and substituting the subsequent solution a n into (3.2), we are able to construct very effective initial guesses for numerical continuation of such patterns. Furthermore, for moderate values of µ it is often sufficient to solve a low-dimension algebraic system, i.e. for N 0 = 1, 2, 3, which can be embedded into an initial guess of the form (3.2) with a higher dimension N 1 , where a n = 0 for all N 0 < n ≤ N 1 . We utilise this approach in order to find small-amplitude localised D m patterns and continue them to larger amplitudes. Three examples of continued localised D m patterns are presented in the Figures 9, 10, and 11 for m = 2, 3, 4, respectively.
The continuation of localised D 6 solutions was extensively covered in [44], and so here we focus on the novel D 2 , D 3 , and D 4 solutions. In Figure 9, we present our numerical results for a localised D 2 solution to (1.1). We begin with an initial guess of the form (3.2) with a 0 = −1, a 1 = − √ 2 and a n = 0 for 1 < n ≤ 10 such that (a 0 , a 1 ) satisfies the D 2 matching equation (2.8) when N = 1. Then, for µ = 0.1, the initial guess converges to a localised solution consisting of two spikes in close proximity, corresponding to the (R 1 ) solution in Figure 2. As µ increases in Figure 9, we observe that the D 2 solution undergoes snaking behaviour, where variations in µ cause solutions to gain extra peaks and grow in width.
Similarly, in Figures 10 and 11 we present our numerical results for a localised D 3 and D 4 solution to (1.1), respectively. In each case we solve the matching equation (2.8) for N = 1 and substitute our solution in the initial guess (3.2); notably, the initial guess for the D 4 solution is identical to the D 2 solution, other than a change of the value of m in (3.2). Then, for µ = 0.1 the initial guess converges to a localised solution consisting to three (m = 3) or four (m = 4) spikes, which is predicted by (T 1 ) and (S 1 ), respectively, in Figure 2. As µ increases we again observe that the D 3 and D 4 solutions exhibit snaking behaviour. This behaviour appears to be quite robust for larger dihedral symmetries as we have similarly observed such snaking at least for D 6 , D 8 , . . . , D 14 solutions as well. However, these results are not presented here for brevity. In one spatial direction this snaking behaviour is known as homoclinic snaking, described by homoclinic cycles in phase space, and is now well-understood [6, 8-10, 14, 18, 19, 42, 55]. However, in two or more spatial directions such snaking bifurcation curves are not well understood, and there continues to be significant interest in this area [2,7,13,21,23,26,54,64]. We note that in Figures 9, 10, and 11 we see that emerging peaks also appear to be subject to hexagonal packing, suggesting that the domain-covering hexagonal lattice may have a pivotal role in the structure of any observed localised patterns.
In the parameter regions of Figures 9-11 solutions continue into the strongly localised regime, where µ is moderately-valued, such that localised solutions resemble the N = 1 patterns in Figure 2. In order to capture more complicated patterns, we numerically solve (1.1) when N = 4 for localised D 6 solutions in the weakly localised regime with γ = 0.4. For this parameter choice, the hysteretic region of the bistable nonlinearity in the SHE is very small and solutions never become strongly localised. We present our results for localised D 6 patterns in Figure 12 where one can observe each of the distinct solutions (H j ) predicted in Figure 5 for j = 1, . . . , 5, and track the associated solution curves in µ-parameter space. Hence, we are able to numerically observe our theoretical solutions for (1.1) when N = 4 in a weakly localised parameter regime.

Localised Solutions to the Galerkin System
The goal of this section is to divide the dynamics of the Galerkin system (2.6) into separate regions over the independent variable r > 0 and then provide the necessary conditions for matching the solutions in each region together. We remark that satisfying the resulting matching conditions is left to the following section. We note that there exists a rescaling of the form for j = 1, 2, such that M 1 has a repeated eigenvalue of λ = −1 and (2.2) remains unchanged. Hence, without loss of generality, we will take λ = −1 throughout since the case when λ = −k 2 c can be recovered by inverting (4.1) at the end. To formulate the problem properly, we express (2.6) as the following first-order system, for each n ∈ [0, N ], where we recall that O 2 and 1 2 are the 2 × 2 zero and identity matrices, respectively.
Recall that the nonlinear sums in F n (U; µ) can equivalently be written as Our goal is to obtain exponentially decaying solutions of (4.2), which give way to the u n (r) in the truncated Fourier expansion (2.5) of the approximate solution u(r, θ).
Since we are interested in solutions which decay as r → ∞, we begin by noting that for each n ∈ [0, N ] we have lim r→∞ A n (r) = A ∞ . So, in the limit as r → ∞, the linearised dynamics of (4.2) decouple for each n and can be understood through the N + 1 identical eigenvalue problems The eigenvalue problem (4.4) reduces to solving the following equation for λ, which, after applying a suitable similarity transformation M → T −1 MT, where for any matrix M ∈ R 2×2 , the determinant (4.5) can be written as Here, we have reintroduced the linearly independent vectorsÛ 0 ,Û 1 ∈ R 2 introduced in Hypothesis 1, defined by following the rescaling (4.1), and equipped with adjoint eigenvectorsÛ * 0 ,Û * 1 such that Û * i ,Û j = δ i,j . The solutions of (4.7) are given by where the coefficient x ∈ R depends on M 2 . Thus, for µ = 0 the linear system has spatial eigenvalues λ = ±i, with double algebraic multiplicity. By Hypothesis 2 we have that Û * 1 , M 2Û0 < 0 and so, as µ increases off of 0, the system undergoes a Hamilton-Hopf bifurcation such that the eigenvalues split off of the imaginary axis and into the complex plane. We then expect localised solutions to bifurcate from the homogeneous state

Core Manifold
Far-field Manifold Exponential Decay as Bounded as ∞ Figure 13: The local invariant manifolds W cu − (µ) and W s + (µ) are constructed over sub-domains 0 ≤ r ≤ r0 and r∞ ≤ r < ∞, respectively. We choose r0 > r∞ and look for intersections between W cu − (µ) and W s + (µ) at the point r = r0.
in the region 0 < µ 1 since in this region each of the N + 1 linearised equations has a two-dimensional eigenspace of solutions that exponentially decay to zero as r → ∞.
In order to construct localised solutions to (4.2), which correspond to localised solutions of (2.6), we utilise local invariant manifold theory for radial systems, as seen in [60]. The general idea is as follows. In §4.1 we construct the set of all small-amplitude solutions to (4.2) that remain bounded as r → 0, which we call the 'core manifold'. This manifold is denoted by W cu − (µ), where we have used the notation for a local centreunstable manifold, and is constructed on the bounded sub-domain r ∈ [0, r 0 ], for some fixed r 0 1. Then, in §4.2 we construct the set of all small-amplitude solutions to (4.2) that decay exponentially as r → ∞, which we call the 'far-field manifold'. This manifold is denoted by W s + (µ), where we have used the notation for a local stable manifold, and can be constructed on r ≥ r ∞ , for some fixed r ∞ > 0. The value of r 0 can be freely chosen as long as it is sufficiently large and so we choose r 0 ≥ r ∞ such that the core and far-field manifolds overlap and can be matched at the coincidental point r = r 0 ; see Figure 13. In § 4.5 we show that this matching can be done under the assumption that we have a nondegenerate solution to a nonlinear system of equations, which we refer to as the 'matching equations'. In Section 5 we provide the properties and some solutions to these matching equations; combining these results with the work in this section leads to our main results in Section 2, since any function that lies on the intersection of both W cu − (µ) and W s + (µ) is, by definition, a localised solution of (4.2).

The Core Manifold
Here we will characterise the core manifold, W cu − (µ), for 0 < µ 1, which contains all small-amplitude solutions to (2.6) that remain bounded as r → 0. This is a local invariant manifold, and so we determine it on some bounded interval r ∈ [0, r 0 ], for a large fixed r 0 > 0. We begin by noting that at the bifurcation point µ = 0 the linearised behaviour about U = 0 of (4.2), given by decouples into (N + 1) distinct systems of the form for each n ∈ [0, N ]. Then, the linear system ∂ r v n = A n (r)v n has solutions of the form and J ν (r), Y ν (r) are ν-th order Bessel functions of the first-and second-kind, respectively. To see how we arrive at the above solutions for v n , we recall that the system d dr v n = A n (r)v n is equivalent to where v n = (u n , ∂ r u n ) T . By decomposing u n by the generalised eigenvectors we arrive at the following radial ODEs The second ODE in (4.14) is just the (mn) th -order Bessel equation, which has solutions of the form v (n) where the factor of 2 is included for future simplicity. Solutions to the first ODE in (4.14) can be written as v (n) where v p n (r) is the particular solution of (4.14) for v and so, after rescaling by a factor of π 2 , we see that the general solution to d dr v n = A n (r)v n is as stated in (4.11). Hence, the full 4(N + 1)-dimensional linear system (4.9) has solutions of the form Here we have used the notation [v] n to denote the n th element of a vector v. Furthermore, from the asymptotic forms in Table 1, V 1,2 (r) remain bounded and V (n) 3,4 (r) blow up as r → 0, for each n ∈ [0, N ]. Hence, we expect the set of solutions to (4.2) that remain bounded as r → 0 to form a 2(N + 1) dimensional manifold in R 4(N +1) for each fixed r > 0. Let us denote by P cu − (r 0 ) the projection onto the subspace of In what follows we will use the Landau symbol O r0 (·) with the meaning of the standard Landau symbol O(·) except that the bounding constants may depend on the value of r 0 .
Furthermore, the right-hand side of (4.20) depends smoothly on (d, µ), and the nonlinear functions Q m n (d 1 ) are defined as Proof. This statement is proven in a similar way to [43, Lemma 1]. We first note that the linear adjoint To see this, we note that d dr W n = −A T n (r)W n can be written as where W n = −r∂ r ( 1 r w n ), w n T . We decompose w n onto the generalised eigenvectors 2 (r)Û * 0 , and we arrive at the following ODEs Hence, we see that 1 r w (n) 1 solve (4.14), and so solutions take the form stated in (4.22). We choose the particular ordering and scaling of our adjoint solutions W (n) j (r) such that the relation such that each W for all n ∈ [0, N ], i, j ∈ {1, 2, 3, 4}, and is independent of r. For a given d = We first check that any solution U ∈ C([0, r 0 ], R 4(N +1) ) of (4.25) gives a solution of (4.2) that is bounded on r ∈ [0, r 0 ]. In the limit as r → 0, we see from Table 1 that the term [W (n) j (s)] 2 , denoting the second block element of W (n) j (s) in (4.22), is bounded by s (mn+3) for j = 3, and s (mn+1) for j = 4. Hence, the integrals multiplying the unbounded solutions V (n) j (r) are bounded by r (mn+4) for j = 3 and r (mn+2) for j = 4. Since |V (n) j (r)| = O r −(mn+1) for both j = 3, 4 as r → 0, we see that the right-hand side of (4.25) is continuously differentiable on its domain whenever U ∈ C([0, r 0 ], R 4(N +1) ). Since (4.25) is a specific case of a variation-of-constants formula, it is straightforward to check that U(r) satisfies (4.2).
We now need to check that any bounded solution U(r) ∈ C([0, r 0 ], R 4(N +1) ) of (4.2) satisfies (4.25). Taking a variation-of-constants formula for small solutions of (4.2), we find (4.26) Then, we introduce for j = 3, 4 so that we can write our small-amplitude core solution as In order to arrive at (4.20), we apply a Taylor expansion to (4.27) about |d 1 | = |d 2 | = µ = 0 and find where ν i,j,n := (2π) Notably, ν |i|,|j|,n is invariant under permutations of its indices and, for the restriction i + j = n, one can always find some a, b ∈ N 0 such that ν |i|,|j|,n = ν a,b,a+b . Computing the explicit value for ν a,b,a+b , we find that giving c (n) where Q m n (d 1 ) is as defined in the statement of the lemma. This completes the proof.
We conclude this section by noting the following. For sufficiently large r 0 , we can use Table 1 to write (4.20) in terms of its elements u n (r 0 ), v n (r 0 ) ∈ R 2 for each n ∈ [0, N ], where for each n ∈ [0, N ]. Here we have written y n := r 0 − mnπ 2 − π 4 , while the O r0 (·) remainders capture the higher order terms when |d| and |µ| are taken to be small.

The Far-Field Manifold
We now turn to characterising the far-field manifold, W s + (µ). To this end, we will introduce the variable σ(r) ≥ 0 to take the place of the 1 r -terms in (4.2). The result is the extended autonomous system with the property that σ(r) = r −1 is an invariant manifold of (4.35), and by definition this invariant manifold recovers the non-autonomous system (4.2). To construct the far-field manifold W s + (µ) we find the set of all small-amplitude solutions (U, σ)(r) to (4.35) such that U(r) decays exponentially as r → ∞. Following this, we evaluate at σ(r 0 ) = r −1 0 , thus restricting to the invariant subspace σ(r) = r −1 , such that U(r) is an exponentially decaying solution of (4.2).
Before attempting to find exponentially decaying solutions to (4.35), it is convenient to transform the system into the normal form for a Hamilton-Hopf bifurcation. We define complex amplitudes A := [ A n ] N n=0 , B := [ B n ] N n=0 , which satisfy the following relations for each n ∈ [0, N ], such that (4.35) becomes  In the coordinates { A n , B n } N n=0 , defined by (4.36), the core manifold (4.34) can be expressed as We apply nonlinear normal form transformations to (4.37), as seen in [43,60,61], to remove the non-resonant terms from the right-hand side. such that (4.37) becomes   Since the linearisation of (4.37) about ( A, B) = 0 decouples for each n ∈ [0, N ], the derivation of these transformations follows in the same way as for the radial problem. We use matched asymptotics in order to calculate the following leading order expansions, where we have used the symmetry of (4.37) with respect to the reverser R : ( A, B, σ, r) → ( A, − B, −σ, −r) in order to write down the higher-order terms containing σ.
Following [61, Lemma 2.6], we note that there exist smooth homogeneous polynomials {Φ n , Ψ n } N n=0 of degree 2 such that the change of coordinateŝ An explicit calculation verifies that every type of quadratic monomial belongs to the range of (D − i), and so we conclude that {Φ n , Ψ n } N n=0 exist. Hence, applying both (4.43) and (4.45) transforms (4.37) into  Having transformed our equations into the radial normal form (4.40) for a Hamilton-Hopf bifurcation, we also introduce the unconstrained variable κ(r) to take the place of √ µ. Then, replacing µ with κ 2 , the normal form (4.40) can be extended to the following system,  The rest of this section is dedicated to finding a parametrisation for exponentially decaying solutions in the far-field region and matching this with the core manifold (4.34) at the point r = r 0 . In order to find solutions to (4.50) where (A, B) decay exponentially as r → ∞, we consider the far-field region to be made up of two distinct sub-regions. These regions are the 'rescaling' region, where the radius r is sufficiently large such that r = O(µ − 1 2 ), and the 'transition' region, where the radius spans the gap between the rescaling region and the core region. We refer the reader to Figure 14 for a visualisation of these different regions. In the rescaling region, we perform a coordinate transformation in order to find exponentially decaying solutions to (4.40) for sufficiently small values of µ; these solutions are tracked backwards through r to the boundary of the rescaling region where they provide an initial condition for solutions in the transition region. Following this, solutions are tracked backwards in r through the transition region, starting at the previously-found initial condition and staying sufficiently close to the linear algebraic flow of (4.50), until they arrive at the matching point r = r 0 . See again, Figure 14 for visual reference. Finally, we perform asymptotic matching to find intersections between the core and far-field parametrisations at the point r = r 0 , thus defining localised solutions to (2.6). The dynamics in the aforementioned regions are described in the following subsection, followed by a subsection devoted entirely to the asymptotic matching.

The Rescaling Chart
We now define rescaling coordinates in order to find exponentially decaying solutions for sufficiently large values of r; that is, we introduce These coordinates are of a similar form to the standard rescaling coordinates seen in [46], except with the higher κ-scaling of A and B observed for axisymmetric spot A solutions in [43]. Then, we can write (4.50) in the rescaling chart as, for all n ∈ [0, N ]. We present the following lemma.
Lemma 4.6. For each fixed choice of r 1 > 0, there is a constant κ 0 > 0 such that the set of exponentially decaying solutions to (4.53) for s ∈ [r 1 , ∞), evaluated at s = r 1 , is given by for all µ < κ 0 , where a ∈ R (N +1) , a = 0, and Y ∈ R is arbitrary.
Proof. To begin, note that κ R acts as a parameter, and the subspace {κ R = 0} is invariant under the flow of (4.53). In this invariant subspace the dynamics reduce to which has an exponentially decaying solution for s ∈ [r 1 , ∞) of the form where a ∈ R (N +1) , a = 0, and Y ∈ R. Evaluating (4.58) at s = r 1 , this becomes Any exponentially decaying solution for κ R (r) ≡ ε 1 remains sufficiently close to the invariant subspace {κ R (r) = 0} since introducing κ R is a regular perturbation of (4.57). Therefore, setting κ R (s) = µ for sufficiently small values of µ. Hence, with the previous lemma we parametrised the far-field manifold W s + (µ) at a transition point r = r 1 µ − 1 2 , after which exponential decay is guaranteed to occur. By tracking the trajectories of (4.50) with initial condition (4.61) backwards in r, we are able to extend this parametrisation to the point r = r 0 , where we can then perform asymptotic matching to find intersections with the core manifold W cu − (µ). This analysis of the transition chart is covered in the following subsection.

The Transition Chart
In the previous subsection we parametrised the set of exponentially decaying solutions to (4.50) for r > r 1 µ − 1 2 . Now, in order to match these exponentially decaying solutions with the core manifold described in §4.1, we must track the trajectories associated with (4.61) backwards through the 'transition' region r 0 ≤ r ≤ r 1 µ − 1 2 . We therefore look to solve the initial value problem d dr ). This leads to the following result.  , and a ∈ R (N +1) and Y ∈ R were introduced in Lemma 4.6.
Proof. We begin by solving (4.62) for κ(r) and σ(r) explicitly, giving We then introduce the following transition coordinates for all n ∈ [0, N ]. Then, integrating over r 0 ≤ r ≤ r 1 µ − 1 2 , (4.66) becomes the integral equation, dp. (4.69) For sufficiently small values of µ, we can apply the contraction mapping principle to show that (4.69) has a unique solution (A T , B T ) in an appropriate small ball centred at the origin in C([r 0 , r 1 µ − 1 2 ], C 2(N +1) ); see, for example, [59]. Furthermore, we can express the unique solution to (4.69), evaluated at r = r 0 , as With the previous lemma we parametrised the far-field manifold W s + (µ) at the matching point r = r 0 . This now allows us to perform the asymptotic matching in order to find intersections of the core manifold W cu − (µ) with W s + (µ) for sufficiently small µ. This matching is covered in the following subsection.

Matching Core and Far Field
We begin by applying (4.39) to the core parametrisation (4.38), such that the core manifold W cu − (µ) can be expressed as for all n ∈ [0, N ]. Then, setting the far-field parametrisation (4.63) and the core parametrisation (4.71) equal to each other, we find that (4.73) Let us briefly summarise the variables and parameters included in (4.73). We recall that the core parametrisation is determined by the linear coefficients d 1 , d 2 ∈ R (N +1) , defined in Lemma 4.1, and the fixed complex phase Y 0 ∈ R. Similarly, the far-field parametrisation is determined by the linear coefficient a ∈ R (N +1) and phase parameter Y ∈ R, both introduced in Lemma 4.6. The other parameters in (4.73) are the matching points r 1 , r 0 > 0, the bifurcation parameter µ ∈ R, and the quadratic coefficient ν = 1 2 π 6 Û * 1 , Q(Û 0 ,Û 0 ) ∈ R; r 1 and r 0 are introduced during the construction of the core and far-field manifolds, whereas µ and Q(u, u) are given by the equation (2.2). Finally, we also note that the nonlinear terms Let us introduce the following coordinate transformations such that the leading order terms of (4.73) scale with the same order in µ. Eliminating common leading factors leads to for all n ∈ [0, N ]. When r 0 1 and 0 < µ, r 1 1, the leading-order system is given by (4.77) after dividing off the common factors of µ. We introduce the real diagonal matrix C ∈ R (N +1)×(N +1) given by C = diag(1, −1, 1, −1, . . . , (−1) N ). We prove in Lemma 5.1 below that Q m N (C m a) = C m Q m N (a) for all a ∈ R N +1 , but for now it suffices to notice this can be verified by direct calculation. Hence, we emphasise that (4.77) is invariant under the transformation (4.78) Then, taking real and imaginary parts, solutions to the leading order expression (4.77) are equivalently expressed as the zeros of the function (4.80) It should be noted that taking the imaginary part of the first equation in (4.77) results in the vector equation sin( Y ) a = 0. However, for a = 0, this reduces to obtaining Y such that sin( Y ) = 0, which is captured by the final component of G, defined in (4.80). The following lemma characterises the roots of G by relating them to fixed points of Q m N . Solving this fixed point problem is the focus of Section 5 which follows, but for the remainder of this subsection we show that these fixed points lead to matched solutions of (4.2) that are defined for all r ≥ 0 and decay exponentially to 0 as r → ∞. Proof. From the definition (4.79), roots of G correspond to having G 1 = 0, G 2 = 0, G 3 = 0, and G 4 = 0, as they are given in (4.80). Then, we immediately find that solving G 4 = 0 gives that Y = απ, for any α ∈ Z. Substituting Y = απ into G 3 , we find that d 2 = 0 is the unique choice that solves G 3 = 0. Then, solving G 1 = 0 results in d 1 = (−1) α a, and so G 2 = 0 implies where we have defined a = (−1) α a * . Hence, after applying the transformation (4.78), we conclude that (4.81) is a zero of G as long as a * = (a 0 , . . . , a N ) T is a fixed point of Q m N , completing the proof.
In order to match solutions from the core manifold to the far-field manifold for small values of µ, r 1 and r −1 0 , we require the Jacobian of G, denoted throughout by DG, to be invertible at the solution (4.81). This follows from the fact that the higher order (µ, r 1 , r −1 0 ) in the matching equations (4.75) constitute a regular perturbation of G when these parameters are taken to be small. Hence, when the Jacobian DG is invertible at V * , we can evoke the implicit function theorem to solve (4.75) uniquely for all 0 < µ, r 1 , r −1 0 1. Inverting the coordinate transformation (4.74), we find that where a * = {a * n } N n=0 is a fixed point of Q m N . Following our chain of arguments here, nondegenerate roots of G are used to solve the matching problem for 0 < µ, r 1 , r −1 0 1, which in turn leads to localised D m solutions to the Galerkin system (2.6). This therefore allows us to determine the leading order profile of u n (r) for these localised solutions in each region of r. For the core region, we substitute (4.83) into (4.25) to see that for r ∈ [0, r 0 ]. For the transition region r 0 ≤ r ≤ r 1 µ − 1 2 , solutions of (4.62) take the form Finally, for the rescaling region r ≥ r 1 µ − 1 2 , leading-order solutions take the form of the connecting orbit (4.58) for A R (s), B R (s) such that when r ≥ r 1 µ − 1 2 . To recover the case when M 1 has a repeated eigenvalue of λ = −k 2 c , we invert the rescaling (4.1) and transform the arbitrary values r 0 , r 1 accordingly. Then, for each n ∈ [0, N ] the radial amplitude u n (r) has the following profile: uniformly as µ → 0 + . Hence, a localised D m solution to the Galerkin system (2.6) has a core profile of the form Hence, it remains to (i) identify fixed points of Q m N to give zeros of G, and (ii) verify that these zeros are nondegenerate to arrive at the main results in Section 2. The following lemma shows that nondegenerate roots of G lie in one-to-one correspondence with fixed points of Q m N , a * ∈ R N +1 , such that the matrix That is, a * is a nondegenerate fixed point of Q m N . As stated above, the focus of Section 5 is to determine such fixed points, while in this section we provide the sufficient conditions that imply the existence of localised D m patches to the Galerkin system (2.6). We therefore conclude this section with the following result. Proof. The Jacobian DG can be written as where 0 ∈ R (N +1) , O N ∈ R (N +1)×(N +1) denote the zero vector and square matrix, respectively, and subscripts on the differential D denote what variable the derivative is taken with respect to. Evaluating at the solution V * , defined in (4.81), we find that (4.93) (4.94) Hence, the proof is complete.

Satisfying the Matching Condition
In Section 4, and in particular §4.5, we showed that to prove the existence of localised D m patch in the N -truncated Galerkin system (2.6) we are required to identify isolated solutions of the matching problem a n = 2 N −n j=1 cos mπ(n − j) 3 a j a n+j + n j=0 cos mπ(n − 2j) 3 a j a n−j , ∀n ∈ [0, N ]. (5.1) To this end, our goal in this section is to identify nontrivial solutions of (5.1) which in turn provide the results from Section 2. We recall the notation of where we recall that 1 N is the identity matrix of size (N + 1) × (N + 1) and DQ m N denotes the Jacobian matrix of the nonlinear function Q m N . In the following subsections we detail the existence and general properties of nondegenerate fixed points of Q m N . We begin in §5.1 with some important properties of the mapping Q m N that allow us to quotient the search for fixed points by important symmetries coming from the system (2.5). Then, in §5.2 we provide explicit solutions up to these symmetries for small-layer patches, which completes the proof of Theorem 2.2. We then provide a theoretical analysis which details the existence of solutions to the matching problem with N 1 in §5.3, which completes the proof of Theorem 2.4. All proofs of solutions for the small-layer patches are left to §B in the Appendix.

Properties of the Matching Equations
We begin by noting some qualitative properties of the Fourier-polar decomposition (2.5), and how these affect the matching problem (5.1). These properties are summarised in the following lemmas and can be used to classify solutions of the matching equations up to symmetry. Transforming a fixed point of Q m N into another one by applying R is a consequence of the fact that one may rotate a localised pattern of (2.6) by π/m to obtain the same localised pattern with a different orientation.
Indeed, using the Fourier-polar decomposition we see that U m (r, θ + π/m) = u 0 (r) + 2 N n=1 u n (r) cos(mnθ + nπ) = u 0 (r) + 2 N n=1 (−1) n u n (r) cos(mnθ). for each ∈ [0, N ], such that any a * that satisfies Q m0 N0 (a * ) = a * also satisfies Q m Proof. Take b * := H i N a * to be as defined in the statement of the lemma. Then, if is not a multiple of i, On the other hand, if = ni for some n = 0, 1, . . . , N , we have where we have used the fact that a * is a fixed point of Q m0 Ni . This concludes the proof.
As with Lemma 5.1, the results of Lemma 5.2 represent a fundamental fact about solutions of (2.6). Indeed, a D mi pattern expanded in N Fourier modes can be seen as a D m pattern expanded in N i Fourier modes. To see this, note that if we may define which is analogous to the definition of H i N a * in Lemma 5.2. This property follows from the fact that D m is a subgroup of D mi for any i ≥ 1. Proof. We begin by noting that with m ∈ N such that 6 | m, we have cos mπ(n − j) 3 = cos mπ(n − 2j) 3 = 1 (5.12) for all integers n and j. Then, taking b * to be as defined in the statement of the lemma, for = 1, . . . , N we have where we have used the fact that a * is assumed to be a fixed point of Q m N . Turning to the case when = 0, we similarly have where we have again used the fact that Q m 0 (a * ) = a * 0 . Hence, we have shown that Q m N (b * ) = b * , proving the claim.
The result of Lemma 5.3 is less intuitive than those of Lemmas 5.1 and 5.2. From Lemma 5.3 we see that if a * is a fixed point of Q 6m N , then b * := a 1 − a * is also a fixed point of Q 6m N , where a 1 = (1, 0, . . . , 0). The fixed point a 1 corresponds to a radial spot U 0 (r) and so Lemma 5.3 implies that, for any localised pattern U 6m (r, θ) of (2.6), there exists a corresponding localised pattern of the form V 6m (r, θ) = U 0 (r) − U 6m (r, θ). Such solutions are analogous to dark solitons in nonlinear optics, where localised solutions appear as 'holes' on a nontrivial uniform background state [38]. If we think of U 6m (r, θ) as a perturbation from the trivial state, which we could call a 'bright' solution, then V 6m (r, θ) can be thought of as an inverse perturbation from the localised radial spot, which we correspondingly call a 'dark' solution.

Small-Layer Patches
Here we will provide the nontrivial solutions of the matching equations (5.1), up to the symmetries described in the previous subsection, for small-layered patches. In particular, we will focus on the cases when N ≤ 4. We note that for any N, m ≥ 1 we have that a = 0 and a = (1, 0, . . . ) are solutions of Q m N (a) = a, representing solutions U m (r, θ) = 0 and U m (r, θ) = u 0 (r), which lead to the well-studied radially-symmetric solutions of the Galerkin system and are therefore deemed to be trivial solutions of the matching equations. With N ≤ 4 the matching equations are of a low enough dimension that we can obtain nontrivial solutions of (5.1) explicitly. These investigations further give us information on the role of the lattice m in finding localised D m patches and the goal is to provide all solutions that are unique to the given value of N . That is, Lemma 5.2 shows that solutions for N = 1 can be extended to solutions with N = 2, 3, 4 and solutions for N = 2 can be extended to solutions with N = 4; hence, we only detail the solutions for the smallest value of N that they appear for. Similarly, Lemma 5.1 demonstrates that solutions map into solutions by applying R, thus giving an equivalence class of solutions. The propositions that follow in this subsection will only provide one representative from each equivalence class.
Finally, for the ease of notation, we will define C m := 2 cos mπ 3 , which gives That is, if m is a multiple of 6 then C m = 2, if it is multiple of 3 and odd then C m = −2, if it is even and not a multiple of 3 then C m = −1, and if it is not divisible by 2 or 3 then C m = 1.
We first present our results for N = 1, 2, 3 and m ∈ N. The proof of this result is left to §B. 1. If N = 1 we have Importantly, there are no nontrivial solutions when 2 m and 3 m.

If N = 2 we have
where the λ j are real roots of the m-dependent polynomial Importantly, there are no nontrivial solutions when 2 m and 3 m.
3. If N = 3 and 3 | m, we have Furthermore, for all m ∈ N, we have and λ j are real roots of the m-dependent polynomial x 5 + 116x 4 − 217x 3 − 113x 2 + 24x + 9, 2 | m, 3 m, , and so on. Hence, Proposition 5.4 lists all solutions up to these equivalences to significantly condense the notation.
We now turn our attention to the case when N = 4; here, we restrict ourselves to the case when 6 | m since the other cases are significantly more complicated. Then, the equation a = Q m 4 (a) is given by

(5.16)
We summarise our findings with the following proposition whose proof is left to §B and completes the proof of Theorem 2.2.
We note that the roots {λ j } 5 j=1 are such that < λ 4 < 1 < λ 5 , and so the ordered set of signs of the elements of a 4 j is distinct for each j = 1, . . . , 5.

Large N Layer Patches with 6 | m
Let us now focus exclusively on the case 6 | m, representing solutions of the N -truncated Galerkin system (2.6) with a D 6d -symmetry, for some d ≥ 1. With 6 | m we have that cos(nmπ/3) = cos(2nπ) = 1 for every integer n, and so the matching equations become a n = 2 N −n j=1 a j a n+j + n j=0 a j a n−j , (5.17) for each n = 0, 1, . . . , N . In this section we will investigate solutions to (5.17) with N 1. To this end, this section will be dedicated to determining a continuous solution α(t) to the integral equation 18) and using this solution to demonstrate the existence of solutions to (5.17) when N 1. Formally, one can see that upon defining a n = 1 N +1 α( n N +1 ) for each n = 0, 1, . . . , N , (5.17) takes the form which can be interpreted as a Riemann sum approximation of (5.18). In what follows we will make the correspondence between solutions of (5.18) and solutions of (5.17) with N 1 explicit.
To begin, we will consider the closure of the space of real-valued step functions on [0, 1], denoted X, with respect to the supremum norm α ∞ := sup for each α ∈ X. Since we are taking the closure, X is a Banach space, and the elements of the space X are referred to as regulated functions. Regulated functions can be characterised by the fact that each α ∈ X and all t ∈ [0, 1] we have that the left and right limits α(t − ) and α(t + ) exist, aside from α(0 − ) and α(1 + ) [27]. We note that the space of continuous functions on for all α ∈ X. Using the fact that elements of X are integrable, one finds that Q ∞ is well-defined, and Q ∞ maps the closed subspace C([0, 1]) back into itself. We present the following theorem.
Proof. The existence of α * ∈ X is achieved by a computer-assisted proof that follows the work of [66], with the details of the implementation left to Appendix C. The general idea is to identify the existence of a fixed point of Q ∞ (α) near a continuous piecewise linear approximation of the solution generated by solving (5.19) numerically. We use interval arithmetic to verify a number of bounds that give way to the fact that Q ∞ is contracting in a neighbourhood of our constructed approximate solution. With this, the Banach fixed point theorem thus gives that generating a sequence from any element in a neighbourhood of our approximate solution by continually applying Q ∞ converges to a locally unique fixed point, giving the existence of α * . Now, since Q ∞ maps C([0, 1]) into C([0, 1]), initiating the sequence with our continuous linear approximation of the solution gives that each successive iterate is again continuous. Since C([0, 1]) is a closed subspace of X, it follows that the locally unique fixed point is also continuous since initiating the sequence with our approximate solution would make it the limit in the · ∞ norm of a sequence of continuous functions. The non-degeneracy of this fixed point is a consequence of the local uniqueness of the solution and is proven by following the arguments in [66,Lemma 5.2]. The proof of positivity comes from the fact that our constructed numerical approximation is positive and the fixed point is close to it, as proven in [66,Lemma 4.3]. This concludes the proof.
Theorem 5.7 gives the existence of a fixed point α * of Q ∞ , which can equivalently be viewed as a root of the nonlinear function G(α) := α − Q ∞ (α). Furthermore, the Fréchet derivative of G about α * , denoted DG(α * ) : X → X, is the bounded linear operator acting on f ∈ X by 22) for which Theorem 5.7 gives that this operator is invertible on X. Similarly, we introduce the functions a j a n−j , n = 0, 1, . . . , N. (5.23) That is, F N (a) = a − (N + 1)Q m N ((N + 1) −1 a), so that roots of F N are exactly solutions of (5.19). In the remainder of this section we will prove the following proposition, which details that the fixed point α * of Q ∞ can be used to closely approximate roots of F N for sufficiently large N ≥ 1, which themselves lie in one-to-one correspondence with solutions of (5.17) after rescaling a → (N + 1) −1 a for each N ≥ 1. In the remainder of this section we will take all finite-dimensional vector norms to be the maximum norm, or ∞-norm, which returns the maximal element in absolute value of the vector in analogy with the norm on X. To emphasise the difference between finite-and infinite-dimensional vector norms, the norm on X will keep the ∞ subscript, while the norms on finite-dimensional spaces will not have a subscript.
Theorem 5.8. Let α * be as given in Theorem 5.7. Then, there exists N 0 ≥ 1 such that for all N ≥ N 0 there exists a positive vector a * To prove Theorem 5.8 we will make use of the following lemma, which is a variant of the Newton-Kantorovich theorem, coming from [12].
Lemma 5.9. If F : R d → R d is smooth and there are constants 0 < κ < 1 and ρ > 0, a vector w 0 ∈ R d , and an invertible matrix A ∈ R d×d so that then F has a unique root w * in B ρ (w 0 ), and this root satisfies We now present the following lemmas which will then allow for the application of Lemma 5.9 to arrive at the proof of Theorem 5.8.
Lemma 5.10. Let a N be as in Theorem 5.8. Then, for any ε > 0, there exists N 1 ≥ 1 such that F N (a N ) < ε for all N ≥ N 1 .
Proof. First note that we can equivalently write F N in (5.23) as F N (a) n = a n − 2 N + 1 a 0 a n − 2 N + 1 N −n j=0 a j a n+j − 1 N + 1 n j=0 a j a n−j , (5.24) for each n = 0, 1, . . . , N , by simply adding and subtracting the j = 0 term from the first sum. Then, the summations in F N can be seen as left Riemann sums for the integrals in G. Since α * (t) is continuous at all t ∈ [0, 1], it follows that α * is bounded and uniformly continuous, thus giving that the error between the integrals and the Riemann sums in F N are bounded uniformly for t ∈ [0, 1]. The rate of convergence depends only on the modulus of continuity of α * , which is fixed and independent of N ≥ 1. Since G(α * ) = 0, this means that the terms a n − 2 N + 1 N −n j=0 a j a n+j − 1 N + 1 n j=0 a j a n−j (5.25) can be made arbitrarily small, uniformly in n = 0, 1, . . . , N . Finally, the terms 2a 0 a n /(N + 1) that appear in (5.24) can be bounded by 2 α * ∞ /(N + 1), which converges to 0 as N → ∞. Hence, the triangle inequality gives us that by taking N sufficiently large we can make F N (a N ) as small as we wish. This therefore proves the lemma. Lemma 5.11. Let a N be as in Theorem 5.8. Then, there exist C > 0 and N 2 ≥ 1 for which the Jacobian Proof. Let us start by showing that for N sufficiently large DF N (a N ) is invertible. We will assume the contrary to arrive at a contradiction. Hence, we may assume that there exists an infinite subsequence of as a nontrivial element of the kernel of DF N k (a N k ) with the property that v k ∞ := max |v n | = 1. We will show that this implies that DG(α * ) is not invertible, which from Theorem 5.7 is a contradiction.
We may consider the elements a N k , v k ∈ R N k as step functions in X which take the value α * ( n N k +1 ) and v n+1 , respectively, on the interval [ n N k +1 , n+1 N k +1 ), for each n = 0, . . . , N k . By abuse of notation we will again denote these step functions as a N k , v k ∈ X for each k ≥ 1. By assumption we have that v k ∞ = 1 as an element in X as well. Furthermore, since the elements a N k and v k are constant on the N k intervals of length 1/(N k + 1) we have that the integrals in DG(a N k ) reduce to discrete sums so that for each t ∈ [ n N k +1 , n+1 N k +1 ), using (5.24), we have for all n = 0, 1, . . . , N k and k ≥ 1. We note that the remaining term in the above equation comes from the j = 0 terms being added and subtracted from the first summation in the definition of each component of F N . Therefore, for each k ≥ 1 we have where we have used the fact that a N k (t) as a function in X takes the value α * (n/(N k + 1)) for each t ∈ [ n N k +1 , n+1 N k +1 ). Since α * (t) is uniformly continuous, it follows that as k → ∞, giving that lim k→∞ DG(ā)v k ∞ = 0, (5.29) since N k → ∞ as k → ∞. We therefore have that the sequence {v k } ∞ k=1 ⊂ X is a Weyl sequence for the linear operator DG(α * ), giving that DG(α * ) is not invertible. This is a contradiction, thus proving that for N sufficiently large we have that DF N (a N ) is invertible.
To show that the operator norm of DF N (a N ) −1 is uniformly bounded for sufficiently large N , we proceed with a nearly identical argument to that which was performed above to prove invertibility. That is, we may assume the contrary, thus giving way to a subsequence {DF N (a N )} ∞ k=1 for which there exists vectors v k ∈ R N k +1 with v k ∞ = 1 such that DF N (a N )v k → 0 as k → ∞. Arguing as above, we can find that DG(a N k )v k ∞ → 0 as k → ∞, after extending the elements a N k and v k to step functions in X. Indeed, following as in (5.27) we get The rightmost term DG(a N k )v k ∞ can be expanded as in (5.26), where it can again be shown that it converges to zero as k → ∞ by following as in (5.27) and replacing the condition that DF N (a N )v k = 0 with DF N (a N )v k → 0. Therefore, we again find that {v k } ∞ k=1 ⊂ X is a Weyl sequence for DG(α * ), contradicting the result of Theorem 5.7 which details that DG(α * ) is invertible. This contradictions means that we must have that the operator norm of DF N (a N ) −1 is uniformly bounded for sufficiently large N . This concludes the proof.
We now use the preceding lemmas to prove Theorem 5.8.
Proof of Theorem 5.8. To prove this result we will apply Lemma 5.9. Let us fix κ = 1 2 . From the continuity of DF N (a) near a N , we may evoke Lemma 5.10 to take N sufficiently large to guarantee the existence of a fixed ρ > 0 sufficiently small to guarantee that for all a ∈ B ρ (a N ). This therefore satisfies condition (1) to apply Lemma 5.9. Then, Lemma 5.11 gives that there exists a C > 0 so that DF N (a N ) −1 ≤ C for all N ≥ N 2 , and so Since F N (a N ) → 0 as N → ∞, we may consider N sufficiently large to guarantee that thus satisfying condition (2) from Lemma 5.9. Hence, we have satisfied the conditions to apply Lemma 5.9, and so for all sufficiently large N we find there exists a * N ∈ R N +1 such that F N (a * N ) = 0 and a * N − a N ≤ 2 DF N (a N ) −1 F N (a N ) . (5.34) Since DF N (a N ) −1 F N (a N ) → 0 as N → ∞, it follows that a * N − a N → 0 as N → ∞. Finally, the non-degeneracy condition DF N (a * N ) is invertible is satisfied for sufficiently large N since these matrices can be made arbitrarily close to the uniformly invertible matrices DF N (a N ) by taking N large enough. This therefore concludes the proof.

Conclusion
In this paper, we have investigated a finite mode truncation of planar two-component RD systems in order to understand the existence of localised dihedral patches. In particular, we showed that through an application of radial normal form theory we could establish precise conditions that guarantee the existence of localised solutions to the reduced Galerkin system (2.6) in a neighbourhood of a Turing bifurcation. These conditions came in the form of a nonlinear algebraic condition, which we were able to solve by hand for small finite mode decompositions. This uncovered new (approximate) localised patterns that can be used to initialise numerical continuation schemes to trace their corresponding bifurcation curves in parameter space. We further demonstrated the existence of localised D 6m , m ∈ N, patches with N 1 that bifurcate off the trivial state at the Turing instability point.
Let us now briefly describe what happens in the limit N → ∞ in our analysis. We find from our 'Core -Transition -Rescaling' approach that the core boundary r 0 must be chosen such that (mN ) 2 r 0 < µ − 1 2 . Although we only provide the existence of a µ 0 , r 0 > 0 in Theorems 2.2 and 2.4, we see from the above that our analysis requires that µ 0 → 0 + and r 0 → ∞ as N → ∞. This means that as N becomes large, the region of validity for our analysis shrinks, thus giving no insight into the existence of these patterns in the case of a full Fourier decomposition, i.e. N = ∞.
It seems natural to use polar coordinates when studying fully localised patterns, especially when using ideas from spatial dynamics, but this creates several additional problems in our analysis. First, to the author's knowledge, there is no established theory for a global reduction of a radial PDE system to a Ginzburg-Landau-type amplitude equation. By using the 'Core -Far-field' matching from radial spatial dynamics, we sidestep this issue and only locally reduce to a complex amplitude in the far-field through our normal form analysis. To do this, we must remove the oscillatory phase of our solutions, which then introduces the above conditions for r 0 and µ. Another issue with using polar coordinates, which makes a direct proof even more difficult, is that the PDE system inherently becomes a small divisor problem due to the 1 r 2 ∂ θθ term in the Laplace operator, ∆. As a result, many of the analytical issues in this problem are reminiscent of recent work on small divisors for quasipatterns [11,36], where the region of validity also shrinks to zero as the level of approximation improves, and the wider literature in this area might provide a possible solution to our problem.
Despite these analytical issues, our numerical results lead to the conjecture that for any µ > 0 taken sufficiently small we expect that true localised dihedral solutions to (2.2), coming from the limit N → ∞, exist and bifurcate from the Turing instability point. In order to establish this existence and persistence of the localised states found in the finite-mode truncation in full planar planar RD systems one could potentially apply similar computer-assisted proofs to [66,67] that are initialised with our finite mode solutions. Extending these finite mode solutions to true solutions of RD systems is an important area of future work that we hope to report on in a follow-up investigation.
Although we have focused on parameter values in a neighbourhood of a Turing instability, numerical investigations have revealed that localised dihedral patterns can further be found bifurcating from curves of localised solutions away from the Turing point. For example, Figure 15 presents numerical continuations of localised D 2 and D 4 patches bifurcating from the spot A and localised hexagon curves. Such branch points manifest themselves as symmetry-breaking bifurcations from the main localised solution curves. It may be possible to characterise these bifurcations using techniques from symmetric bifurcation theory, but predicting where these symmetry-breaking bifurcations take place along the curves remains an open problem.
Beyond just extending our approximate solutions to true solutions of RD PDEs, this work opens up several interesting future research directions. First, note that we have only investigated the simplest possible type of localised radial solutions in (2.6). We expect there to also exist the equivalent of radial spot B and radial ring solutions that have been established for the SHE [43,46]. It would be interesting to see whether the existence of spot B or ring localised dihedral solutions reduces to satisfying a matching condition akin to (2.8  the investigation of our spot A solutions could potentially be extended to fluid dynamic systems such as the ferrofluid problem [34,42], integral equations [29,55], and more general, multi-component reaction diffusion systems [60] where a centre-manifold reduction will have to be first carried out. We emphasise that while we have shown that the matching condition (2.8) turns up in the wide class of RD systems satisfying our hypotheses, it may be possible to have other matching conditions for different systems. and follow the same procedure presented above. The Galerkin dynamical system would be solvable using an identical approach to the dihedral patterns studied in this paper, however it is unclear whether one could solve the resultant matching equations without fixing additional restrictions on the solution.
Finally, one could follow our ideas to examine the existence of localised structures in three spatial dimensions. In this case our Fourier decomposition would have to be changed to a spherical or cylindrical harmonic decomposition, which may provide insight into pattern formation in phase field crystal models [62]. In summary, with this work we have opened the door to the investigation of existence and persistence of localised solutions beyond just those in one spatial dimension or those that are axisymmetric.

Code availability
Codes used to produce the results in this paper are available at: https://github.com/Dan-Hill95/Localised-Dihedral-Patterns.

Data availability statement
All data used to produce the results in this paper will be made available upon reasonable request.

A Numerical Implementation
We now briefly describe our numerical procedure for finding localised D m solutions to the (N +1)-dimensional Galerkin system (2.6), given some fixed m, N ∈ N. Recall that all computations use the SHE (1.1) and so we precisely target our discussion on the Galerkin system derived from it. We discretise the radial variable r using finite difference methods and fix r * > 0 so that we define our radial domain to be 0 ≤ r ≤ r * . The result is a decomposition of the domain r ∈ [0, r * ] into T mesh points {r i } T i=1 . In order for the operator (1 + ∆ n ) 2 to be well-defined in (2.7), the Galerkin system is also equipped with the following boundary conditions        ∂ r u n = (mn)u n = 0, r = 0, ∂ r ∆ n u n = (mn)∆ n u n = 0, r = 0, ∂ r u n = ∂ r ∆ n u n = 0, r = r * , for each n ∈ [0, N ]. We introduce differentiation matrices D 1 , D 2 using the central finite difference formula for first and second order differentiation, respectively. By defining the matrix R as we can express ∆ n as the finite-difference matrix D n , where Here the first and last rows of the matrices D 1 , D 2 , and R have been chosen such that (A.1) is satisfied.
for all n ∈ [0, N ], where • denotes the Hadamard, or element-wise, product and V = (v 0 , . . . , v N ) ∈ R T (N +1) . For computational speed, we do not evaluate (A.2) in this form; rather, we introduce the block matrices We also introduce a block exchange matrix J k ∈ R k(N +1)×k(N +1) , defined as so that we can express the nonlinear terms in (A.2) as where ⊗ denotes the Kronecker product of two matrices. We use MATLAB's Newton trust-region solver fsolve, starting with an initial guess of the form (3.2). If an initial guess (V, µ * ) converges to a localised solution for (A.2) with a trust-region accuracy of 10 −7 , we then define (V 0 , µ 0 ) = (V, µ * ) and use this output as initial data for our numerical continuation scheme. To numerically approximate the branch of solutions in parameter space that this initial data belongs to we employ secant continuation. See [5] for details on the implementation of this continuation technique. The result is an approximation of a continuous branch of solutions ( u(r), µ)[s] which solves the Galerkin SHE for s ∈ R. This allows us to map out solution curves for localised solutions in µ-parameter space, resulting in the figures from Section 3.2.

B Proof of Propositions 5.4 and 5.6
Let us begin by proving Proposition 5.4. Throughout we will fix m ∈ N and we recall that C m = 2 if 6 | m, C m = 1 if 2 m and 3 m, C m = −1 if 2 | m and 3 m, and C m = −2 if 2 m and 3 | m. We also note that we can use trigonometric identities to find that C 2m = (−1) m C m , which will be used frequently throughout the proof to simplify expressions. First note that if a 1 = 0, we are left with only the first and third equations, which give We note that these solutions are guaranteed by the symmetries outlined in Lemma 5.2. Now, assuming that a 1 = 0, we can divide it off of the second equation in (B.2) to give Substituting this expression for a 0 into the third equation of (B.2) and rearranging gives Then, putting these expressions for a 0 and a 2 1 into the first equation of (B.2) and rearranging results in the quadratic equation involving only a 2 : The right hand side of the above expression is exactly the function p m defined in the statement of the proposition, and so the conclusion follows by solving for a 2 and back-substituting to find a 0 and a 1 . Moreover, in the case that 2 m and 3 m the quadratic equation above reduces to a 2 2 = 0, which gives only the radial spot solution to the matching problem. For both the N = 1 and N = 2 cases, one can verify that each solution a * is nondegenerate by explicitly computing 1 N − DQ m N (a * ). For each non-trivial solution we find that det[1 N − DQ m N (a * )] = 0, and hence we conclude that these solutions are isolated roots of the matching condition.
We now turn to the case N (B.7) We note that a 1 = 0 implies that either a 2 = 0 or a 3 = 0. Taking a 1 = a 3 = 0, we see that (B.7) reduces to the N = 1 problem for a 0 , a 2 whose solution has already been found. Taking a 1 = a 2 = 0 instead gives that the remaining equations reduce to the N = 1 problem with m a multiple of 3. This situation was previously solved, and so we consider now the case a 1 = 0 by introducing the variable λ ∈ R, defined by We note that λ = 0 would imply that a 3 = 0, which in turn gives a 2 = 0. But, this would then imply that a 1 = 0, and since we are assuming a 1 = 0, we must have that at least one of a 2 and a 3 is non-zero, giving λ = 0. Substituting λ and α into (B.7), we find that a 0 = C m a 2 1 + (−1) m C m a 2 2 + (−1) m 2C 2 m λ 2 a 2 1 + a 2 0 1 = 2a 2 + C 2 m λa 2 + C m a 0 a 2 = C 2 m λa 2 1 + (−1) m C m a 0 a 2 + a 2 1 λ = (−1) m 2λa 0 + a 2 . (B.9) We input (B.9) into Wolfram Alpha's GroebnerBasis function [70], which reduces our system to a set of ideal polynomials equipped with the same roots. Then, finding a solution to this reduced system yields a solution to (B.9), which in turn gives us a solution to (B.7) after we recover a 3 = C m λa 1 . In order to find explicit solutions, we consider each case of m separately.
In the above examples (B.11),(B.13),(B.14) and (B.16), one can verify that each solution a 0 is a simple root of the respective polynomial p m (x), as defined in Proposition 5.4. Therefore each solution a 0 = λ j is isolated in R and, since a 0 corresponds to an element of a 3 j that solves (B.7), we conclude that each root a 3 j is also isolated in R N +1 . when a 1 = 0 and a 3 = 0. These systems are identical to the systems (B.1) and (B.2), respectively, when 6 | m, and so we can preclude this case since solutions to the matching equations are consequences of Lemma 5.2.
The set up is as follows. We consider the function G : X → X, defined in Section 5.3, given by Our goal is to find a solutionw with G(w) = 0, which by definition givesw = Q ∞ (w). Let us define a mesh ∆ M = {0 = t 0 < t 1 < · · · < t M = 1} such that t k+1 − t k = δt, ∀k ∈ [0, M − 1] and consider the subspace of X given by the linear splines S M ⊂ X with base points in ∆ M . Clearly, S M R M +1 since each element of S M can be uniquely determined by its value at the base points, and so the S M is a closed subspace of X. such that f (t k ) ∈ E k (f ) and f (t) ∈ E [t k ,t k+1 for all t ∈ [t k , t k+1 ]. It will also be useful to define the vector of intervals E(f ) such that (E(f )) k := E k (f ), for k = 0, . . . , M .
In what follows we demonstrate how one computes the Y and Z bounds. All bounds follow a similar derivation to [66] except the Z ∞ bound which is slightly different in our case since our mapping lacks some of the stronger differentiability properties to that studied in [66]. Sinceŵ M is a linear interpolant h 1 = h 2 ≡ 0, giving way to the proof of the lemma.
Remark C. 3. We note that in the numerical calculation of Y ∞ , we use the fact that such that the Y ∞ bound becomes order δt.