\section{Introduction}
\label{sec:introduction}

Pattern formation most commonly occurs with a single wavelength, as in for
example zebra stripes, Rayleigh--B\'{e}nard convection and the Taylor--Couette
flow \citep{Hoyle2006}. In these examples, there is a featureless basic state
that loses stability to waves with a nonzero wavelength as a control parameter
is increased. Typically the marginal stability curve, which separates stable
from unstable waves depending on their wavelength and the control parameter,
has a quadratic minimum.

In recent years, it has been recognised that pattern formation with two length
scales can lead to a wide variety of complex and interesting patterns, such as
superlattice patterns, quasipatterns and spatio-temporal chaos \cite[see for
example][and
references therein]{Castelino2020}. Having two length scales can arise in
different ways: in the Faraday wave problem with multi-frequency forcing, for
example, patterns with the two length scales arise in response to different
components of the forcing
\citep{Edwards1994,Topaz2002,Rucklidge2009,Skeldon2015}. Another possibility is
that the quadratic minimum in the marginal stability curve can change to a
quadratic maximum with two nearby quadratic minima at the two length scales.
This transition can occur via a quartic minimum, and is found in the
magnetised Taylor--Couette experiement \citep{Stefani2009,Mamatsashvili2019},
Lapwood--Prats convection \citep{Rees2013}, binary phase field crystals
\citep{Holl2020}, surface waves in ferrofluids \citep{Raitt1997} and nonlinear
optics \citep{Kozyreff2009}. This paper is concerned first with developing and
analysing a model that contains this transition in as simple a form as
possible, and second with investigating localised patterns in the model.

Problems with a single length scale where the pattern-forming bifurcation is
subcritical can have parameter intervals where both the featureless and the
patterned solutions are stable. In this case, it is possible to find localised
solutions consisting of a region of a spatially periodic pattern embedded in a
spatially homogeneous background
\cite[see][for reviews]{Dawes2010,Knobloch2015}. With two length scales,
there is a wider variety of possibilities, including having patterns with one
wavelength embedded in a pattern with a different wavelength. This phenomenon
has been observed in Rayleigh--B\'{e}nard convection in a long, thin channel,
or slot \citep{Hegseth1992} and has been explored in the context of
generalised Ginzburg--Landau models
\citep{Riecke1990,Raitt1995,Raitt1997,Kozyreff2009}.

A useful mathematical tool for studying pattern formation is the construction
of model equations that display qualitatively similar behaviour as the physical
system under consideration, but whose analysis is more tractable. Perhaps the
most ubiquitous of such model equations is the Swift--Hohenberg (SH) equation
\citep{Swift1977}, originally introduced as a model of thermal
fluctuations near the onset of Rayleigh--B\'{e}nard convection. It has been
used extensively in the study of localised patterns, starting with the work of
\citet{Hilali1995}, \citet{Crawford1999b} and \citet{Woods1999}. The equation
(in one dimension) is
 \begin{equation}
 \label{eq:SHE}
 \frac{\partial u}{\partial t} = \mu u -
      \left( 1 + \frac{\partial^{2}}{\partial x^{2}} \right)^{2} u +
      n_{2} u^{2} + n_{3} u^{3},
 \end{equation}
where $u(x,t)\in\mathbb{R}$ represents the pattern, $\mu$~is the driving
parameter, and $n_2$ and $n_3$ are parameters controlling the nonlinear terms
(typically $n_3=-1$). We consider equation~\eqref{eq:SHE} subject to periodic
boundary conditions on a domain $x\in[0,L]$.

The featureless (or trivial) solution $u=0$ is stable for~$\mu<0$. Small
amplitude perturbations of the form $e^{\sigma t+ikx}$ grow as
$\sigma=\mu-\left(1-k^{2}\right)^{2}$, so for $\mu>0$ the maximum growth rate is
at critical wavenumber $k=1$ (independent of~$\mu$), and if $\mu>0$, a range of
wavenumbers will grow exponentially until nonlinear effects become important.
The marginal stability curve is found by setting $\sigma=0$, so
 \begin{equation*}

 \mu = \left( 1 - k^{2} \right)^{2},
 \end{equation*}
which has a quadratic minimum at $k=1$.
At $\mu=0$ the trivial solution undergoes a pitchfork bifurcation, creating a
branch of spatially periodic solutions, which is stable if it is
supercritical, and unstable if not.

In large domains ($L\gg1$), and with small amplitude solutions ($u={\mathcal
O}(\epsilon)$), standard weakly nonlinear theory can be applied. The pattern
is written, with scaled space, time and parameter, as
 \begin{equation}
 \label{eq:GLE_scaling}
 u(x,t) = \epsilon A(X,T) e^{ix} + \text{c.c.} + \text{h.o.t.},
 \quad
 X = \epsilon x,
 \quad
 T = \epsilon^2 t,
 \quad
 \mu = \epsilon^2 \mu_2,
 \end{equation}
where $\text{c.c.}$ refers to the complex conjugate and $\text{h.o.t.}$
refers to higher-order terms. In this limit, the solvability condition for~$A$ at
third-order in~$\epsilon$ results in the Ginzburg--Landau (GL) equation
\citep{Cross1993}:
 \begin{equation}
 \label{eq:GLE}
 A_T = \mu_2 A + 4 A_{XX} + \left(3 n_3 + \tfrac{38}{9}n_2^2\right) |A|^2A,
 \end{equation}
where subscripts $T$ and $X$ refer to partial derivatives.
This equation governs the long-wavelength slow evolution of
the amplitude of solutions of~\eqref{eq:SHE}.
We define the coefficient of the nonlinear term
 \begin{equation*}

 n_{SH} = 3 n_{3} + \tfrac{38}{9} n_{2}^{2},
 \end{equation*}
which determines the criticality of the bifurcation at $\mu = 0$: if $n_{SH}<0$
the bifurcation is supercritical, and if $n_{SH}>0$, the bifurcation is
subcritical. In the supercritical case, the GL equation gives nonlinear
stability of striped patterns to long-wavelength perturbations \citep[the
Eckhaus instability, see][]{Eckhaus1965}. In the subcritical case, the GL
equation allows localised $\sech$-profile solutions to
equation~\eqref{eq:SHE} that can be continued in~$\mu$, leading to the
well known homoclinic ``snaking'' structure of localised solutions of the
Swift--Hohenberg and many other pattern forming
problems~\citep{Woods1999,Burke2006,Burke2007,Beck2009,Chapman2009}.

The GL equation~\eqref{eq:GLE} is one standard tool useful in the analysis of
the SH equation~\eqref{eq:SHE}. There are three others that we mention briefly.
First, the SH equation as written above is variational in time, and admits
a Lyapunov functional:
 \begin{equation}
 \label{eq:SH_Lyaponov}
 \mathcal{F}[u]
      = \int_{0}^{L} \left( - \tfrac{1}{2} \mu u^{2}
                            + \tfrac{1}{2} \left( \left( 1 + \partial_{x}^{2} \right) u \right)^{2}
                            - \tfrac{1}{3}n_{2} u^{3}
                            - \tfrac{1}{4}n_{3} u^{4} \right) dx.
 \end{equation}
With $n_3<0$, it can easily be shown that $\mathcal{F}[u]$ is bounded below and
that it is a decreasing function of time:
 \begin{equation*}

 \frac{d \mathcal{F}}{d t} =
      - \int_{0}^{L} \left( \frac{\partial u}{\partial t} \right)^{2} dx
      \leq 0.
 \end{equation*}
Equilibrium states correspond to stationary points of~$\mathcal{F}$, and those
coinciding with local minima of $\mathcal{F}$ must necessarily be stable. A
front connecting two patterns with different values of~$\mathcal{F}$ will tend
to move towards (and so eliminate) the pattern with the larger value. Localised
solutions are found near the Maxwell point, where the pattern and the zero
state have the same value of~$\mathcal F$, and near this point the difference
in~$\mathcal{F}$ is small enough that the front becomes pinned
\citep{Pomeau1986} to the underlying pattern.

The second tool is the observation that the steady
Swift--Hohenberg equation~\eqref{eq:SHE} admits a first integral in space.
Multiplying the time-independent version of~\eqref{eq:SHE}
through by $-u_{x}$, and integrating with respect to~$x$, yields
 \begin{align*}

 \mathcal{H} &= - \int \left( \mu u - u - 2 u_{xx} - u_{xxxx} + n_{2} u^{2} + n_{3} u^{3} \right) u_{x} dx \nonumber \\
 &= - \tfrac{1}{2} \left( \mu - 1 \right) u^{2} + u_{x}^{2}
    - \tfrac{1}{2} u_{xx}^{2} + u_{x} u_{xxx}
    - \tfrac{1}{3}n_{2} u^{3} - \tfrac{1}{4}n_{3} u^{4},
 \end{align*}
 and so $\frac{d\mathcal{H}}{dx}=0$. The quantity $\mathcal{H}$ is sometimes
referred to as the Hamiltonian for the steady version of~\eqref{eq:SHE}, since
there is a change of coordinates under which the system has Hamiltonian
structure. If there is a steady front connecting two patterns, the condition
$\frac{d\mathcal{H}}{dx}=0$ means that the two patterns must have the same value
of~$\mathcal{H}$.

The third useful tool is to note, again for the time-independent version
of~\eqref{eq:SHE}, that there is a Hamiltonian--Hopf bifurcation in space as
$\mu$~crosses zero. At the bifurcation point, there is a pair of double spatial
eigenvalues~$\pm{i}$, and the normal form can be written as a pair of
first-order ODEs in~$x$ for two complex variables. A bifurcation analysis
performed by \citet{Iooss1993} and \citet{Iooss1998a}, and extended by
\citet{Woods1999}, provides a geometrical interpretation of the solutions
of the normal form. In particular, there are parameter values where there are
solutions of the SH equation that are homoclinic to the origin as
$x\to\pm\infty$ \citep{Burke2007a}: these homoclinic orbits represent
localised solutions.

The existence and bifurcation structure of localised solutions in the
Swift--Hohenberg equation is now well understood: see
\citet{Dawes2010} and \cite{Knobloch2015} for reviews. In this paper, we
consider a generalised version of the Swift--Hohenberg equation that allows a
quartic minimum of the marginal stability curve ($\S$\ref{sec:model_eqn}).
Unfolding this quartic minimum, and using tools such as generalised versions of
the Ginzburg--Landau equation ($\S$\ref{sec:wnlt}), the Lyapunov function and
the first integral introduced above ($\S$\ref{sec:Lyapunov_first_int}), allows
us to identify parameter regimes where we can find patterns of one wavenumber
localised in a background of patterns with a different wavenumber
($\S$\ref{sec:localised}). We discuss the significance of our results in
$\S$\ref{sec:conclusion}. We include normal form calculations in
Appendix~\ref{sec:appendix}, but we have not found first integrals of the
normal form, and so we have not been able to put it to immediate use.

\section{The model equation}
\label{sec:model_eqn}

In this section we build a model equation to explore the unfolding of the
quartic minimum in the marginal stability curve. We start with the
SH~equation~\eqref{eq:SHE}, modified to allow a marginal
stability that can change from having one to having two minima, and then
add a selection of nonlinear terms.

\subsection{Linear terms}
\label{subsec:linear_terms}

As a starting
point, we consider a linear part of the PDE based on the polynomial
 \begin{equation*}

 p_0(K) = \left( 1 - K \right)^{4},
 \end{equation*}
where $K=k^{2}$, so $1-K$ will become $1 + \partial_{xx}$ in the model equation.
The quartic minimum at $K=1$ can be unfolded by adding two small terms to the
equation, yielding
 \begin{equation}
 \label{eq:polynomial}
 p(K) = p_0(K) + f_1 K + f_2 K^2,
 \end{equation}
with $f_1=f_2=0$ at the quartic minimum.
In principle, small terms $f_0$, $f_3K^3$ and $f_4K^4$ could be added as well,
but $f_0$~can be absorbed into the bifurcation parameter~$\mu$, $f_3$~can be
eliminated by making a small shift in $K$ by~$-\frac{1}{4}f_3$, and $f_4$ can be
absorbed by an overall scaling.

Before writing down the PDE model, we consider the transition from having one
minimum to two. The condition that $p(K)$ has two minima with a maximum in
between is the same as the condition that the derivative $p'(K)$ has three
distinct real
roots. Now
 \begin{equation*}

 p'(K) = 4K^3 - 12K^2 + (2f_2+12)K + f_1-4,
 \end{equation*}
and the condition that a cubic polynomial has three distinct real roots is that
its discriminant should be positive. The boundary, where the discriminant is
zero, occurs where
 \begin{equation*}
 f_1^2  + 4f_1f_2 + 4f_2^2 + \tfrac{8}{27}f_2^3  = 0,
 \quad\text{or}\quad
 f_1 = -2f_2 \pm \left(-\tfrac{2}{3}f_2\right)^{\frac{3}{2}}.
 \end{equation*}
In addition, when the discriminant is positive, $p(K)$~has two minima at $K_1$ and $K_3$
with a maximum at $K_2$ in between, with $K_1<K_2<K_3$, $K_1+K_2+K_3=3$ and
$K_1K_2K_3=1-\frac{1}{4}f_1$ (found
from the relationship between the roots and coefficients of the cubic $p'(K)=0$).
Manipulating the conditions $p'(K_1)=0$, $p'(K_3)=0$ and $p(K_1)=p(K_3)$ leads
to the conclusion that the two minima are equal when
 \begin{equation*}

 K_1 + K_3 = 2,\quad
 K_2 = 1, \quad
 \text{and}\quad
 f_1 + 2f_2 = 0,
 \end{equation*}
that is, the intermediate maximum is at $K=1$ and the two minima are equally
spaced on each side.
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_01
  \begin{tikzpicture}
  \begin{axis}[width=14cm,height=8cm,xmin=-0.12,xmax=0.12,ymin=-0.36,ymax=0.07,xtick={-0.1,0,0.1}]
  \addplot[color=black, mark=none, line width=1pt]
          table [x index=4,y index=0] {paper-cusp.dat};
  \addplot[color=black, mark=none, line width=1pt]
          table [x index=5,y index=0] {paper-cusp.dat};
  \addplot[color=black, mark=none, line width=1=0.5pt, dashed]
          table [x index=6,y index=0] {paper-cusp.dat};
  \addplot[color=blue,mark=*] coordinates{(0.095,0.022)};
  \addplot[color=blue,mark=*] coordinates{(0.00, 0.00)};
  \addplot[color=blue,mark=*] coordinates{(0.0340,-0.2269)};
  \addplot[color=blue,mark=*] coordinates{(0.000,-0.35)};
  \addplot[color=blue,mark=*] coordinates{(-0.0680414,-0.25)};
 \end{axis}
\draw (10.1,5.5) node[above] {$(i)$};
\draw (5.1,5.5) node[above] {$(ii)$};
\draw (7.0,1.75) node[above] {$(iii)$};
\draw (5.1,0.95) node[above] {$(iv)$};
\draw (1.6,2.25) node[above] {$(v)$};
\draw (9.0,-0.1) node[below] {$f_{1}+2f_{2}$};
\draw (-0.2,3.2) node[left] {$f_{2}$};
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,xshift=9.4cm,yshift=4.75cm,width=3cm,height=3cm]
 \addplot[mark=none,samples=20,color=blue,line width=1pt] table[x index=0,y index=1] {paper-cusp_curve_i.dat};
 \addplot[mark=none,color=black,line width=1=0.5pt] coordinates{(-0.1,1.1) (-0.1,-0.1) (2.1,-0.1) (2.1,1.1) (-0.1,1.1)};
\end{axis}
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,xshift=4.4cm,yshift=4.75cm,width=3cm,height=3cm]
 \addplot[mark=none,samples=20,color=blue,line width=1pt] table[x index=0,y index=1] {paper-cusp_curve_ii.dat};
 \addplot[mark=none,color=black,line width=1=0.5pt] coordinates{(-0.1,1.1) (-0.1,-0.1) (2.1,-0.1) (2.1,1.1) (-0.1,1.1)};
\end{axis}
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,xshift=6.3cm,yshift=1cm,width=3cm,height=3cm]
 \addplot[mark=none,samples=20,color=blue,line width=1pt] table[x index=0,y index=1] {paper-cusp_curve_iii.dat};
 \addplot[mark=none,color=black,line width=1=0.5pt] coordinates{(-0.1,1.1) (-0.1,-0.1) (2.1,-0.1) (2.1,1.1) (-0.1,1.1)};
\end{axis}
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,xshift=4.4cm,yshift=0.2cm,width=3cm,height=3cm]
 \addplot[mark=none,samples=20,color=blue,line width=1pt] table[x index=0,y index=1] {paper-cusp_curve_iv.dat};
 \addplot[mark=none,color=black,line width=1=0.5pt] coordinates{(-0.1,1.1) (-0.1,-0.1) (2.1,-0.1) (2.1,1.1) (-0.1,1.1)};
\end{axis}
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,xshift=0.9cm,yshift=1.5cm,width=3cm,height=3cm]
 \addplot[mark=none,samples=20,color=blue,line width=1pt] table[x index=0,y index=1] {paper-cusp_curve_v.dat};
 \addplot[mark=none,color=black,line width=1=0.5pt] coordinates{(-0.1,1.1) (-0.1,-0.1) (2.1,-0.1) (2.1,1.1) (-0.1,1.1)};
\end{axis}
   \end{tikzpicture}
 \endpgfgraphicnamed}
   \end{center}
 \vspace{-3ex}
 \caption{\label{fig:cusp} Five examples of the marginal stability
 curves $\mu=p(K)$ in different regions of the
 $(f_{1}+2f_{2},f_{2})$ parameter space.
 The curves (blue in black boxes) shown are:
 $(i)$~single minimum,
 $(ii)$~quartic minimum,
 $(iii)$~double minima, with the left minimum being the lower,
 $(iv)$~double minima, with both minima at the same height, and
 $(v)$~a single minimum and an inflexion point. The specific parameters for each
example are shown as blue dots.
The solid lines indicate where
the discriminant is zero (the derivative has a double zero), and the dashed
line ($f_{1}+2f_{2}=0$, $f_2<0$) indicates where the two minima exist and have the same height.}
\end{figure}

The next step is to convert the polynomial to the linear operator of the model
equation:
 \begin{equation*}

 \frac{\partial u}{\partial t} = \mu u - \left( 1 + \frac{\partial^{2}}{\partial x^{2}} \right)^{4} u + f_{1} \frac{\partial^{2} u}{\partial x^{2}} - f_{2} \frac{\partial^{4} u}{\partial x^{4}}.
 \end{equation*}
A mode $u=e^{\sigma t+ikx}$ has growth rate
$\sigma=\mu-(1-k^2)^4-f_1k^2-f_2k^4=\mu-p(k^2)$, connecting the dispersion
relation to the polynomial~$p(K)$ discussed above. Marginal
stability, when $\sigma=0$, occurs when $\mu=p(k^2)$, and Figure~\ref{fig:cusp}
shows examples of the marginal stability curves in the
$(f_1+2f_2,f_2)$ parameter plane. The discriminant is positive, and there is a
double minimum, within the cusp-shaped region below the solid curves, and the
two minima are equal on the dashed line.
The cusp, the point at which there is a quartic minimum with $f_{1}=f_{2}=0$,
represents a codimension-three bifurcation as $\mu$~crosses zero.

However, a model based solely on this would have a wavenumber for maximum growth
rate for solutions that did not depend on~$\mu$.
Regions of secondary instability of patterns
are organised around the curve of maximum growth rate, and it is important that
this curve is modelled correctly.
In pattern-forming problems with a quadratic minimum, the
wavenumber for maximum growth rate typically depends only linearly on~$\mu$, so
with $\mu \ll 1$, this lack of dependence on~$\mu$ in the SH~equation is
reasonable. However, with a quartic minimum, \citet{Proctor1991} argues that
additional terms should be included in the unfolding of a quartic minimum in
order to allow the wavenumber of maximum growth rate to depend on~$\mu$. For
this reason,  we propose the modified linear operator
 \begin{equation}
 \label{eq:PDE_linear_Proctor}
 \frac{\partial u}{\partial t} =
   \mu \left( 1 + \mu_{p} \left( 1 + \frac{\partial^{2} }{\partial x^{2}} \right) \right) u - \left( 1 + \frac{\partial^{2}}{\partial x^{2}} \right)^{4} u + f_{1} \frac{\partial^{2} u}{\partial x^{2}} - f_{2} \frac{\partial^{4} u}{\partial x^{4}},
 \end{equation}
where we will call the extra term, proportional to~$\mu_{p}$, the Proctor term,
although the form of this extra term differs from that proposed by
\citet{Proctor1991}.
The dispersion relation is now
 \[
 \sigma = \mu \left( 1 + \mu_{p} \left( 1 - k^{2} \right) \right) - \left( 1 - k^{2} \right)^{4} - f_{1} k^{2} - f_{2} k^{4},
 \]
and so
 \[
 \frac{d \sigma}{d k^{2}} = - \mu \mu_{p} + 4 \left( 1 - k^{2} \right)^{3} - f_{1} - 2 f_{2} k^{2}.
 \]
The curve of maximum (or minimum) growth rate is then defined by
$\frac{d\sigma}{dk^{2}}=0$, \textit{i.e.},
 \begin{equation*}

 \mu = \frac{1}{\mu_{p}} \left( 4 \left( 1 - k^{2} \right)^{3} - f_{1} - 2 f_{2} k^{2} \right).
 \end{equation*}
For $\mu_{p}\neq0$ and $f_1=f_2=0$ (Figure~\ref{fig:Proctor_curves}$a$),
this curve of maximum growth rate is a cubic and so is tangent to the neutral
stability curve, while for non-zero $f_1$ and $f_2$
(Figure~\ref{fig:Proctor_curves}$b$), the curves of maximum growth rate
intersect the two minima in the marginal stability curve transversally.

 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_02
  \begin{tikzpicture}
   \begin{axis}[width=6.5cm,height=5cm,
                scaled ticks=false,
                xmin=0.8,xmax=1.2,
                xtick={0.8,1,1.2},
                ymin=-0.0001,ymax=0.002,
                ytick={0,0.001,0.002},
                yticklabel style={anchor=east,
                                  /pgf/number format/precision=3,
                                  /pgf/number format/fixed,
                                  /pgf/number format/fixed zerofill},
               ]
    \addplot[color=blue, mark=none, line width=1pt]
             table [x index=0,y index=1] {paper-proctor.dat};
    \addplot[color=blue, mark=none, line width=1pt, dashed]
             table [x index=0,y index=2] {paper-proctor.dat};
   \end{axis}
\draw (4.0, 0.0) node[below] {$k^2$};
\draw (0, 2.75) node[left] {$\mu$};
\draw (0.25, 3.5) node[above] {(a)};
  \end{tikzpicture}
\hskip 0.5truecm
  \begin{tikzpicture}
   \begin{axis}[width=6.5cm,height=5cm,
                scaled ticks=false,
                xmin=0.8,xmax=1.2,
                xtick={0.8,1,1.2},
                ymin=0.0259,ymax=0.028,
                ytick={0.026,0.027,0.028},
                yticklabel style={anchor=east,
                                  /pgf/number format/precision=3,
                                  /pgf/number format/fixed,
                                  /pgf/number format/fixed zerofill},
               ]
    \addplot[color=blue, mark=none, line width=1pt]
             table [x index=0,y index=3] {paper-proctor.dat};
    \addplot[color=blue, mark=none, line width=1pt, dashed]
             table [x index=0,y index=4] {paper-proctor.dat};
   \end{axis}
\draw (4.0, 0.0) node[below] {$k^2$};
\draw (0, 2.75) node[left] {$\mu$};
\draw (0.25, 3.5) node[above] {(b)};
  \end{tikzpicture}
 \endpgfgraphicnamed}
 \end{center}
\caption{Marginal stability curve (solid) and curve of maximum (or minimum)
 growth (dashed) for~\eqref{eq:PDE_linear_Proctor}, for
 (a)~$f_{1} = f_{2} = 0$, and
 (b)~$f_{1}=0.05$ and $f_{2}=-0.0235$.
 The coefficient $\mu_{p} = -0.1$ in both cases.
 Two of the three dashed lines in~(b) join together above the top of the frame.
 \label{fig:Proctor_curves}}
\end{figure}

Including the Proctor term influences the shape of the marginal stability
curve. Throughout most of this paper we will assume that $\mu_{p}$ is small
enough not to influence the behaviour of solutions, apart from the asymptotic
analysis in~$\S$\ref{subsec:Proctor_term} (where the Proctor term necessarily
appears as a higher order term), and a numerical consideration
in~$\S$\ref{subsec:snaking}.

The linear part of the model derived so far is symmetric under spatial
reversals, \textit{i.e.}, $x$~and $-x$ are equivalent. Some systems break
reflection symmetry, for example, the Taylor--Couette system in the presence of
an azimuthal and axial magnetic field \citep{Stefani2009,Mamatsashvili2019};
this can be modelled by including terms such as $\left(f_{d0} -
f_{d1}\partial_{xx}\right)u_x$, leading to drifting solutions. However, the
reflection symmetry is useful in the analysis of the model equation, so we do
not include drift terms here.

\subsection{Nonlinear terms}
\label{subsec:nonlinear_terms}

Nonlinear terms in the model will saturate exponentially growing solutions at
finite amplitude, and should respect the symmetry (or lack thereof) of any
underlying physical systems. The usual SH~nonlinear term
is~$-u^3$, allowing $u\rightarrow-u$ symmetry, but here we break this symmetry
and chose
 \[
 n_{2} u^{2} + n_{3} u^{3}.
 \]
as nonlinear terms for the model. There are many more possibilities for
nonlinear terms that have been used by many authors
\citep{Knobloch1990d,Crawford1999b,Burke2012} for similar equations in
different contexts. For the purposes of this paper, we consider
only the nonlinear terms above, retaining the coefficients $n_2$ and $n_3$ as
parameters. We focus on the case where the bifurcation is supercritical, which
limits the range of values that $n_{2}$ and $n_{3}$ can take. In particular, we
must have
 \begin{equation}
 \label{eq:nA}
 n_A \equiv 3 n_{3} + \tfrac{326}{81} n_{2}^{2} < 0,
 \end{equation}
derived in section~\ref{subsec:derivation_GLE}. For numerical examples, we will
take $n_{2} = 0.1$ and $n_{3} = -1$ throughout.

We now have our complete model equation, written here with the Proctor term
 \begin{equation}
 \label{eq:PDE_nonlinear_Proctor}
 \frac{\partial u}{\partial t} =
   \mu \left( 1 + \mu_{p} \left( 1 + \frac{\partial^{2} }{\partial x^{2}} \right) \right) u - \left( 1 + \frac{\partial^{2}}{\partial x^{2}} \right)^{4} u + f_{1} \frac{\partial^{2} u}{\partial x^{2}} - f_{2} \frac{\partial^{4} u}{\partial x^{4}},
   + n_{2} u^{2} + n_{3} u^{3},
 \end{equation}
and without
 \begin{equation}
 \label{eq:PDE_nonlinear}
 \frac{\partial u}{\partial t} = \mu u - \left( 1 + \frac{\partial^{2}}{\partial x^{2}} \right)^{4} u + f_{1} \frac{\partial^{2} u}{\partial x^{2}} - f_{2} \frac{\partial^{4} u}{\partial x^{4}} + n_{2} u^{2} + n_{3} u^{3}.
 \end{equation}
For the remainder of this paper, we will concentrate mostly on the parameter
values where the marginal stability curves have two minima and there is
bistability between patterns of different wavelength, leading to the
possibility of these patterns coexisting in separate parts of the domain.
Although these equations are related to the two length scale models of
\citet{Lifshitz1997} and \citet{Rucklidge2012}, our model cannot be derived
from these simply by setting the two length scales to be equal.

\section{Weakly nonlinear analysis}
\label{sec:wnlt}

In this section we compute weakly nonlinear solutions for the
model~\eqref{eq:PDE_nonlinear} by deriving a generalised version of the
Ginzburg--Landau equation, and use it to establish where one-dimensional
periodic patterns are stable to long-wave perturbations. For most of the
calculations we set $\mu_{p}=0$, but we do consider the effect of the Proctor
term briefly in~$\S$\ref{subsec:Proctor_term}. We also compare our derivation
to previously published derivations for related problems
in~$\S$\ref{subsec:Kozyreff}.
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_03
 \mbox{\begin{tikzpicture}[xscale=4.0,yscale=4.5]
  \draw[->] (0,0) -- (2,0) node[below] {$k^{2}$};
  \draw[->] (0,0) -- (0,1.1) node[left] {$\mu$};
  \draw[<->] (0.43,0.34) -- (1.7,0.34);
  \draw[<->] (0.3,0.255) -- (0.3,0.335) node[left] {$\mathcal{O}(\epsilon^{2})$};
  \draw[<->] (1,0) -- (1,0.32);
  \draw (1,0.16) node[left] {$\mathcal{O}(\epsilon)$};
  \draw[line width=1pt,color=blue] plot[id=pde_scales2, domain=0.0:2, samples=50, smooth]
          function{1 - 4*x + 6*x**2 - 4*x**3 + x**4 + 0.7*x - 0.38*x**2};
  \draw (0.61,-1pt) -- (0.61,1pt);
  \draw (0.92,-1pt) -- (0.92,1pt);
  \draw (1.47,-1pt) -- (1.47,1pt);
  \draw (0.61,-1pt) node[below] {$k_{1}^{2}$};
  \draw (0.92,-1pt) node[below] {$k_{2}^{2}$};
  \draw (1.47,-1pt) node[below] {$k_{3}^{2}$};
  \draw (1,0.34) node[above] {$\mathcal{O}(\epsilon^{1/2})$};
\end{tikzpicture}}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:marginal_scalings} Scalings used for the asymptotic analysis.
 The deviations from a quartic minimum are contained in a
 $\mathcal{O}(\epsilon^{1/2}) \times \mathcal{O}(\epsilon^{2})$ box, and the
 whole curve is shifted by an $\mathcal{O}(\epsilon)$ amount.}
 \end{figure}

\subsection{Derivation of a generalised Ginzburg--Landau equation}
\label{subsec:derivation_GLE}

For the model PDE~\eqref{eq:PDE_nonlinear}, we consider only situations
where we have a small perturbation of the quartic marginal stability curve. At
the quartic minimum, an appropriate scaling is that if the solution $u$ is
of~$\mathcal{O}(\epsilon)$, then the slow time~$T$ and long length~$X$ should
be $\mathcal{O}(\epsilon^2)$ and $\mathcal{O}(\epsilon^{1/2})$ respectively:
the long length~$X$ is even longer than the $\mathcal{O}(\epsilon)$ scaling in
the case of a quadratic minimum in~\eqref{eq:GLE_scaling}. We therefore
restrict $f_{1}$ and~$f_{2}$ such that the two minima in the marginal stability
curve lie in an $\mathcal{O}(\epsilon^{1/2}) \times \mathcal{O}(\epsilon^{2})$
box, as illustrated in Figure~\ref{fig:marginal_scalings}. Other scalings are
possible.

To simplify the form of the marginal stability curve, we take $\mu_{p}=0$ and
consider the two minima and intermediate maximum of the function $\mu=p(k^2)$
in~\eqref{eq:polynomial}. The condition for an extremum is $p'(K)=0$, with
$K=k^2$: writing $\tfrac{1}{4}p'(K)$ as $(k^2-k_1^2)(k^2-k_2^2)(k^2-k_3^2)$,
with $k_1<k_2<k_3$, leads to the conclusion that
 \begin{equation}
 \label{eq:root_conditions}
 k_{1}^{2} + k_{2}^{2} + k_{3}^{2} = 3,
 \qquad
 k_{1}^{2} k_{2}^{2} + k_{1}^{2} k_{3}^{2} + k_{2}^{2} k_{3}^{2} = 3 + \tfrac{1}{2}f_{2},
 \qquad
 k_{1}^{2} k_{2}^{2} k_{3}^{2} = 1 - \tfrac{1}{4}f_{1}.
 \end{equation}
To satisfy these three equations,
we introduce new parameters $\gamma$ and~$\delta$, defined by
 \begin{equation}
 \label{eq:k_gamma_delta}
 k_{1}^{2} = 1 - \gamma - \tfrac{1}{2}\delta,
 \qquad
 k_{2}^{2} = 1 + 2 \gamma,
 \qquad
 k_{3}^{2} = 1 - \gamma + \tfrac{1}{2}\delta,
 \end{equation}
such that the first equation in~\eqref{eq:root_conditions} is satisfied.
Requiring $k_1<k_2<k_3$ implies \hbox{$\delta>6|\gamma|\geq0$}.
Substituting~\eqref{eq:k_gamma_delta} into the second and third equations
in~\eqref{eq:root_conditions} yields
 \begin{align*}
 k_{1}^{2} k_{2}^{2} + k_{1}^{2} k_{3}^{2} + k_{2}^{2} k_{3}^{2}
          &= 3 - 3\gamma^2 - \tfrac{1}{4}\delta^2
          = 3 + \tfrac{1}{2}f_{2}, \\
 k_{1}^{2} k_{2}^{2} k_{3}^{2}
          &= 1 - 3\gamma^2 - \tfrac{1}{4}\delta^2 - \tfrac{1}{2}\gamma\delta^2 +
             2\gamma^3
          = 1 - \tfrac{1}{4}f_{1}.
 \end{align*}
These are rearranged to result in expressions for $f_{1}$, $f_{2}$ and
$f_{1}+2f_{2}$ in terms of $\gamma$ and~$\delta$:
 \begin{equation*}
 f_{1} = 12\gamma^2 + \delta^2 + 2\gamma\delta^2 - 8\gamma^3,
 \qquad
 f_{2} = -6\gamma^2 - \tfrac{1}{2}\delta^2,
 \qquad
 f_{1}+2f_{2} = 2\gamma\left(\delta^2 - 4\gamma^2\right).
 \end{equation*}
The restriction $\delta>6|\gamma|$ means that $\delta^2 - 4\gamma^2>0$,
and so $f_{1}+2f_{2}=0$ implies $\gamma=0$ and $k_2=1$.

Now, if the minima and maxima of the marginal stability curve are to satisfy the
scalings in Figure~\ref{fig:marginal_scalings}, we need
$\delta=\mathcal{O}(\epsilon^{1/2})$ and $\gamma$~no larger than this. We also
need $p(k_{2}^2)-p(k_{1}^2)$ and $p(k_{2}^2)-p(k_{3}^2)$ to be
$\mathcal{O}(\epsilon^2)$:
 \begin{equation*}
 p(k_{2}^2)-p(k_{1}^2) = \tfrac{1}{16}
                         \left(\delta-2\gamma\right)
                         \left(\delta+6\gamma\right)^3,
 \qquad
 p(k_{2}^2)-p(k_{3}^2) = \tfrac{1}{16}
                         \left(\delta+2\gamma\right)
                         \left(\delta-6\gamma\right)^3,
 \end{equation*}
so with $\delta=\mathcal{O}(\epsilon^{1/2})$ and $\gamma$~no larger, these
differences fit within~$\mathcal{O}(\epsilon^2)$. The overall marginal stability
curve is shifted up by an amount
$p(1)\approx\tfrac{1}{2}\left(12\gamma^2+\delta^2\right)=\mathcal{O}(\epsilon)$,
again within the scalings indicated in Figure~\ref{fig:marginal_scalings}. The
scaling of $\gamma$ and $\delta$ imply $f_{1}$ and $f_{2}$ are both
$\mathcal{O}(\epsilon)$, but that $f_{1}+2f_{2}=\mathcal{O}(\epsilon^{3/2})$.

We are now in a position to perform the multiple scales analysis. The full set
of scalings used are
 \begin{align}
 \label{eq:wnlt_scaling}
 u &= \sum_{n=2}^{6} \epsilon^{n/2} u_{n/2},
 &
 \mu &= \epsilon \mu_{1} + \epsilon^{2} \mu_{2}, \nonumber \\
 \partial_{t} &= \epsilon^{2} \partial_{T},
 &
 f_{1} &= \epsilon f_{1},\\
 \partial_{x} &= \partial_{x} + \epsilon^{1/2} \partial_{X},
 &
 f_{2} &= \epsilon f_{2},
 \qquad\text{with $f_{1}+2f_{2}=\mathcal{O}(\epsilon^{3/2})$},
 \nonumber
 \end{align}
where $\epsilon\mu_{1}$ is the amount by which the marginal stability curve is
shifted. We also define a singular linear operator $\mathcal{L}$ to be
$(1+\partial_{x}^2)^4$, with $\mathcal{L}e^{\pm ix}=0$.

Inserting these scalings into~\eqref{eq:PDE_nonlinear}, we obtain at leading
order $\mathcal{O}(\epsilon)$:
 \[
 0 = - \mathcal{L} u_{1},
 \]
which is satisfied by taking
 \[
 u_{1} = A(X,T) e^{i x} + \text{c.c.},
 \]
where \hbox{c.c.} represents the complex conjugate. It is useful to observe that
$\partial_{x}^2u_1=-u_1$.

Proceeding to $\mathcal{O}(\epsilon^{3/2})$, we have
 \[
 0 = - \mathcal{L} u_{3/2} - 8 \partial_{xX} \underbrace{(1 + \partial_{x}^{2})^{3} u_{1}}_{=0}.
 \]
The most convenient way to solve this is to set $u_{3/2} = 0$.

At $\mathcal{O}(\epsilon^{2})$, we have
 \[
 0 = \mu_{1} u_{1} - \mathcal{L} u_{2}
    - 4 \partial_{X}^{2} \underbrace{\left( 1 + 9 \partial_{x}^{2} + 15 \partial_{x}^{4} + 7 \partial_{x}^{6} \right) u_{1}}_{=0}
    + f_{1} \partial_{x}^{2} u_{1}
    - f_{2} \partial_{x}^{4} u_{1}
    + n_{2} u_{1}^{2}.
 \]
The terms proportional to $e^{ix}$ have a prefactor of $\mu_{1}-\left(f_{1}+f_{2}\right)$
We need to eliminate these terms in order to solve for $u_{2}$, which we do by
setting $\mu_{1}  = f_{1} + f_{2}$. With this, the remaining terms are
$0=-\mathcal{L} u_{2}+n_{2} u_{1}^{2}$, which can be solved to give
  \[
  u_{2} = n_{2} \left( \tfrac{1}{81}A^{2} e^{2ix} + \tfrac{1}{81}\bar{A}^{2} e^{-2ix} + 2 |A|^{2} \right).
 \]
The factors $\tfrac{1}{81}$ come from $\mathcal{L}e^{2ix}=81e^{2ix}$ when
$\mathcal{L}$ is inverted.

At $\mathcal{O}(\epsilon^{5/2})$ we have
 \[
 0 = - \mathcal{L} u_{5/2}
     - 8 \partial_{xX} (1 + \partial_{x}^{2})^{3} u_{2}
     - 8 \partial_{x}^{3} \partial_{X} \underbrace{(3 + 10 \partial_{x}^{2} + 7 \partial_{x}^{4}) u_{1}}_{= 0}
     + \underbrace{2 f_{1} \partial_{x} \partial_{X} u_{1} - 4 f_{2} \partial_{x}^{3} \partial_{X} u_{1}}_{=2(f_{1}+2f_{2})\partial_{xX}u_1}.
 \]
The only terms involving $e^{\pm i x}$ on the \hbox{RHS} are
the two on the end that combine to have a prefactor of $f_{1}+2f_{2}$. As
discussed above, this combination is a factor of $\epsilon^{1/2}$ smaller than
$f_{1}$ and $f_{2}$ separately, so this term can be pushed
to~$\mathcal{O}(\epsilon^{3})$ and dropped from this equation. With this, the
two remaining terms in the equation can be solved for~$u_{5/2}$:
 \[
 u_{5/2} = \frac{32}{243}n_2\left(
                         i A A_{X} e^{2ix} - i \bar{A} \bar{A}_{X} e^{-2ix}
                         \right).
 \]
There is no constant term in~$u_{5/2}$.

Finally, continuing to $\mathcal{O}(\epsilon^{3})$, and including the term
pushed from $\mathcal{O}(\epsilon^{5/2})$, we have
\begin{align*}
 \partial_{T} u_{1} = \mu_{1} u_{2}
 &  + \mu_{2} u_{1} - \mathcal{L} u_{3}
    - 8 \partial_{xX} \left(1 + \partial_{x}^{2} \right)^{3} u_{5/2} \\
 &{}- 4 \partial_{X}^{2} \left( 1 + 9 \partial_{x}^{2} + 15 \partial_{x}^{4} + 7 \partial_{x}^{6} \right) u _{2}
    - 2 \partial_{X}^{4} \left( 3 + 30 \partial_{x}^{2} + 35 \partial_{x}^{4} \right) u_{1} \\
 &{}+ f_{1} \partial_{x}^{2} u_{2} - f_{2} \partial_{x}^{4} u_{2}
    + f_{1} \partial_{X}^{2} u_{1} - 6 f_{2} \partial_{x}^{2} \partial_{X}^{2} u_{1}
    + 2 n_{2} u_{1} u_{2} + n_{3} u_{1}^{3}\\
 &{}+ 2(f_{1}+2f_{2})\partial_{xX}u_1.
 \end{align*}
The solvability condition requires the elimination of all terms proportional to
$e^{\pm ix}$. There are no contributions from $\mathcal{L} u_{3}$ and from the
terms linear in~$u_2$ and~$u_{5/2}$. Recalling that $\partial_{x}^{2}u_1=-u_1$,
the terms proportional to~$e^{ix}$ result in the solvability condition:
 \begin{equation}
 \label{eq:generalised_GLE}
 \frac{\partial A}{\partial T} = \mu_{2} A
        + 2 i \nu_{1} \frac{\partial A}{\partial X}
        + \nu_{2} \frac{\partial^{2} A}{\partial X^{2}}
        - 16 \frac{\partial^{4} A}{\partial X^{4}}
        + \left( 3 n_{3} + \tfrac{326}{81} n_{2}^{2} \right) |A|^{2} A,
 \end{equation}
where $\nu_{1} = f_{1} + 2 f_{2}$ and $\nu_{2} = f_{1} + 6 f_{2}$.
The nonlinear term is $n_A|A|^2A$, with $n_A$ defined in~\eqref{eq:nA}.
This is a generalisation of the Ginzburg--Landau equation~\eqref{eq:GLE},
and each term in the equation is, in terms of the original unscaled variables,
of order~$\mathcal{O}(\epsilon^{3})$.
Similar amplitude equations have been proposed before
\citep{Riley1989,Riecke1990,Raitt1995}, although not formally derived via
an asymptotic expansion, to model other problems with very flat marginal
stability curves.

\subsection{Nonlinear stability of rolls}
\label{subsec:nonlinear_stability}

We can find roll solutions to~\eqref{eq:generalised_GLE} easily, and
use the equation to examine their stability to long wavelength Eckhaus
instabilities. We restrict ourselves to the supercritical case, where the
coefficient of the nonlinear term is negative. We consider a roll solution at
slightly off-critical wavenumber, \textit{i.e.},
 \begin{equation}
 \label{eq:GLE_offcritical_rolls}
  A(X,T) = R e^{i q X},
 \end{equation}
which corresponds to a solution $u = R \exp\left(i \left( 1 + \epsilon^{1/2} q
\right) x\right)$ of~\eqref{eq:PDE_nonlinear} since $X = \epsilon^{1/2} x$, and $q =
\mathcal{O}(1)$. These are also known as phase-winding solutions.
Substituting~\eqref{eq:GLE_offcritical_rolls} into~\eqref{eq:generalised_GLE},
and rearranging, we obtain
 \begin{equation}
 \label{eq:wnlt_roll_amplitude}
 R^{2}
       = - \frac{1}{n_A}\left( \mu_{2} - 2 \nu_{1} q - \nu_{2} q^{2} - 16 q^{4} \right),
 \end{equation}
where $n_A<0$ for a supercritical bifurcation.
The existence boundary of rolls, which is also the marginal stability curve,
is where $R^{2} = 0$, or equivalently
 \begin{equation}
 \label{eq:wnlt_marginal}
 \mu_{2} = 2 \nu_{1} q + \nu_{2} q^{2} + 16 q^{4},
 \end{equation}
with $\nu_1$ and $\nu_2$ here playing the roles of unfolding parameters
for the quartic minimum of the marginal stability curve.

To determine the stability of these roll solutions, we
perturb~\eqref{eq:GLE_offcritical_rolls}, writing
 \begin{equation*}

 A(X,T) = R \left( 1 + r(X,T) \right) e^{i \left( q X + \phi(X,T) \right)},
 \end{equation*}
where $|r|$,$|\phi| \ll 1$, following \citet{Hoyle2006}.
We substitute this expressions into~\eqref{eq:generalised_GLE},
linearise and separate the real and imaginary parts to obtain two linear
constant coefficient PDEs for $r$ and~$\phi$. To solve these, we seek solutions
of the form $e^{\sigma T + i m X}$, with $m\ll q$ and derive a quadratic
equation for the growth rate~$\sigma$. One root of this is always negative for a
supercritical solution, and the other root is
 \begin{equation*}
  \sigma = - m^{2} \left( \nu_{2} + 96q^{2} + \frac{2}{n_A R^{2}} \left( \nu_{1} + \nu_{2} q + 32 q^{3} \right)^{2} \right) + \mathcal{O}(m^{3}).
 \end{equation*}
See \citet{Bentley2012} for details.
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_04
 \begin{tikzpicture}[xscale=25, yscale=350]
  \draw[line width=1pt,color=blue] plot[id=fig6_eck, domain=-0.2:0.2, samples=100, smooth]
          function{-0.28*x**2 + 16*x**4};
  \draw[line width=1pt,color=blue, dashed] plot[id=fig6_eck1, domain=-0.18:-0.055, samples=100, smooth]
          function{-0.28*x**2 + 16*x**4 + (-0.28*x + 32*x**3)**2 / (-0.28 + 96*x**2)};
  \draw[line width=1pt,color=blue, dashed] plot[id=fig6_eck2, domain=0.055:0.18, samples=100, smooth]
          function{-0.28*x**2 + 16*x**4 + (-0.28*x + 32*x**3)**2 / (-0.28 + 96*x**2)};
\draw[->] (-0.21,-0.0015) -- (0.2,-0.0015);
\draw[->] (-0.21,-0.0015) -- (-0.21,0.016);
\draw (0.2,-0.0015) node[below] {$q$};
\draw[-] (0,-0.0017) -- (0,-0.0013);
\draw (0,-0.0015) node[below] {$q=0$};
\draw (-0.21,0.016) node[left] {$\mu_{2}$};
\end{tikzpicture}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:marg_and_eck} Marginal stability curve (solid) and Eckhaus
 curve (dashed) for $\nu_{1} = 0$ and $ \nu_{2} = -0.28$. The regions above the
 dashed curve are Eckhaus stable.}
\end{figure}

For the rolls to be unstable, we require $\sigma > 0$ for some $m$,
giving a stability boundary defined by
 \[
 R^{2} = - 2 \frac{\left( \nu_{1} + \nu_{2} q + 32 q^{3} \right)^{2}}{n_A \left( \nu_{2} + 96 q^{2} \right)},
 \]
and, recalling~\eqref{eq:wnlt_roll_amplitude}, we find the Eckhaus stability
boundary
 \begin{equation}
 \label{eq:wnlt_eckhaus}
 \mu_{2} = 2 \nu_{1} q + \nu_{2} q^{2} + 16 q^{4} + \frac{2}{\nu_{2} + 96 q^{2}} \left( \nu_{1} + \nu_{2} q + 32 q^{3} \right)^{2}.
 \end{equation}
Figure~\ref{fig:marg_and_eck} shows an example of the
marginal~\eqref{eq:wnlt_marginal} and Eckhaus~\eqref{eq:wnlt_eckhaus} curves. Inside the
Eckhaus (dashed) curve, patterns are stable to long wavelength disturbances.
Note that $\nu_{2} < 0$ when there are two distinct minima, so~\eqref{eq:wnlt_eckhaus}
has two vertical asymptotes at $q = \pm \sqrt{-\nu_{2}/96}$.
Between these asymptotes, $\mu_2$~is below the marginal stability curve.

\subsection{The Proctor term}
\label{subsec:Proctor_term}

The presence of these asymptotes means that the Eckhaus boundary does not close
up in the middle, as might be expected. However, \citet{Proctor1991} considered
an amplitude equation similar to~\eqref{eq:generalised_GLE} but with the first
term on the RHS replaced by $\mu_2(1\pm{i}\partial_X)A$, and found that
(depending on parameters) the inner edges of the left and right stability
boundaries can meet, or the stability region can close with increasing~$\mu$
with two separate stability regions, and the stability boundaries can be
non-monotonic -- see \citet{Proctor1991} and \citet{Bentley2012} for examples.

However, when we considered the model
equation~\eqref{eq:PDE_nonlinear_Proctor}, with the Proctor term included, we
encountered difficulties. Writing $\mu=\epsilon\mu_1+\epsilon^2\mu_2$, as
in~\eqref{eq:wnlt_scaling}, turns out to be inadequate, so we modified the
scaling to include an additional
$\epsilon^{3/2}\mu_{3/2}\left(1+\partial_{x}^{2}\right)$ term, aiming to relate
$\mu_{3/2}$ to~$\mu_{p}$. It turns out that this doesn't work either, and
possibly further (potentially higher order) terms would need to be considered
for a consistent scaling \citep{Bentley2012}.

\subsection{The Lugiato--Lefever and complex Swift--Hohenberg equations}
\label{subsec:Kozyreff}

As an aside, we consider the
Lugiato--Lefever equation
 \begin{equation}
 \label{eq:lugato_lefever}
 \frac{\partial \psi}{\partial t} = S - \left( 1 + i \bigtriangleup \right) \psi + i | \psi |^{2} \psi - i B_{2} \frac{\partial^{2} \psi}{\partial \tau^{2}} + i B_{4} \frac{\partial^{4} \psi}{\partial \tau^{4}},
 \end{equation}
which governs the envelope of the complex electromagnetic field~$\psi(t,\tau)$
inside a photonic crystal fibre cavity. In this equation, $S$~represents an
injected field, $\bigtriangleup$~is a cavity detuning, and $B_{2}$ and $B_{4}$
incorporate chromatic dispersion. The two time variables $t$ and $\tau$
represent, respectively, the average evolution of~$\psi$ over one cavity round
trip and the fast variations of~$\psi$. We note the similarities
between~\eqref{eq:lugato_lefever} and~\eqref{eq:PDE_nonlinear}: the independent
variable~$\tau$ is equivalent to~$x$, there is a cubic nonlinearity, the
inhomogeneous $S$ term breaks the $\psi \rightarrow -\psi$ symmetry and implies
a quadratic nonlinearity, and a fourth order complex equation is equivalent to
an eighth order real equation.

This system also allows marginal stability curves with a double minimum,
and \citet{Kozyreff2009} derived a similar amplitude equation for this case.
Their
analysis concentrates only on the case when the two minima
occur at the same height, equivalent to $\gamma=0$ in~\eqref{eq:k_gamma_delta}.
This is a degenerate situation however; to recover the
generic situation it is necessary to include a $B_{3} \partial^{3} \psi /
\partial \tau^{3}$ term in~\eqref{eq:lugato_lefever}. The degenerate case is
chosen by \citet{Kozyreff2009} both as a means of simplifying the analysis and as a
situation easily achievable experimentally.

With $\gamma=0$, \citet{Bentley2012} showed that the appropriate amplitude
equation for the Lugiato--Lefever equation is~\eqref{eq:generalised_GLE} but with
$\nu_1=0$. This special case is interesting: with $\nu_1=0$ and for $\nu_2<0$,
equation~\eqref{eq:generalised_GLE} is equivalent to the \emph{complex} Swift--Hohenberg
equation:
 \begin{equation}
 \label{eq:complex_SHE}
 \frac{\partial A}{\partial T} = \lambda A - \left( 1 + \frac{\partial^{2}}{\partial X^{2}} \right)^{2} A + n_{A} |A|^{2} A,
 \end{equation}
having scaled and changed notation. This differs from the \emph{real}
SH~equation~\eqref{eq:SHE} in that there is no quadratic term
and that the cubic nonlinearity is $|A|^2A$ rather than~$u^3$.
It also differs from the more usual complex SH equation,
which has complex coefficients \citep{Sakaguchi1997}.
The complex SH~equation with real coefficients has been investigated by
\citet{Gelens2010} and is equivalent to to the equation studied by
\citet{Raitt1995}.
However, \citet{Kozyreff2009} find the real SH~equation~\eqref{eq:SHE} as the
amplitude equation for the Lugiato--Lefever problem, rather than the complex
Swift--Hohenberg equation with real coefficients~\eqref{eq:complex_SHE}. We
believe this is an error, and it means that the interpretation by
\citet{Kozyreff2009} of localised solutions of~\eqref{eq:lugato_lefever} as
localised solution of~\eqref{eq:SHE} is incorrect; instead they should be
interpreted in terms of localised solutions of~\eqref{eq:complex_SHE}, which are
different. We investigate these localised solutions in more detail below.

\section{Lyapunov functional and the first integral}
\label{sec:Lyapunov_first_int}

In the Swift--Hohenberg equation, the Lyapunov functional and first integral
are useful for finding localised solutions: if two solutions are to be
connected by a stationary front, they should have the same values of the first
integral and similar values of the Lyapunov functional. In this section, we
generalise the SH results to the model equation~\eqref{eq:PDE_nonlinear}.

Multiplying the steady version of~\eqref{eq:PDE_nonlinear} by $u_x$ and
integrating by parts gives us a first integral:
 \begin{align*}

 \mathcal{H} &= - \tfrac{1}{2} \left( \mu - 1 \right) u^{2}
                - \tfrac{1}{3}n_{2}u^{3}
                - \tfrac{1}{4}n_{3}u^{4}
                + \tfrac{1}{2} \left( 4 - f_{1} \right) u_{x}^{2}
                - \tfrac{1}{2} \left( 6 + f_{2} \right) \left( u_{xx}^{2} - 2 u_{x} u_{xxx} \right) \\ \nonumber
    & \qquad {} + 2 u_{xxx}^{2}
                - 4 u_{xx} u_{xxxx}
                + 4 u_{x} u_{xxxxx}
                - \tfrac{1}{2} u_{xxxx}^{2}
                + u_{xxx} u_{xxxxx}
                - u_{xx} u_{xxxxxx}
                + u_{x} u_{xxxxxxx}.
 \end{align*}
Any steady solution of~\eqref{eq:PDE_nonlinear} must have
$\frac{d\mathcal{H}}{dx}=0$. It is a straight-forward modification to include
the Proctor term from~\eqref{eq:PDE_nonlinear_Proctor}.

Using the Lyapunov functional for the Swift--Hohenberg
equation~\eqref{eq:SH_Lyaponov} as a starting point, we define a similar
functional for~\eqref{eq:PDE_nonlinear}, namely
 \begin{equation*}

 \mathcal{F}[u] = \int_{0}^{L} \left(
                       - \tfrac{1}{2} \mu u^{2}
                       + \tfrac{1}{2} \left( \left( 1 + \partial_{x}^{2} \right)^{2} u \right)^{2}
                       + \tfrac{1}{2} f_{1} \left( \partial_{x} u \right)^{2}
                       + \tfrac{1}{2} f_{2} \left( \partial_{x}^{2} u \right)^{2}
                       - \tfrac{1}{3} n_{2}u^{3}
                       - \tfrac{1}{4} n_{3}u^{4} \right) dx.
 \end{equation*}
One can readily show that
 \[
 \frac{\partial u}{\partial t} = - \frac{\delta \mathcal{F}}{\delta u}
 \qquad \text{and} \qquad
 \frac{d \mathcal{F}}{d t} = - \int_{0}^{L} \left( \frac{\partial u}{\partial t} \right)^{2} dx \leq 0.
 \]
It is possible to show that $\mathcal{F}[u]$ is bounded below provided $f_2>-4$
and~$n_3<0$
\citep{Bentley2012},
and so, as in the SH equation, that stable solutions are local minima
of~$\mathcal{F}[u]$. Similarly, it is a straight-forward modification to
include the Proctor term from~\eqref{eq:PDE_nonlinear_Proctor}.

Small-amplitude solutions of~\eqref{eq:PDE_nonlinear} can be found using the
generalised GL equation~\eqref{eq:generalised_GLE}, and these can be used to
find weakly nonlinear estimates of $\mathcal{H}$ and~$\mathcal{F}$.
Alternatively, fully nonlinear solutions of~\eqref{eq:PDE_nonlinear}  can be
found numerically, for example by using AUTO \citep{Doedel2007a}. Examples of
$\mathcal{H}$ and~$\mathcal{F}$ computed numerically in this way are shown in
Figure~\ref{fig:H_and_F_vs_k}, in the cases where the minima in the marginal
stability curve are at the same height, and where there is a single minimum
with an almost-minimum just outside the cusp. These were computed using an
initial solution on a domain of size $L=2\pi$, which was then continued in~$L$,
increasing and decreasing to cover the range of wavenumbers for which a pattern
solution exists. We note that $\mathcal{H}$ and~$\mathcal{F}$ are both zero at
the extremities of the existence region, and that the extrema of~$\mathcal{H}$
correspond to the Eckhaus stability boundaries.
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_05ab
 \begin{tikzpicture}
 \node [inner sep=0pt,outer sep=0pt] at (0,0) {\includegraphics[trim=5mm 15mm 10mm 10mm, clip, width=6.0cm, height=4cm]{ch6_1_hm1.pdf}};
\draw[->] (-3.0,-2.0) -- (3.0,-2.0);
\draw[->] (-3.0,-2.0) -- (-3.0,2.0);
\draw (0,2.0) node(above) {$(a)$};
\draw (-3.0,2.0) node[left] {$\mathcal{H}$};
\draw ( 3.0,-2.0) node[below] {$k$};
\draw[dashed] (-3.0,0.45) -- (3.0,0.45);
\draw (-3.0,0.45) node[left] {$0$};
  \end{tikzpicture}
  \begin{tikzpicture}
  \node [inner sep=0pt,outer sep=0pt] at (0,0) {\includegraphics[trim=5mm 15mm 10mm 10mm, clip, width=6.0cm, height=4cm]{ch6_1_lf1.pdf}};
\draw[->] (-3.0,-2.0) -- (3.0,-2.0);
\draw[->] (-3.0,-2.0) -- (-3.0,2.0);
\draw (0,2.0) node(above) {$(b)$};
\draw (-3.0,2.0) node[left] {$\mathcal{F}$};
\draw ( 3.0,-2.0) node[below] {$k$};
\draw[dashed] (-3.0,0.45) -- (3.0,0.45);
\draw (-3.0,0.45) node[left] {$0$};
\end{tikzpicture}
\endpgfgraphicnamed}%
 \end{center}
  \vspace{-5ex}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_05cd
 \begin{tikzpicture}
 \node [inner sep=0pt,outer sep=0pt] at (0,0) {\includegraphics[trim=5mm 15mm 10mm 5mm, clip, width=6.0cm, height=4cm]{ch6_3_hm1.pdf}};

\draw[->] (-3.0,-2.0) -- (3.0,-2.0);
\draw[->] (-3.0,-2.0) -- (-3.0,2.0);
\draw (0,2.0) node(above) {$(c)$};
\draw (-3.0,2.0) node[left] {$\mathcal{H}$};
\draw ( 3.0,-2.0) node[below] {$k$};
\draw[dashed] (-3.0,1.30) -- (3.0,1.30);
\draw (-3.0,1.30) node[left] {$0$};
  \end{tikzpicture}
  \begin{tikzpicture}
  \node [inner sep=0pt,outer sep=0pt] at (0,0) {\includegraphics[trim=5mm 15mm 10mm 5mm, clip, width=6.0cm, height=4cm]{ch6_3_lf1.pdf}};

\draw[->] (-3.0,-2.0) -- (3.0,-2.0);
\draw[->] (-3.0,-2.0) -- (-3.0,2.0);
\draw (0,2.0) node(above) {$(d)$};
\draw (-3.0,2.0) node[left] {$\mathcal{F}$};
\draw ( 3.0,-2.0) node[below] {$k$};
\draw[dashed] (-3.0,1.30) -- (3.0,1.30);
\draw (-3.0,1.30) node[left] {$0$};
\end{tikzpicture}
\endpgfgraphicnamed}%
 \end{center}
  \vspace{-5ex}
 \caption{\label{fig:H_and_F_vs_k} Plot of $(a)$ $30 \mathcal{H}$ and $(b)$ $3
 \times 10^{3} \mathcal{F}$, against $k$, for $f_{1}=0.14$, $f_{2}=-0.07$ and
 $\mu=0.07$, computed using AUTO.
 $(c,d)$: $15 \mathcal{H}$ and $0.5\times10^{3}\mathcal{F}$,
 for $f_{1}=0.0488$, $f_{2}=-0.0227$ and $\mu=0.0306$.
 The red crosses correspond to Eckhaus unstable
 wavenumbers, and the blue circles to Eckhaus stable wavenumbers. The marginal
 stability curve (solid black) and Eckhaus curves (black crosses) are also shown
 -- the change from red crosses to blue circles does not exactly match the
 Eckhaus boundary owing to the scaling of $\mathcal{H}$ and~$\mathcal{F}$.
 Note that in $(c,d)$ there is only one minimum in the marginal stability
 curve, but there is a region of Eckhaus-stable patterns above the
 almost-minimum.}
 \end{figure}

The third tool, mentioned in $\S$\ref{sec:introduction}, is the normal form of
this variant of the Hamiltonian--Hopf bifurcation with four-fold degenerate
eigenvalues~$\pm{i}$. We derive the normal form for this bifurcation in
Appendix~\ref{sec:appendix} (see equations~\eqref{eq:app:normal_form_version_1}
and~\eqref{eq:app:normal_form_version_2}), but as we haven't found any first integrals of the
normal form, we don't see a way to use it at this point.

\section{Localised solutions}
\label{sec:localised}

The first integral~$\mathcal{H}$ and the Lyapunov functional~$\mathcal{F}$ are
two imporant tools for identifying where a pattern of one type can be localised
within a background of a pattern of another type: the values of the first
integral for the two patterns must be the same (since
$\frac{d\mathcal{H}}{dx}=0$ on any steady solution) and the Lyapunov functional
for the two patterns should be approximately the same.

We see from the example in Figure~\ref{fig:H_and_F_vs_k}$(a,b)$ that there is a
range of possible wavenumbers which satisfy the requisite criteria for patterns
with two different wavenumbers coexisting, namely there are wavenumbers that
have the same value of $\mathcal{H}$, and there are (different) wavenumbers with
the same value of $\mathcal{F}$. To narrow down the allowable wavenumbers, we
look for wavenumber pairs $(k_{-} < 1, k_{+} > 1)$ such that $\mathcal{H}(k_{-})
= \mathcal{H}(k_{+})$ and $\mathcal{F}(k_{-}) = \mathcal{F}(k_{+})$. We do this
by looking for intersections of contour lines plotted in the $(k_{+},k_{-})$
plane. For the parameter values in Figure~\ref{fig:H_and_F_vs_k}$(a,b)$, a pair of
wavenumbers that satisfy this condition are $(k_{+},k_{-})=(1.0778,0.8839)$. We
view such a point as an extension of the Maxwell point for the Swift--Hohenberg
equation, though it plays a different role: in the Swift--Hohenberg the
localised solutions are organised about the Maxwell point, whereas we use the
extension merely as a starting point to look for localised solutions.
 %
 \begin{figure}
 \begin{center}
\begin{tabular}{c}
 \mbox{\beginpgfgraphicnamed{paper_figure_06a
  \begin{tikzpicture}

\begin{axis}[width=14cm,height=4.5cm,xtick={0,400},ytick={-0.04,0.04},xmin=0, xmax=401.7312]
 \addplot[mark=none, color=blue] table[x index=0, y index=1]{\tbl};
\end{axis}
\draw (6.2,0) node[below] {$x$};
\draw (0,1.5) node[left] {\hspace{1.5em} $u$};
\draw (0.25,0) node[above] {$(a)$};
\end{tikzpicture}
\endpgfgraphicnamed}%
\\
 \mbox{\beginpgfgraphicnamed{paper_figure_06b
  \begin{tikzpicture}

\begin{axis}[width=14cm,height=4.5cm,xtick={0,400},ytick={-0.04,0.04},scaled ticks=false,ymin=0.8, ymax=1.2, xmin=0, xmax=401.7312]
 \addplot[color=blue,mark=o,only marks,mark size=0.5pt] table[x index=0, y index=2]{\tbl};
\addplot[mark=none,dashed] coordinates{(0,0.9015) (401.7312,0.9015)};
\addplot[mark=none,dashed] coordinates{(0,1.0895) (401.7312,1.0895)};
\end{axis}
\draw (6.2,0) node[below] {$x$};
\draw (0,0.75) node[left] {\hspace{1em} $k_{-}$};
\draw (0,2.1) node[left] {$k_{+}$};
\draw (0.25,0) node[above] {$(b)$};
\end{tikzpicture}
\endpgfgraphicnamed}%
  \end{tabular}
 \end{center}
  \vspace{-2ex}
\caption{\label{fig:PDE_example} Numerical simulation of mixed pattern initial
condition with wavenumbers $k=58/64=0.90625$ and $k=70/64=1.09375$. The
parameter values are: $f_{1}=0.14$, $f_{2}=-0.07$ and $\mu = 0.07$,
$n_{2}=0.1$, $n_{3}=-1$, $L=64\times2\pi$ and timestep~$0.01$. $(a)$~final
solution profile $u(x)$. $(b)$~approximation to the local wavenumber, defined
in \eqref{eq:local_wavenumber}. The values indicated are $k_{-}=0.9015$ and
$k_{+}=1.0895$.}
 \end{figure}

On a periodic domain of length~$L$, wavenumbers are restricted to  integer
multiples of $k=2\pi/L$. We therefore construct an initial condition consisting
of a region of pattern with wavenumber close to~$1.0778$ embedded in a
background of pattern with wavenumber close to~$0.8839$. Fixing a domain size
$L=64\times2\pi$, we choose $k_{-} = 58/64$ and $k_{+} = 70/64$, and
solve~\eqref{eq:PDE_nonlinear} using a second-order numerical scheme based on
Exponential Time Differencing \citep{Cox2002}. One example solution after
transients can be seen in Figure~\ref{fig:PDE_example}, which demonstrates that
localised solutions to~\eqref{eq:PDE_nonlinear} exist and are stable. This
solution is made up of a high wavenumber patch ($k\approx1.0895$) in the
centre of the domain, surrounded by low wavenumber regions ($k\approx0.9015$).
The approximation to the local wavenumber in Figure~\ref{fig:PDE_example}$(b)$
is found via
 \begin{equation}
 \label{eq:local_wavenumber}
 \textrm{local wavenumber} = \sqrt{\frac{- u_{xx}}{u}},
 \end{equation}
and matches (at least approximately) the expected values. The oscillations seen
in the amplitude and in the local wavenumber represent a beating between the
two constituent wavenumbers $k_{-}$ and~$k_{+}$ as the pattern adjusts from one
wavenumber to the other.

We can continue this and other solutions we have found in AUTO, continuing
in~$\mu$ to obtain the bifurcation diagram shown in Figure~\ref{fig:bifurcation_diag}.
Localised solutions lie on distinct branches that do not join up, created in
saddle-node bifurcations. The solution branches corresponding to a pattern of
single wavelength $k_{-}$ or $k_{+}$ are included for reference. On localised
solutions branches closer to the $k_{-}$ branch, more of the domain is filled by
the $k_{-}$ pattern than the $k_{+}$ pattern, and \textit{vice versa}.
The localised solution branches extend to values of $\mu$ below the value at the
local maximum of the marginal stability curve~($\mu(k_2)=0.0699$).
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_07
\begin{tikzpicture}[spy using outlines={rectangle, magnification=4,connect spies}]
\begin{axis}[axis lines=none,xtick=\empty,ytick=\empty,scaled ticks=false,width=10cm,height=6cm,xmin=-0.0005,xmax=0.01,ymin=0,ymax=0.35]
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b0.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=black,line width=1.0pt] table[x index=0,y index=1]{b1.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=black,line width=1.0pt] table[x index=0,y index=1]{b14.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b2.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b3.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b4.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b5.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b6.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b7.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b8.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b9.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b10.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b11.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b12.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b13.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b15.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b16.dat};
 \addplot[restrict x to domain=0:0.01,mark=none,color=blue] table[x index=0,y index=1]{b17.dat};
\coordinate (spypoint) at (axis cs:0.00097,0.078);
\coordinate (spyviewer) at (axis cs:0.0015,0.26);
\spy[width=2cm,height=2cm] on (spypoint) in node [fill=white] at (spyviewer);
\end{axis}
 \draw[->] (0.375,0) -- (0.375,4.5);
 \draw[->] (0.375,0) -- (9.5,0);
 \draw[dashed] (1.55,0) -- (1.55,1.6);
 \draw (9.5,0) node[below] {$\mu$};
 \draw (0.375,0) node[below] {$\mu_{c}$};
 \draw (0.375,0) node[left] {$0$};
 \draw (8.375,0) node[below] {$0.0787$};
 \draw (8.375,-0.1) -- (8.375,0.1);
 \draw (0.375,4.5) node[left] {$L_{2}$};
 \draw (8.375,4.2) node[right] {$u(k_{+})$};
 \draw (8.375,2.15) node[right] {$u(k_{-})$};
 \draw (1.55,0) node[below] {$\mu(k_{2})$};
 \end{tikzpicture}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:bifurcation_diag} Bifurcation diagram for parameter values
 $f_{1}=0.14$, $f_{2}=-0.07$. The critical value of the bifurcation parameter
 $\mu$ is $\mu_{c}=0.0687$, and the local maximum of the marginal stability
 curve occurs at $(k_{2}, \mu(k_{2})) = (1,0.0699)$. The thick black branches
 correspond to periodic patterns with wavenumbers $k_{-}$, $k_{+}$, and are
 given for reference. The blue branches correspond to localised solutions. Inset
 is a magnification of the saddle-nodes.}
 \end{figure}

We find similar disconnected branches of localised solutions even when the local
minima in the marginal stability curve are not at the same height, and even just
outside the cusp, where there is a single minimum and a second almost-minimum,
as in Figure~\ref{fig:H_and_F_vs_k}$(c,d)$. Solutions in this region rely on there
being a region of Eckhaus-stable patterns still present above the
almost-minimum, disconnected from the marginal stability curve.

\subsection{Interpretation of localised solutions via the amplitude equation}

The amplitude equation~\eqref{eq:generalised_GLE} has phase-winding solutions
$A=Re^{iqX}$, with $R$ and~$q$ related by~\eqref{eq:wnlt_roll_amplitude}.
Previous work on this amplitude equation with $\nu_1=0$ \citep{Raitt1995,
Raitt1997, Gelens2009, Gelens2010} -- the complex
Swift--Hohenberg equation with real coefficients -- has identified solutions
that are combinations of two phase-winding
solutions with positive and negative values of~$q$: these are precisely the
localised patterns we found in the model PDE~\eqref{eq:PDE_nonlinear} and shown in
Figure~\ref{fig:PDE_example} (with $\nu_1=f_{1}+2f_{2}=0$). Here we extend this
interpretation to the case~\hbox{$\nu_1\neq0$}.

In order to find localised solutions, we could develop a first integral and a
Lyapunov function for~\eqref{eq:generalised_GLE} and look for pairs of solutions with the
same values of the quantities. We reserve this for future work, and instead
locate localised solutions of~\eqref{eq:generalised_GLE} by starting with a mixture of
two phase-winding solutions with constituent wavenumbers $q_{-}$ and~$q_{+}$,
and timestepping the~\hbox{PDE}. Two examples are shown in
Figures~\ref{fig:EGL_example_1} (with $\nu_1\neq0$ and a small value of~$\nu_2$)
and~\ref{fig:EGL_example_2} (with $\nu_1=0$ and a larger value of~$\nu_2$).
 %
 \begin{figure}
 \begin{center}
  \mbox{\beginpgfgraphicnamed{paper_figure_08
  \begin{tabular}{cc}
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={-0.04,0.04},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=1]{ch6_amp_004.dat};
    \addplot[color=blue,dashed,mark=none] table[x index=0,y index=2]{ch6_amp_004.dat};
   \end{axis}
   \draw (0.3,0) node[above] {$(a)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$A$};
  \end{tikzpicture}
 &
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={0.035,0.038},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=3]{ch6_amp_004.dat};
   \end{axis}
   \draw (0.3,0) node[above] {$(b)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$|A|$};
  \end{tikzpicture}
\\
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={-0.06,0.06},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=4]{ch6_amp_004.dat};

   \end{axis}
   \draw (0.3,0) node[above] {$(c)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$u$};
  \end{tikzpicture}
&
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={0.9,1.1},xmin=0,xmax=401.12,ymin=0.9,ymax=1.1]
    \addplot[color=blue,mark=o,only marks,mark size=0.5pt] table[x index=0,y index=5]{ch6_amp_004.dat};
    \draw[dashed] (axis cs:0,0.946) -- (axis cs:401.12,0.946);
    \draw[dashed] (axis cs:0,1.051) -- (axis cs:401.12,1.051);
   \end{axis}
   \draw (0.3,0) node[above] {$(d)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$k$};
   \draw (0,2.15) node[left] {$k_{+}$};
   \draw (0,0.65) node[left] {$k_{-}$};
  \end{tikzpicture}
  \end{tabular}
  \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:EGL_example_1} Solution to~\eqref{eq:generalised_GLE}, with $\nu_{1}=0.0034$, $\nu_{2}=-0.0874$, $\mu_{2}=0.004$, $n_{2}=0.1$ and $n_{3}=-1$. $(a)$: the real (solid) and imaginary (dashed) parts of the amplitude $A$, $(b)$: the absolute value $|A|$, $(c)$: the reconstructed solution $u = A e^{i x} + \bar{A} e^{-i x}$, and $(d)$: an approximation to the local wavenumber of the reconstructed solution in $(c)$. The dashed lines correspond to $(k_{-},k_{+}) = (0.946,1.053)$.}
 \end{figure}
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_09
 \begin{tabular}{c}
  \begin{tabular}{cc}
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={-0.04,0.04},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=1]{ch6_amp_0025.dat};
    \addplot[color=blue,dashed,mark=none] table[x index=0,y index=2]{ch6_amp_0025.dat};
   \end{axis}
   \draw (0.3,0) node[above] {$(a)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$A$};
  \end{tikzpicture}
&
  \begin{tikzpicture}
   \begin{axis}[width=7cm, height=4.5cm,xtick={0, 400},ytick={0.028,0.037},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=3]{ch6_amp_0025.dat};
   \end{axis}
   \draw (0.3,0) node[above] {$(b)$};
   \draw (2.8,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$|A|$};
  \end{tikzpicture}
  \end{tabular}
\\
  \begin{tikzpicture}
   \begin{axis}[width=14cm, height=4.5cm,xtick={0, 400},ytick={-0.07,0.07},xmin=0,xmax=401.12]
    \addplot[color=blue,mark=none] table[x index=0,y index=4]{ch6_amp_0025.dat};
    \addplot[color=black,mark=none] table[x index=0,y expr=2*\thisrowno{3}]{ch6_amp_0025.dat};
   \end{axis}
   \draw (0.3,0) node[above] {$(c)$};
   \draw (6.2,0) node[below] {$x$};
   \draw (0,1.55) node[left] {$u$};
  \end{tikzpicture}
  \end{tabular}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:EGL_example_2} Solution to~\eqref{eq:generalised_GLE}, with $\nu_{1}=0$, $\nu_{2}=-0.28$, $\mu_{2}=0.0025$, $n_{2}=0.1$ and $n_{3}=-1$. $(a)$: the real (solid) and imaginary (dashed) parts of the amplitude $A$, $(b)$: the absolute value $|A|$, and $(c)$: the reconstructed solution $u = A e^{i x} + \bar{A} e^{-i x}$, with the amplitude $2 |A|$ plotted also.}
\end{figure}
 %
Figure~\ref{fig:EGL_example_1}$(a)$ shows the real (solid) and imaginary (dashed)
parts of the solution $A$, $\Re(A)$ and $\Im(A)$ respectively. We can clearly
see the transition from $q_{-}$ to $q_{+}$ from $\Im(A)$, as the peaks of
$\Im(A)$ shift from being on the right of the peaks of $\Re(A)$ to the left, and
then back again. Figure~\ref{fig:EGL_example_1}$(b)$ shows~$|A|$, (half) the
amplitude of the reconstructed pattern $ u = A e^{i x} + \bar{A} e^{-i x}$ shown
in Figure~\ref{fig:EGL_example_1}$(c)$. Figure~\ref{fig:EGL_example_1}$(d)$ shows the
approximation to the local wavenumber of the reconstructed pattern.

The example in Figure~\ref{fig:EGL_example_2}, with~$\nu_1=0$, is similar, but
shows more pronounced beating between the two wavenumbers, evident in the
reconstructed solution and in Figure~\ref{fig:PDE_example}. In this case,
the localised patterns found in the model PDE are seen to be interpreted
correctly in terms of localised combinations of phase-winding solutions of the
complex Swift--Hohenberg equation, rather than as localised solutions of the
real Swift--Hohenberg equation, as suggested by \citet{Kozyreff2009}. The complex
Swift--Hohenberg equation does (since the coefficients are real) admit real
solutions, but these solutions appear to be unstable.

\subsection{Addition of the Proctor term}
\label{subsec:snaking}

The branches of localised solutions shown in Figure~\ref{fig:bifurcation_diag}
(without the Proctor term) do not close. We now investigate the addition of the
Proctor term into the model, considering~\eqref{eq:PDE_nonlinear_Proctor} with
a rather large value of $\mu_{p}=-0.65$, chosen as to make the effects of this
term more pronounced. We also fix $f_{1}=0.2814$ and $f_{2}=-0.0721$ so that
the marginal stability curve has two minima at different heights.

Seeking localised solutions, we plotted (as before) $\mathcal{H}$ and
$\mathcal{F}$ (modified to include the Proctor term), looking for zero contours
of $\mathcal{H}(k_{+})-\mathcal{H}(k_{-})$ and
$\mathcal{F}(k_{+})-\mathcal{F}(k_{-})$, but we found no intersection of
contours for the two minima at different heights for this choice of parameters.
Notwithstanding this, we returned to wavenumbers $k_{-}=58/64$ and
$k_{+}=70/64$ and used a localised solution constructed from these constituent
wavenumbers as a starting point for timestepping and continuation. Part of the
resultant localised solution branch is shown in
Figure~\ref{fig:bifurcation_diag_Proctor}. We see that including the Proctor
term allows localised solutions of different widths on the same branch, rather
than lying on distinct branches as in Figure~\ref{fig:bifurcation_diag}. The
existence of the localised solutions is also limited to a finite range of $\mu$
values; we expect the upper limit is introduced owing to the band of Eckhaus
stable wavenumbers closing.

Figure~\ref{fig:solutions_Proctor} shows solutions at the saddle-nodes
indicated in Figure~\ref{fig:bifurcation_diag_Proctor}. We notice that at each
of these saddle-nodes the proportion of each pattern in the domain varies.
Lower down the branch the pattern with the smaller wavenumber fills more of the
domain, and conversely higher up the branch the  pattern with the larger
wavenumber fills more of the domain. This behaviour is qualitatively similar to
the snaking behaviour of localised solutions in the subcritical
Swift--Hohenberg equation \citep{Burke2006}, in that moving up the snaking
branch adds to the width of the spatially periodic part of the localised
solution -- but the details, with saddle-node bifurcations appearing at many
different places along the branch, are considerably more complicated, typical
of snaking in more complicated situations such as hexagons \citep{Lloyd2008} or
quasipatterns \citep{Subramanian2018} in two dimensions.
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_10
\begin{tikzpicture
\begin{axis}[xtick={0.1,0.2},ytick={0,1},xticklabels={$-0.214$,$-0.114$},axis x line=bottom,axis y line=left,width=15cm,height=12cm,xmin=-0.0,xmax=0.275,ymin=0,ymax=1.82]
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p1.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p2.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc3.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc4.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p5.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc6.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p7.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc8.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc9.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc10.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc11.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc12.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc13.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc14.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p15.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc16.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc17.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p18.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p19.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc20.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_proc21.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p22.dat};
 \addplot[mark=none,color=blue] table[x index=0,y index=1]{ch6_p23.dat};
  \draw[-] (axis cs:0.192,1.46) -- (axis cs:0.202,1.46) -- (axis cs:0.202,1.62) -- (axis cs:0.192,1.62) -- cycle;
  \draw[-] (axis cs:0.175,1.455) -- (axis cs:0.179,1.455) -- (axis cs:0.179,1.475) -- (axis cs:0.175,1.475) -- cycle;
  \draw[-] (axis cs:0.08775,1.021) -- (axis cs:0.0879,1.021) -- (axis cs:0.0879,1.022) -- (axis cs:0.08775,1.022) -- cycle;
  \draw[-] (axis cs:0.0519,0.767) -- (axis cs:0.0521,0.767) -- (axis cs:0.0521,0.77) -- (axis cs:0.0519,0.77) -- cycle;
  \draw[-] (axis cs:0.048878,0.6868) -- (axis cs:0.048894,0.6868) -- (axis cs:0.048894,0.687) -- (axis cs:0.048878,0.687) -- cycle;

\end{axis}

\setlength{\fboxsep}{0pt}%
\setlength{\fboxrule}{0.5pt}%
\node [inner sep=0pt,outer sep=0pt] at (12,6) {\fbox{\includegraphics[trim=5mm 5mm 5mm
5mm, clip, width=3cm, height=4cm]{ch6_procsn1}}};
\node [inner sep=0pt,outer sep=0pt] at (8,4) {\fbox{\includegraphics[trim=5mm 5mm 5mm
5mm, clip, width=2.5cm, height=2.5cm]{ch6_procsn2}}};
\node [inner sep=0pt,outer sep=0pt] at (4.5,8.5) {\fbox{\includegraphics[trim=5mm 5mm 5mm
5mm, clip, width=2.5cm, height=2.5cm]{ch6_procsn3}}};
\node [inner sep=0pt,outer sep=0pt] at (1.5,6) {\fbox{\includegraphics[trim=5mm 5mm 5mm
5mm, clip, width=2.5cm, height=2.5cm]{ch6_procsn4}}};
\node [inner sep=0pt,outer sep=0pt] at (4,2) {\fbox{\includegraphics[trim=5mm 5mm 5mm
5mm, clip, width=2.5cm, height=2.5cm]{ch6_procsn5}}};

 \draw (0.5,5.5) node[below] {$(a)$};
 \draw (4,7.8) node[below] {$(b)$};
 \draw (12.9,7.6) node[right] {$(c)$};
 \draw (12.3,4.8) node[right] {$(d)$};
 \draw (9,5.1) node[below] {$(e)$};
 \draw (3,1.1) node[above] {$(f)$};

 \draw[-] (1.5,4.75) --(2.55,4.39);
 \draw[-] (4.5,7.25) --(4.27,5.85);
 \draw[-] (12,8) --    (9.8,8.35);
 \draw[-] (8,5.25) --  (8.6,8.33);
 \draw[-] (4,3.25) --  (2.4,3.93);

\draw (0,10.4) node[left] {$L_{2}$};
\draw (12.4,-0.1) node[below] {$\mu$};
\draw (0,-0.1) node[below] {$\mu_{c}$};
 \end{tikzpicture}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:bifurcation_diag_Proctor} Branch of localised solutions
 with $\mu_{p}=-0.65$, $f_{1}=0.2814$, $f_{2}=-0.0721$. The critical value
 of the bifurcation parameter $\mu$ is $\mu_{c}=-0.314$. The interior
 saddle-nodes are magnified and the labels $(a)-(e)$ correspond to
 the solution profiles shown in Figure~\ref{fig:solutions_Proctor}.}
 \end{figure}
 %
 \begin{figure}
 \begin{center}
 \mbox{\beginpgfgraphicnamed{paper_figure_11
\begin{tabular}{cc}
\begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.3,0.3},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msna.dat};
\end{axis}
\draw (0,1.5) node[left] {\hspace{1em} $u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(a)$};
\end{tikzpicture}
&
  \begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.3,0.3},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msnb.dat};
\end{axis}
\draw (0,1.5) node[left] {\hspace{1em} $u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(b)$};
\end{tikzpicture}
\\
\begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.5,0.5},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msnc.dat};
\end{axis}
\draw (0,1.5) node[left] {$u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(c)$};
  \end{tikzpicture}
&
  \begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.5,0.5},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msnd.dat};
\end{axis}
\draw (0,1.5) node[left] {$u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(d)$};
\end{tikzpicture}
\\
\begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.5,0.5},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msne.dat};
\end{axis}
\draw (0,1.5) node[left] {$u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(e)$};
  \end{tikzpicture}
&
\begin{tikzpicture}
\begin{axis}[xtick={0,0.995},xticklabels={0,400},ytick={-0.3,0.3},width=7cm,height=4.5cm,xmin=0,xmax=1]
\addplot[mark=none,color=blue] table[x index=0, y index=1]{ch6_p_msnf.dat};
\end{axis}
\draw (0,1.5) node[left] {$u$};
\draw (2.75,0) node[below] {$x$};
\draw (0.3,0) node[above] {$(f)$};
  \end{tikzpicture}
  \end{tabular}
 \endpgfgraphicnamed}%
 \end{center}
 \caption{\label{fig:solutions_Proctor} Solutions at the saddle-nodes
 indicated in Figure~\ref{fig:bifurcation_diag_Proctor}.}
 \end{figure}

\section{Discussion and conclusions}
\label{sec:conclusion}

The aim of this paper was to develop a new model
equation~\eqref{eq:PDE_nonlinear_Proctor} that captures qualitatively the
behaviour of pattern-forming problems with a quartic marginal stability curve,
and then explore the existence (and snaking) of localised solutions within this
model equation. There has been much progress in recent years developing a
framework for the understanding of localised solutions in the Swift--Hohenberg
equation \citep{Dawes2010,Knobloch2015}, where the existence of localised
patterns is interpreted in terms of a stable pattern existing at the same
parameter values as the stable trivial solution. Here we have shown that this
scenario holds in the case of the unfolding of a quartic minimum, where there
are coexisting stable patterns with two similar wavenumbers, and localised
solutions consisting of combinations of these. This work fits in with other
recent efforts that use Swift--Hohenberg-based models to explore problems with
localisation and snaking with multiple coexisting patterns
\citep{Knobloch2019,Alriheili2020,Subramanian2020}.

We computed the weakly nonlinear amplitude equation for the
model~\eqref{eq:PDE_nonlinear} and recovered a generalised Ginzburg--Landau
equation~\eqref{eq:generalised_GLE}. This allowed us to compute small-amplitude
solutions with different wavelengths and to classify the Eckhaus stable
patterns. Our work generalises that of \citet{Raitt1995} and \citet{Gelens2010}
(and corrects that of~\citet{Kozyreff2009}), to the case where the heights of
the minima in the marginal stability curve are different. We made use of a
first integral and a Lyapunov function to indentify candidate initial
conditions for finding localised solutions. Once the Proctor term, which allows
the wavenumber of maximum growth rate to depend on the bifurcation parameter,
is included in~\eqref{eq:PDE_nonlinear_Proctor}, the branches of localised
solutions join up, but the snaking we find is considerably more complicated
that the standard Swift--Hohenberg scenario.

An alternative approach to looking for localised solutions is to use spatial
dynamics, seeking only steady solutions of the model equation. For the
Swift--Hohenberg equation, this entails performing a normal form analysis of
the Hamiltonian--Hopf bifurcation, which describes the bifurcation from the
basic state, and contains the familiar homoclinic snaking of localised
solutions \citep{Woods1999}. In this framework, the existence of localised
solutions is determined by means of a geometric argument, whereby two integrals
of the normal form define a space which can be divided into regions which allow
or preclude the existence of localised solutions. In
Appendix~\ref{sec:appendix}, we derive a normal form for the bifurcation
occurring at the quartic minimum of the model equation, following the
derivation of the normal form for the Hamiltonian--Hopf bifurcation
\citep{Iooss1998a,Woods1999}. We use two different methods (following
\citet{Iooss1998a} and \citet{Burke2007a}) to find the coefficients in the
normal form. The two methods do not produce the same values for the
coefficients in the normal form, and in either case, there were terms appearing
in the normal form that one would have expected to be of higher order. As a
consequence of these terms we were unable to find normal form integrals,
including those one might expect by extension from the Hamiltonian--Hopf
analysis. This means we can not employ a similar geometric analysis to find
localised solutions. This normal form analysis merits further investigation. In
particular, it would be interesting to ascertain the reasons for the
difference between the normal form coefficients calculated by the two different
methods, and also to determine whether any variant of the normal form is
integrable; the variant we have derived does not seem to be.

Having found numerical localised solutions in the model equation, we could
ask next whether examples more closely connected to reality might also have
them, for example the magnetised Taylor--Couette system \citep{Stefani2009} and
rotating magnetoconvection \citep{Cox2001, Chandrasekhar1961}, where there are
suggestions that marginal stability curves can change from having one to having
two minima. This last example offers the interesting possibility of exploring
two-dimensional patterns, potentially with regions of small hexagons embedded
in a background of large hexagons, for example. This would be a natural
extension to the study of localised solutions in the two-dimensional
Swift--Hohenberg equation \citep{Lloyd2008,Subramanian2018}.

The coalescence of the two minima is in itself an interesting problem. In one
dimensional systems with marginal stability curves with two minima far apart,
the natural approach is to reduce the problem to two coupled second-order
Ginzburg--Landau equations \citep{Dawes2008a}, and there are model equations
based on the Swift--Hohenberg equation that allow pattern formation on two
length scales \citep{Muller1994,Lifshitz1997,Rucklidge2012}. It would thus be
interesting to investigate how one would transition from the situation of two
well separated length scales to the unfolding of a quartic minimum.

\section*{Acknowledgements}

Thomas Wagenknecht (1974--2012) co-supervised much of the work presented in
this paper: DCB and AMR gratefully acknowledge his perceptive insight and
infectious enthusiasm. This project was inspired by conversations with Rainer
Hollerbach. We also acknowledge useful discussions and correspondence with
Jonathan Dawes, G\'erard Iooss, Edgar Knobloch, Gregory Kozyreff, Eric
Lombardi, Michael Proctor and Priya Subramanian. This research was supported by
a PhD studentship from the Science and Technology Facilities Council (STFC).