Snaking without subcriticality: grain boundaries as non-topological defects

Non-topological defects such as grain boundaries abound in pattern forming systems, arising from local variations of pattern properties such as amplitude, wavelength, orientation, etc. We introduce the idea of treating such non-topological defects as spatially localised structures that are embedded in a background pattern, instead of treating them in an amplitude-phase decomposition. Using the two-dimensional quadratic-cubic Swift--Hohenberg equation as an example we obtain fully nonlinear equilibria that contain grain boundaries which are closed curves containing multiple penta-hepta defects separating regions of hexagons with different orientations. These states arise from local orientation mismatch between two stable hexagon patterns, one of which forms the localised grain and the other its background, and do not require a subcritical bifurcation connecting them. Multiple robust isolas that span a wide range of parameters are obtained even in the absence of a unique Maxwell point, underlining the importance of retaining pinning when analysing patterns with defects, an effect omitted from the amplitude-phase description.


Introduction and Motivation
Defects are of fundamental importance in the study of patterns as well as in materials science (Mermin, 1979;Cross & Hohenberg, 1993;Chaikin & Lubensky, 1995) and may take the form of domain walls, grain boundaries, dislocations and disclinations. These structures can be stationary or glide or climb or otherwise move through the pattern, and are present in one, two and three dimensions. One traditional approach to the study of defects is using an amplitude-phase decomposition in which a periodic pattern is described by a homogeneous amplitude and defects are associated with zeros of the amplitude. Since the phase at these locations is undefined these defects are known as topological defects and these have attracted the greatest attention (Mermin, 1979). Such defects are characterised by a non-zero 2 of 16 SUBRAMANIAN, ARCHER, KNOBLOCH, RUCKLIDGE FIG. 1: When neighbours keep their distance: gannets nesting at Muriwai beach, New Zealand. The image on the left shows the overall hexagonal ordering of the birds while that on the right shows a pentahepta defect (PHD). The green (pink) markers identify gannets with five (seven) neighbours instead of the usual six. Photo credit: Barbora Knobloch. topological charge, proportional to the integral of the phase around the defect. However, many defects are non-topological in character, but are nonetheless of fundamental significance (Chaikin & Lubensky, 1995). For these defects the phase is defined everywhere and the topological charge vanishes. Nontopological defects are frequently believed to be less important since isolated defects of this type may heal or otherwise disappear without interaction with another defect. However, this is not always the case, and if the defect has internal structure it may persist as a result of frustration or locking to the background state.
In the present work we adopt a different point of view and think of defects as spatially localised structures embedded in a background pattern. In general, defects represent a local disruption of the background state, and so can consist of holes where its amplitude is locally reduced, or local changes in the background wavenumber or its orientation. In most cases the background state is a periodic or crystalline state, but defects in quasicrystals or even disordered media also arise (Korkodi et al., 2013). Spatial localisation has been the subject of numerous recent studies (see Knobloch (2015) for a review) leading to the following picture of such structures. Defects may be described as holes in a pattern and these bifurcate from a pattern branch near its folds and are initially unstable. Both of these facts are a consequence of the presence of a Floquet multiplier +1 at such folds. In variational systems the resulting hole branch undergoes homoclinic snaking (Woods & Champneys, 1999) in the vicinity of a Maxwell point, the parameter value at which the free energy associated with the pattern is the same as that of a competing state. Homoclinic snaking is responsible for stabilising some of the hole states and is a consequence of the locking of a front between the two competing states to the heterogeneity or spatial structure on one or both sides of the front. On the basis of this picture we expect stable defect states to be present in regions of bistability, and associated with the presence of a Maxwell point, i.e., in the so-called snaking or pinning region. In particular, we expect holes to come in the form of infinite families of defect structures of specific symmetry and ever-increasing size, all of which coexist stably within the pinning region. In general, these structures either snake or lie on a stack of isolas (Beck 3 of 16 et al., 2009). Much of the above picture has been developed on the basis of detailed studies of the bistable Swift-Hohenberg equation (Burke & Knobloch, 2006), a prototypical pattern-forming equation exhibiting bistability between a pattern state and a spatially homogeneous state.
We point out that the traditional approach to analyse defects using an amplitude-phase decomposition ignores the possibility of pinning thereby collapsing the family of defects to a stationary heteroclinic front at the Maxwell point. The present work shows that pinning must be retained in studies of defects and that its effects in general extend over a broad range of parameter values. To demonstrate the importance of pinning and the associated multiplicity of defect structures we adopt the two-dimensional quadratic-cubic Swift-Hohenberg equation (SH23) and focus on one of the most important defect structures in two dimensions, the penta-hepta defect (PHD) in a hexagonal pattern. This classical structure consists of a bound state of a peak with five neighbours and a peak with seven neighbours, instead of the usual six neighbours in a defect-free hexagonal pattern. Figure 1 shows such a defect in a largely hexagonal arrangement of gannet nests at Muriwai beach in New Zealand.
In particular, we compute here three distinct defect structures and the associated isolas, thereby providing, compelling evidence for the coexistence of multiple defect structures in two dimensions. Earlier work in one dimension examined domain boundaries between two distinct spatially periodic states in SH357, the cubic-quintic-septic Swift-Hohenberg equation (Knobloch et al., 2019), and confirmed the presence of defect snaking in one spatial dimension. Snaking of grain boundaries between stripe patterns with different orientations in 2D has been investigated in the Swift-Hohenberg equation with cubic nonlinearity (Lloyd & Scheel, 2017). In our work, we investigate the Swift-Hohenberg equation with quadratic and cubic nonlinearities in the supercritical regime where hexagonal patterns are stable. The PHDs we study are the result of a local change in the orientation of the hexagonal pattern. Thus, these defects are the result of bistability between a hexagonal pattern and a rotated version of the same pattern. We expect similar structures in mass-conserving systems described by the conserved Swift-Hohenberg equation, or equivalently the phase-field crystal (PFC) model, since stationary states of this model satisfy the same stationary equation as that studied below (Thiele et al., 2013).
The bistable two-dimensional SH23 equation describes competition between hexagonal patterns and a trivial, spatially homogeneous state. The equation is variational with a Maxwell point between the hexagonal state and the trivial state. For this reason this equation admits both holes in a hexagonal pattern and hexagon patches embedded in the trivial background (Lloyd et al., 2008). Here we are interested in defects that are not associated with the trivial state and that are found in the supercritical regime. These defects are associated with grain boundaries between two hexagonal patterns, obtained by rotating the hexagonal pattern in part of the domain through 30 degrees relative to the rest. This procedure generates closed rings of PHDs. These are found to lie on distinct isolas despite the absence of a discrete Maxwell point between the two competing hexagonal states (both of which have the same free energy). We mention that an isolated PHD in a hexagonal pattern can be described using an amplitudephase description (Tsimring, 1995). This description shows that such defects can be thought of as a bound state of two dislocations, of opposite signs, in two of the three Fourier modes involved in the basic hexagonal structure. Such defects can be stationary or can climb (Tsimring, 1995). However, the description of such motion is expected to be strongly influenced by the presence of pinning between of the PHD and the surrounding hexagonal pattern. It is this effect that is responsible for the multiplicity of PHD structures in this system.
As already mentioned, equilibria of the Swift-Hohenberg equation are also equilibria of the PFC model, which is a simple theory for the crystallisation of matter. In this context, the results presented below apply to the defects observed at grain boundaries in polycrystalline materials. In particular, since we are considering a two-dimensional system, our results apply to the grain boundaries in two-dimensional solids such as graphene (Hirvonen et al., 2016(Hirvonen et al., , 2017, in which the defects are of great importance near to melting (Zakharchenko et al., 2011) and indeed affect the very nature of melting in two dimensions (Halperin & Nelson, 1978). This paper is organised as follows. In the next section we introduce the model equation we study, and the demodulation technique we use to isolate the defect state. In section 3 we describe our procedure for generating rings of PHDs. These structures are shown in section 4 and found to lie on distinct isolas. The spatial extent of the influence of the resulting PHD structures on the background hexagonal pattern is quantified. The paper concludes with a brief summary and some questions concerning the implications of multiple equilibria with defects for both pattern formation and material science that will, we hope, inform future work.

Model and Methodology
We consider the pattern-forming system SH23, a partial differential equation for a scalar field u(x,t), For some background on this equation, see e.g. Burke & Knobloch (2006). The parameters in this model are the linear growth rate µ and the strength of the second order nonlinear interactions Q. The domain chosen in all the calculations in this work, unless specified otherwise, is a square with side length L x = L y = 60π (30 wavelengths), discretised with 256 grid points in each direction. Time-stepping from a given initial condition is accomplished pseudospectrally using second-order exponential time differencing (ETD2) (Cox & Matthews, 2002). Numerical continuation is done using a pseudo arclength continuation similar to that employed in Subramanian et al. (2018). Note that Eq. (2.1) describes non-conserved gradient dynamics ∂ u/∂t = −δ F [u]/δ u, with the free energy functional Figure 2(a) shows the bifurcation diagram of extended hexagonal states H ± (so-called up-hexagons and down-hexagons, respectively) as a function of the chosen bifurcation parameter µ when Q = 0.75. A transcritical bifurcation occurs at µ = 0, generating H + for µ < 0 and H − for µ > 0. The H + states undergo a saddle-node bifurcation at µ = −0.05. Shown are two superposed branches differing in the orientation of their wavevectors in their Fourier spectrum. The blue line has solutions that include the wavevector (0, 1), whose real space pattern is shown in the inset with a blue border and labeled u 01 . Solutions along the dashed red line have wavevectors including (1, 0) whose real space pattern is shown in the inset with a red border and labeled u 10 . We observe that the root mean square values of both these extended states are equal, as expected since the patterns are just rotations of each other.
We can differentiate between the two extended hexagonal patterns if we consider one of them to be the reference state, say u 01 , and measure the norm of the difference between a given state and u 01 . In order to identify the parts of a given field u that differ from a chosen reference state, we create a filter matrix P in spectral space from the reference state u 01 in the following way.
First, we identify all the Fourier components that contribute nontrivially to the spectrum of the reference state by identifying the reciprocal lattice vectors (RLVs) of the reference stateû 01 , including wavevectors that are obtained via higher order interactions (see e.g. the discussion in chapter 4 of Chaikin & Lubensky (1995)). Figure 2 as blue circles. The dashed-dotted circle has been added to help identify the order 1 RLVs of magnitude |k k k| = 1, the radius of the circle. In order to calculate the distance of a given state from this hexagonal pattern, we seek to set all contributions from the above Fourier components to be zero in the pattern of interest. In this way we define the projection P whose application to the Fourier transform of a given field u(x) sets all the Fourier components identified in Fig. 2(b) to zero. The resulting field IFFT(Pû), where IFFT denotes the inverse fast Fourier transform, is now spectrally filtered and has no components along the Fourier components of the reference field, in this case u 01 . If we consider the same extended hexagonal patterns as in Fig. 2(a) and replot them using the root mean square of the spectrally filtered state u f il , calculated as we obtain Fig. 3. Since the new measure rms(u f il ) represents the departure of a pattern u from another chosen reference pattern, in this case u 01 , the state u 01 corresponds to the blue horizontal line with zero amplitude, while the state u 10 , shown in red, remains unmodified at the scale of the figure (compare the red dashed lines in Figs. 2(a) and 3). This is because in the case of extended periodic patterns, the higher order contributions to the Fourier spectrum of the pattern u 10 at wavevectors that are part of the higher order spectrum of the pattern u 01 are very small. The largest change arises from the removal of the bulk mode at |k k k| = 0. We mention that if the reference state is chosen to be the flat state, the measure (2.3) reduces to the usual rms(u).
In this section we have introduced the idea of spectral filtering to help visualize states with defects. In the rest of this work we focus on the formation of spatially localised patterns that involve pentahepta defects separating the two hexagonal states u 01 and u 10 . However, the idea of spectral filtering to visualize defects is more general and can be used to describe more complex grain boundaries.

Penta-hepta defects and grain boundaries
In the SH23 equation introduced in the previous section, we look for equilibria that involve penta-hepta defects. In order to promote the formation of hexagonal patterns, we set Q = 0.75 and choose a positive value of the growth rate, µ = 0.25. At these parameters, we use a combination of time stepping and numerical continuation to compute multiple coexisting states that involve spatial localisation of a patch of u 10 within a u 01 background. Figure 4 shows several such examples that are dynamically stable as confirmed by time-evolving small perturbations using Eq. (2.1).
Since it is difficult to identify visually the location of the defects in this state, we also look at the corresponding spectrally filtered fields u f il in Fig. 5. Here we see that the regions in Fig. 4 that are purely hexagonal with orientation u 01 are transformed into flat (green) regions in the spectrally filtered plot. Additionally, peaks in the u 10 patch coinciding with those of u 01 are also removed, see e.g. the flat/green regions in the middle of the hexagonal ring of peaks in Fig. 5(a). Note also that as we move across the three panels, the size of the u 10 patch in the middle occupies a larger part of the domain. However, all of the cases shown here possess symmetry with respect to reflections in the xand y-axes.
In order to correlate the location of the grain boundaries and associated penta-hepta defects in the u fields shown in Fig. 4 with the corresponding peaks in Fig. 5, we compare the two fields in panels (a): we zoom in to the centre of the domain in Fig. 4(a) and show the results in Fig. 6(a). Figure 6(b) shows the corresponding zoomed-in view of Fig. 5(a). In order to combine the information in both these figures, we estimate contours from panel (a) where the value of u is 65% of the maximum value over the domain. Then we superpose these contours (as black solid lines) over the u f il field as shown in Fig. 6(c). We see that the black contours line up with the peaks in the patch of u 10 while valleys between them, indicated in blue, indicate the locations of missing u 01 peaks. It is also possible to discern distortions in the background u 01 state in the vicinity of the grain boundary.
We seek to identify peaks that have five (seven) neighbours in order to locate penta-hepta defects automatically. In order to do this, we start with the data from the contours determined at 65% of max(u) as seen in Fig. 6(c). We determine the grid G comprised of the centres of each closed black  Fig. 4(c). The locations of penta-hepta defects are indicated using green star markers for peaks with 5 neighbours and pink star markers for peaks with 7 neighbours. circle (representing a single peak) in the contour. For points in the grid G , we create a 2D Delaunay triangulation and use it to determine the number of neighbours which surround each point in G . From this, we can immediately identify the locations of peaks that have five (seven) neighbours. Figure 7 reproduces Fig. 6(c) with additional green star markers indicating peaks with five neighbours and pink star markers indicating peaks with seven neighbours. This same procedure is adopted to create figures comprising the contours of the field u superposed on the image of u f il and overlaid with the locations of the penta-hepta defects in the movies associated with this paper.

To snake or not to snake
Having obtained multiple states with penta-hepta defects as equilibria of the SH23 equation at a given parameter value (as shown in Fig. 4), we seek to understand if these states are connected via homoclinic snaking, i.e., if they lie on the same solution branch. In order to determine this, we perform numerical continuation (Doedel et al., 1991) of the solutions of the equation Here, as with time stepping, the field u is discretised in space and nonlinear terms are calculated pseudospectrally. Pseudo-arclength continuation is performed with variable stepsize along the solution branch using a biconjugate gradient-stabilised method (van der Vorst, 1992) to approximate the action of the Jacobian of the system. We start with the state in Fig. 4(a) and vary the bifurcation parameter µ. Figure 8 shows part of the solution branch starting from this state. With decreasing µ, the measure rms(u f il ) decreases along the solution branch and the branch turns around at a saddle-node bifurcation at µ = 0.032, followed by two 9 of 16 further turning points at µ = 0.189 and µ = 0.173 before returning to the starting value of µ = 0.25. In order to compare the states on the solution branch before and after the first saddle-node bifurcation, we pick two nearby values µ = 0.1539 and µ = 0.1534, labelled as (b) and (c) in Fig. 8(a) and display the corresponding fields u f il in Fig. 8(b) and (c), respectively. Between these two solutions, we observe that the localised patch of penta-hepta defects has lost four peaks and has become slightly smaller. Such a change, involving addition/removal of part of the pattern across a saddle-node bifurcation, is similar to changes observed during homoclinic snaking of localised hexagon states in a u = 0 background in the subcritical regime µ < 0 (Lloyd et al., 2008). In addition the larger amplitude state (b) distorts the background state over a larger distance.
We track the full solution curve connected to the state in Fig. 4(a) and we find that for the chosen parameters we uncover an isola, which is displayed in Fig. 9. Panel (a) shows the complete isola in terms of rms(u f il ) plotted as a function of the bifurcation parameter µ over the range 0.032 < µ < 0.41. Supplementary material Movie 1 shows how rms(u f il ) changes as we move over this isola. Note that the identification of peaks with five or seven neighbours is based on the 2D Delaunay triangulation, which in turn is based on the contour curves for the field u. This means that any asymmetry in the location of the associated penta-hepta defects we observe at multiple locations along the isola could be a consequence of two levels of discretisation: one at the level of determining the centre of a contour peak and the second at the level of the automated determination of the number of nearest neighbours from the resulting 2D Delaunay triangulation.
In the next panel, Fig. 9(b), we see another version of the same bifurcation diagram with a different norm on the y-axis. In order to understand the relative cost in terms of free energy that is needed to create the observed ring of penta-hepta defects, i.e., to calculate the grain boundary free energy (Hirvonen et al., 2016(Hirvonen et al., , 2017, we first determine the free energy, F [u 01 ], of the reference state u 01 from Eq. (2.2). We subtract this free energy from the free energy, F [u], corresponding to the field in Fig. 4(a). With this difference in free energy costs as our norm, the resulting isola is as shown in Fig. 9(b). As shown enlarged in the inset, this curve displays swallow-tail behaviour at locations where we observe two consecutive saddle-node bifurcations in Fig. 9(a). This result is similar to the occurrence of swallow-tail structures in the snaking of two-dimensional quasicrystals in Subramanian et al. (2018). The difference in free energy, i.e., F [u] − F [u 01 ], gives the free energy of the interface and shows that this increases with µ, since the perimeter of the grain hardly changes as we go around the isola. Of course, the grain boundary free energy depends on the relative orientation (Hirvonen et al., 2016) of the neighbouring grains, but this also does not change around the isola. Figure 10 shows the consolidated version of the previous figure with the solution branches associated with all the three states shown in Fig. 5 for µ = 0.25. Panel 10(a) shows the variation of rms(u f il ) for the three isolas that result, all of which exist over a large range of the bifurcation parameter µ. From this panel, we deduce a direct correlation between increased rms(u f il ) values and the size of the central patch. Figure 10(b) shows the variation of the free energy difference between the extended hexagonal pattern and a state with a closed ring of PHDs, for the three different isolas. All of these free energy differences are positive, and increasing with the perimeter of the central patch. We observe the presence of multiple swallow-tail like structures along the isolas and we include an enlargement of one such swallow-tail from the isola connected to solution in Fig. 4(b). In this figure we also show additional solutions from each isola focusing on solutions at smaller values of µ. The locations of each solution is marked with a red circle in panels Figs. 10(a-b) and identified with the letters A-D. Compared to Fig. 7, we see that the variation of u f il in these panels is reduced, as reflected in the smaller difference between the maximum and minimum value of u f il in the associated colour bars. For a comprehensive look at how these solutions vary along different isolas, see the supplemental movie files. In order to understand the length scale over which the strain field generated by the defect decays, we look in Fig. 11 at the cross-section of the filtered field u f il along two different directions and do so for the field shown in Fig. 10, panel A. Cuts of abs(u f il ), i.e., the absolute value of the spectrally filtered field u f il , along the x− and y−directions both show large values at the centre of the domain with long oscillating exponential tails. From the slope of the decay we estimate the decay length L d to be 50, i.e., several times the pattern length scale (2π). We mention that exponential localisation is expected of the subcritical regime µ < 0, even in sheared situations ; ) but not necessarily for the localised axisymmetric spots that are present in SH35 even in the supercritical regime µ > 0 .
We next turn to an example where defect states do snake, instead of lying on a stack of isolas, as in Fig. 10. Our example comes from a related but mass-conserving phase field crystal (PFC) model that forms patterns at two length scales. The model describes crystallisation of soft matter into complex patterns and is closely related to the Swift-Hohenberg equation (2.1) suitably modified to generate instability at two length scales instead of just one. The model has been shown to produce both extended and spatially localised dodecagonal quasicrystals in 2D and icosahedral quasicrystals in 3D. See Subramanian et al. (2016) and Subramanian et al. (2018) for details. The model describes the evolution of a scalar density-like field u(x) and contains a parameter σ 0 that controls length scale selectivity, with more negative values of σ 0 corresponding to sharper length scale selectivity. As we follow the branch of fully nonlinear dodecagonal quasicrystal solutions (Fig. 12, leftmost inset) for increasing σ 0 the system starts to form localised hexagonal inclusions within the quasicrystal representing defects. Once such defects are present, the solution branch undergoes snaking as shown in Fig. 12 where the norm ||u|| is plotted as a function of σ 0 .
The RLVs for a quasicrystalline pattern are drawn from two circles of wavevectors and the higher order RLVs introduce new length scales at both ever larger wavevectors and at wavevectors close to zero, including wavevectors close to the first order RLVs (Rucklidge & Rucklidge, 2003;Subramanian et al., 2020). Therefore, the construction of the correct filter matrix P in the case of a quasicrystalline reference pattern is rather more complex. We therefore choose a fixed reference state at σ 0 = −8.94 and use this state to visualise the changes that occur in the field u with increasing σ 0 . We then identify points along different parts of the snaking curve (identified by red crosses in Fig. 12) where we plot the field u (left panels) side-by-side with the spectrally filtered rms(u f il ) calculated with respect to the reference state.
The left inset panels in each pair show that as we move higher up the snaking branch more of the region is covered by hexagonal patches, while the right inset panels display prominent peaks at locations where there is a change from quaiscrystalline to hexagonal local ordering. We note that these right panels are not as flat as the rms(u f il ) states that we obtained in the SH23 system for PHDs. This is due to several reasons: the complex way in which RLVs add new length scales for quasicrystalline patterns at higher orders, the choice of the reference state, etc. As a result this analysis remains incomplete but we use it here to highlight the fact that fully nonlinear equilibria with defects can indeed snake, forming multiple coexisting defect structures over a parameter interval, and leave a more detailed study of this system to future work.

Conclusion and future work
In this paper we have demonstrated by explicit computation on a model problem, the SH23 equation in the plane, that defects in a background hexagonal pattern come in discrete families, as expected of systems in which pinning plays an important role. In particular, we have shown that these defects lie on closed solution curves called isolas and so do not snake, in contrast to the behaviour of localised states in the subcritical parameter regime in SH23 (Burke & Knobloch, 2006;Lloyd et al., 2008). In other cases, such the two-scale PFC model, the defects are associated with the presence of classical snaking behaviour.
These results suggest that defects in patterns, whether periodic or quasiperiodic, behave much like localised structures on a spatially uniform background and that they can be studied using the same techniques from pattern formation theory as already used very effectively in the study of other states. Careful demodulation, that is, removal of the background pattern, appears to be key to identifying defects in simulations or experiments, although the nature of the background state brings its own properties into the behaviour of these states: pinning to background microstructure. This effect must be retained in any theoretical description of these defect states, precluding the use of the commonly invoked amplitude-phase decomposition.
The defects studied here are associated with grain boundaries and we have introduced these into the system by rotating the hexagonal pattern in a small part of the domain through 30 • . This angle determines the number of PHDs associated with the grain boundary, and hence the structure of this boundary (Hirvonen et al., 2016). Since this structure is associated with the inclusion of one pattern within another, both of which have the same free energy when separately covering the plane, there is no bulk phase free energy difference to drive the growth or shrinkage of the grain boundary between the two regions of different orientation. Instead, all the energy difference between the state with two grains and a defect-free pattern is associated with the boundary itself, and in particular the PHDs introduced into the pattern. It is noteworthy that this energy is positive ( Fig. 9(b)) indicating that such defect structures will not form spontaneously but that they require finite amplitude perturbations to initiate them.
The fact that the grain boundaries are associated with a positive free energy implies that the system will seek to eliminate such boundaries as much as possible -this is, by grain rotation. Since many of the observed states discussed above are dynamically stable, this implies that for the parameter values where these isolas exist, there is a free energy barrier that the system must surmount to decrease the inner grain size. In the materials science context these results show that in this parameter regime grain shrinkage is a thermally activated process. Note that PFC models have previously been applied to study grain rotation (Radhakrishnan et al., 2012;Hüter et al., 2017) and also other, more complex grain boundary and defect dynamics that occur, for example, when a material is sheared (Chan et al., 2010). Our results identify the parameter regime where grains in a polycrystalline material are dynamically stable and so identifies the regime in which grain rotation is an activated process. In addition, the different states identified by our approach will be attractors in the dynamics of a polycrystalline material equilibrating over time (Trautt & Mishin, 2012).
Related to this, we emphasise a point that can clearly be seen from in Fig. 9: the localised grain defect structures persist over a large range of values of the bifurcation parameter µ. We ascribe this fact to the absence of a unique Maxwell point between the two competing patterns (recall that both have the same free energy), thereby reducing the tendency towards depinning. In related systems possessing a Maxwell point, such as SH23 in the subcritical regime, the pinning region is limited by the onset of depinning that leads to the time-dependent invasion of the higher energy state by the lower energy state. We interpret Fig. 9 as suggesting that the presence of adjacent PHDs, and in particular the observed ring structure, stabilises individual PHDs and suggest that this fact may explain our inability to generate isolated PHDs in this system.
Our results also reveal that the localised defect structures influence the background pattern over a region larger than the defect itself. For example Fig. 8(b) clearly shows that the perturbations in u f il persist over a large part of the domain. This indicates that the perturbations due to the presence of the defects decay over length scales that are several times larger than the typical distance between peaks in the pattern. Figure 11 allows us to see precisely the decay length L d of these correlations, showing that 15 of 16 they falls off exponentially with distance from the defect structure.
This paper represents a first step towards understanding the multiplicity of defect states in two spatial dimensions. A number of questions remain, including • What is the smallest stationary and/or dynamically stable defect structure that can exist?
• Can we develop other topological measures to determine local order in a pattern with defect (Motta et al., 2018).
• How do the defect structures depend on the angle between the two hexagonal patterns (Hirvonen et al., 2017) and how does this angle affect the snaking properties of these states? For example, is it possible to uncover angles for which the defect states snake, rather than lying on disjoint isolas?
• What happens if non-variational terms are incorporated in the model equation? Does this generate motion of the defect state or do the defects remain stationary and stable?
• More generally, what are implications of these structures for conserved systems and materials science in general?
We leave these questions for future work on this rich and fruitful class of problems.