# Collisions of localized patterns in a nonvariational Swift-Hohenberg
equation

Mathi Raja  Adrian van Kan  Benjamin Foster  Edgar Knobloch  Department of
Physics, University of California at Berkeley, Berkeley, California 94720, USA

###### Abstract

The cubic-quintic Swift-Hohenberg equation (SH35) has been proposed as an
order parameter description of several convective systems with reflection
symmetry in the layer midplane, including binary fluid convection. We use
numerical continuation, together with extensive direct numerical simulations
(DNSs), to study SH35 with an additional nonvariational quadratic term to
model the effects of breaking the midplane reflection symmetry. The
nonvariational structure of the model leads to the propagation of asymmetric
spatially localized structures (LSs). An asymptotic prediction for the drift
velocity of such structures, derived in the limit of weak symmetry breaking,
is validated numerically. Next, we present an extensive study of possible
collision scenarios between identical and nonidentical traveling structures,
varying a temperature-like control parameter. These collisions are inelastic
and result in stationary or traveling structures. Depending on system
parameters and the types of structures colliding, the final state may be a
simple bound state of the initial LSs, but it can also be longer or shorter
than the sum of the two initial states as a result of nonlinear interactions.
The Maxwell point of the variational system, where the free energy of the
global pattern state equals that of the trivial state, is shown to have no
bearing on which of these scenarios is realized. Instead, we argue that the
stability properties of bound states are key. While individual LSs lie on a
modified snakes-and-ladders structure in the nonvariational SH35, the multi-
pulse bound states resulting from collisions lie on isolas in parameter space,
disconnected from the trivial solution. In the gradient SH35, such isolas are
always of figure-eight shape, but in the present non-gradient case they are
generically more complex, although the figure-eight shape is preserved in a
small subset of cases. Some of these complex isolas are shown to terminate in
T-point bifurcations. A reduced model is proposed to describe the interactions
between the tails of the LSs. The model consists of two coupled ordinary
differential equations (ODEs) capturing the oscillatory structure of SH35 at
the linear level. It contains three parameters: two interaction amplitudes and
a phase, whose values are deduced from high-resolution DNSs using gradient
descent optimization. For collisions leading to the formation of simple bound
states, the reduced model reproduces the trajectories of LSs with high
quantitative accuracy. When nonlinear interactions lead to the creation or
deletion of wavelengths the model performs less well. Finally, we propose an
effective signature of a given interaction in terms of net attraction or
repulsion relative to free propagation. It is found that interactions can be
attractive or repulsive in the net, irrespective of whether the two closest
interacting extrema are of the same or opposite signs. Our findings highlight
the rich temporal dynamics described by this bistable nonvariational SH35, and
show that the interactions in this system can be quantitatively captured, to a
significant extent, by a highly reduced ODE model.

##  I Introduction

Spatially localized structures (LSs) are observed in a wide variety of
physical systems, from solitary water waves [1] to neurons [2], fluid
convection [3, 4], shear flows [5, 6] and reaction-diffusion systems [7, 8],
to name only a few. Generically, these systems are subject to dissipation and
require forcing to maintain the structure, see [9] for a review of spatial
localization in such systems. A simple model of pattern formation in forced
dissipative systems is provided by the bistable Swift-Hohenberg equation,
originally suggested in the context of pattern formation in Rayleigh-Bénard
convection [10, 11]. This equation supports well-known localized solutions
that are organized in a snakes-and-ladders bifurcation structure [12, 13, 14].

When the Swift-Hohenberg equation has gradient structure, solutions with
nontrivial time dependence are precluded. However, nongradient generalizations
of the Swift-Hohenberg equation arise frequently in applications [15] and
these permit both time dependence, see e.g. [16], and persistent propagation,
e.g. [17]. In this paper, we consider a specific instance of such models,
namely the one-dimensional cubic-quintic Swift-Hohenberg equation with broken
reflection symmetry,

$\partial_{t}u=ru-(1+\partial_{x}^{2})^{2}u+b_{3}u^{3}-u^{5}+\epsilon(\partial_{x}u)^{2}.$
(1)

Here the parameter $\epsilon$ controls both the nongradient structure of the
equation (the equation has gradient structure when $\epsilon=0$) and the
breaking of the symmetry $u\rightarrow-u$ (the equation is invariant under
$u\rightarrow-u$ when $\epsilon=0$). Since the equation is also symmetric
under spatial reflections, $x\rightarrow-x$, both effects are required for
spontaneous propagation of LSs in this system.

Equation (1) was suggested in [18] as a model of binary fluid convection with
broken midplane reflection symmetry and its properties are indeed in
qualitative agreement with direct numerical simulations of the Navier-Stokes
equations describing this system [19]. The present work extends significantly
the collision studies undertaken in [18] and clarifies a number of key issues.

Figure 1: Bifurcation diagram for Eq. (1) with $\epsilon=0$. The patterned
state branch is shown in blue with a sample solution profile in red (lower
left inset). The snaking branches of symmetric ($L_{0},L_{\pi}$) and
antisymmetric ($L_{\pi/2},L_{3\pi/2}$) LSs are shown in black (sample profiles
in yellow and green, respectively). For clarity, only three of the
interconnecting rung states are shown (detail in upper right inset), cf. [18].
Larger norm indicates longer LSs. Stable solutions are found on branches with
positive slope.

In the following, we refer to Eq. (1) as SH35. When $\epsilon=0$, the problem
reduces to variational form such that, on a domain of any size $L$, there
exists a free energy functional

$\mathcal{F}[u(x)]=\int_{0}^{L}\left(-\frac{1}{2}ru^{2}+\frac{1}{2}[(1+\partial_{x}^{2})u]^{2}-\frac{b_{3}}{4}u^{4}+\frac{u^{6}}{6}\right)dx,$
(2)

with the property that $\partial_{t}u=-\frac{\delta\mathcal{F}}{\delta u}$.
Thus $\mathcal{F}$ decreases along trajectories towards local minima as
$t\rightarrow\infty$, and these represent stable steady states of the system.
The free energy of spatially periodic patterns passes through zero at a
$r=r_{M}$, known as the Maxwell point. In the vicinity of the Maxwell point
one can create a variety of localized structures involving both the pattern
state and the trivial state $u=0$ at little or no cost in energy. These are of
two types, localized even solutions, denoted here as $L_{0}$ ($L_{\pi}$) if
their maximum (minimum) is located at their center, and localized odd
solutions, denoted as $L_{\frac{\pi}{2}}$ ($L_{\frac{3\pi}{2}}$) if they have
a negative (positive) slope at the center.

When $0<\epsilon\ll 1$, the symmetry $u\to-u$ as well as the variational
structure is broken but similar solutions continue to exist, albeit with
modified properties, as described in [18]. The symmetric solution branches
$L_{0},L_{\pi}$ remain symmetric and stationary. Odd solution profiles cease
to be odd and hence propagate. In this paper, we first study the propagation
of isolated structures and then go on to investigate in detail the collisions
that can result. In contrast to the collisions familiar from studies of
integrable partial differential equations on the real line, here the
collisions are inelastic and can lead to annihilation and sticking as well as
scattering. One question of interest concerns the role, if any, played by the
Maxwell point of the variational system in such collisions when $\epsilon$ is
small: for example, is the collision process accompanied by nucleation or
annihilation of new wavelengths according to the free energy minimization
principle valid at $\epsilon=0$ or, if this principle is not followed, what
other mechanism determines collision outcomes?

The remainder of this paper is structured as follows. In section II, we
describe the general bifurcation structure of Eq. (1) obtained from numerical
continuation. Next, in section III we present an asymptotic computation of the
drift speed of asymmetric LS in the limit of weak symmetry breaking
$0<\epsilon\ll 1$, whose accuracy is confirmed by comparison with direct
numerical simulations (DNSs) and numerical continuation. In section IV, we
present the results of extensive DNSs of all possible collision scenarios in
this system, and describe the dependence of the collision outcome on the
control parameter $r$, with the stability of multi-pulse bound states playing
a key role. We also show that the bound-state solutions arising from
collisions are in many cases a part of non-trivial isolas, whose structure
depends on the symmetry breaking parameter $\epsilon$. In section V, we
present a reduced model of the interactions between colliding patterns, which
is based on the linear structure of SH35 and accurately captures the
trajectories of LSs so long as no significant nonlinear interactions creating
or destroying wavelengths occur. The paper concludes in section VI with a
discussion of our results. All our results are obtained for the choice
$b_{3}=2$ as used in [18], employing periodic boundary conditions on a domain
of size $L=40\pi$.

Figure 2: Free energy $F[u]$, defined in Eq. (2), of the periodic pattern
state versus $r$. Thick (thin) lines correspond to stable (unstable)
solutions. The dotted vertical red line indicates the Maxwell point
$r_{M}\approx-0.675$, where the free energy changes sign. Figure 3:
Bifurcation diagram for Eq. (1) with $\epsilon=0.03$. Stable Z branch states
with at least two wavelengths exist in the range $-0.70\lesssim
r\lesssim-0.63$. Inset: sample solution profiles at color-coded locations in
ascending order, cf. [18]. Larger norm indicates longer LSs. Stable solutions
are found on branches with positive slope.

##  II Bifurcation analysis

In this section we provide an overview of the solution structure of Eq. (1),
first in the case $\epsilon=0$, and then in the case $\epsilon>0$, obtained
from numerical continuation using AUTO [20] and pde2path [21]. The results are
presented in terms of the $L_{2}$ norm of the solutions $u(x)$ given by

$\|u\|_{2}\equiv\sqrt{\frac{1}{L}\int_{0}^{L}u^{2}(x)dx}$ (3)

and provide important background information for subsequent sections.

###  II.1 The variational case: $\epsilon=0$

We first consider the variational problem with $\epsilon=0$ and employ
numerical continuation to construct the bifurcation diagram shown in Fig. 1,
cf. [13]. A trivial flat state $u=0$ exists at all $r$, and undergoes a
subcritical bifurcation to a periodic pattern at $r=0$ (red profile in lower
left inset in Fig. 1). Four snaking branches bifurcate from the periodic
pattern branch in secondary bifurcations, corresponding to the symmetric LSs
$L_{0}$, $L_{\pi}$ and the antisymmetric LSs $L_{\pi/2}$, $L_{3\pi/2}$
mentioned in the introduction. The branches overlap in pairs owing to the
symmetry $u\to-u$, and both sets display characteristic snaking behavior in
the vicinity of the Maxwell point $r_{M}$ [13, 9]. In addition, rung states
connect the symmetric and antisymmetric snaking branches, arising in pitchfork
bifurcations close to every saddle-node bifurcation on these branches [18].
Each rung actually corresponds to four branches of unstable asymmetric
localized solutions, related by the symmetries $x\to-x$ and $u\to-u$ [14].
Examples of symmetric and antisymmetric solution profiles are shown in the
lower left inset in Fig. 1, together with the periodic pattern state. When the
LSs fill the available domain the snaking branches reconnect to the periodic
pattern.

Figure 2 shows the free energy of the periodic pattern from Fig. 1 versus $r$.
The stable part of the periodic pattern state (thick blue line) has a free
energy which decreases monotonically with $r$ and changes sign at the Maxwell
point $r_{M}\approx-0.675$.

###  II.2 The nonvariational case: $\epsilon>0$

Here we describe how the bifurcation structure changes in the nonvariational
case, specifically for $\epsilon=0.03$ (Fig. 3). As in the case $\epsilon=0$,
there is a trivial branch $u=0$, a periodic pattern branch emerging
subcritically from it, and two snaking branches. However, the two snaking
branches in Fig. 3 now correspond to $L_{0}$ and $L_{\pi}$, since these states
are no longer related by symmetry, and consequently snake in phase before
reconnecting to the periodic state. Moreover, the $L_{\pi/2}$, $L_{3\pi/2}$
states and the rung states reconnect, forming a sequence of ’Z’-shaped
branches consisting of asymmetric solutions.

As a consequence of the nonvariational structure of the problem when
$\epsilon>0$, these asymmetric solutions drift and hence may collide. The
structure of several Z branch solutions at different $r$ values is shown in
the inset in Fig. 3. It is important to observe that the Z branch states may
be stable or unstable. Stable solutions are present on the “diagonal” part of
each Z, i.e. in the range $-0.70\lesssim r\lesssim-0.63$ (except for the
lowest branch, which is stable only between $-0.65\lesssim r\lesssim-0.62$).
The corresponding drift velocity $c$ depends on $r$, as shown in Fig. 4 for
several of the states depicted in Fig. 3. In Fig. 4, the stable part of any
given Z branch corresponds to the segment between $-0.70\lesssim
r\lesssim-0.63$, where the drift speed is largest. We emphasize that longer
asymmetric LSs drift more slowly than short ones.

Figure 4: Drift velocity $c$ versus $r$ for various Z branch states at
$\epsilon=0.03$; colors correspond to the branches in Fig. 3. The curves shown
correspond to the Z branches for which profiles are displayed in Fig. 3.

##  III Drift velocity of isolated structures

###  III.1 Asymptotic theory

To compute the drift velocity of asymmetric LSs using perturbation theory in
the limit of weak symmetry breaking, $0<\epsilon\ll 1$, we introduce the fast
and slow time variables

$\tau=\epsilon^{\frac{1}{2}}t,\qquad T=\epsilon t$ (4)

and denote their spatial phase by $\theta(T)$. Following the calculation
presented in [17], we posit the expansion

$\displaystyle u(x,t)=\hphantom{l}$ $\displaystyle
U_{0}[x-\theta(T)]+\epsilon^{\frac{1}{2}}u_{1}[x-\theta(T),\tau]$
$\displaystyle+\epsilon u_{2}[x-\theta(T),\tau]+o(\epsilon),$ (5)

where $U_{0}$ is a known moving pattern whose propagation velocity we seek to
determine. At leading order, Eq. (1) implies

$rU_{0}-(1+\partial_{x}^{2})^{2}U_{0}+b_{2}U_{0}^{3}-U_{0}^{5}=0,$ (6)

while at $O(\epsilon^{\frac{1}{2}})$ we obtain

$\mathcal{L}u_{1}=ru_{1}-(1+\partial_{x}^{2})^{2}u_{1}+3b_{2}U_{0}^{2}u_{1}-5U_{0}^{4}u_{1}=0,$
(7)

where the operator $\mathcal{L}\equiv
r-(1+\partial_{x}^{2})^{2}+3b_{2}U_{0}^{2}-5U_{0}^{4}$ is self-adjoint.

Differentiating Eq. (6) with respect to $x$, one finds that
$\mathcal{L}U_{0}^{\prime}=0$, i.e. if $U_{0}$ solves (6), then
$U_{0}^{\prime}$ solves (7). In [17], the authors considered $r$ to be
asymptotically close to the edge of the snaking interval, with a focus on the
dynamics of depinning. This required taking into account two additional
solutions which lie in the null space of $\mathcal{L}$: one symmetric and one
asymmetric mode. Here, we instead consider asymmetric states which are located
within the snakes-and-ladders structure, away from the saddle-nodes of the
snaking branches. Consequently $U_{0}^{\prime}$ is the only function in the
null space of $\mathcal{L}^{\dagger}=\mathcal{L}$. Since this translation mode
is already included in the Ansatz (5) we can set $u_{1}\equiv 0$. Finally, at
$O(\epsilon)$, we obtain

$-U^{\prime}_{0}\theta_{T}=\mathcal{L}u_{2}+\epsilon(U^{\prime}_{0})^{2}.$ (8)

Multiplying by $U_{0}^{\prime}$ and integrating over the domain, using the
fact that $\mathcal{L}$ is self-adjoint and that
$\mathcal{L}U_{0}^{\prime}=0$, we obtain

$-\theta_{T}\int(U^{\prime}_{0})^{2}dx=\int
U^{\prime}_{0}\mathcal{L}u_{2}+\epsilon(U^{\prime}_{0})^{3}dx=\int\epsilon(U^{\prime}_{0})^{3}dx.$
(9)

Rearranging, we arrive at the drift velocity,

$\displaystyle
c\equiv\theta_{T}=-\epsilon\frac{\int_{0}^{L}\left(U_{0}^{\prime}(x)\right)^{3}dx}{\int_{0}^{L}\left(U_{0}^{\prime}(x)\right)^{2}dx},$
(10)

a special case of a more general result [22].

Equation (10) is invariant under $(\epsilon,U_{0})\to(-\epsilon,-U_{0})$, a
symmetry that is inherited from Eq. (1). For symmetric profiles $U_{0}(x)$,
the numerator vanishes; symmetric states therefore remain at rest even when
$\epsilon>0$. However, for the asymmetric profiles shown in Fig. 3 the
velocity $c$ is nonzero. This is due to the nonsinusoidal nature of these
profiles and in particular the contribution from the fronts at either end of
the LS profile. We define the number $n$ of significant extrema in any given
LS as the number of extrema whose amplitude is larger than $1/e$ times the
maximum amplitude in the LS, and the number of significant wavelengths as
$1/2$ the number of significant extrema. The adjective significant will be
dropped in the following when there is no ambiguity. For solutions with many
significant wavelengths, $n\gg 1$, the numerator converges to a constant,
while the denominator grows approximately linearly with $n$. This implies that
$c\propto 1/n$ at large $n$, which is quantitatively consistent with the
results from numerical continuation, as shown in Fig. 5: longer structures
move more slowly, as already highlighted in the discussion of Fig. 4.

Figure 5: Log plot showing $c(r=-0.67)$ (blue) vs. number of wavelengths $n$
overlaid with $1/n$ (black). A satisfactory agreement is observed. Figure 6:
Comparison of the theoretically predicted drift velocity $c$ with the
corresponding DNS results for stable states on the second lowest Z branch.
Blue solid line indicates the theoretical prediction from Eq. (10), with the
integrals evaluated numerically. The green dash-dotted line shows the velocity
obtained from numerical continuation. Red circles indicate the DNS results.
The DNS and numerical continuation are in excellent agreement, including the
deviation from asymptotic behavior at larger $\epsilon$. The asymptotic theory
is in close quantitative agreement with both DNS and numerical continuation at
small $\epsilon$. The Z branch states cease to exist above $\epsilon=0.13$.
Inset shows the Z branch solution profiles at $\epsilon=0$ (blue,
antisymmetric) and $\epsilon=0.13$ (red, asymmetric) at $r=-0.66$.

###  III.2 Numerical verification

To determine the range of validity of the prediction in Eq. (10), we perform
direct numerical simulations (DNSs) of Eq. (1) for selected values of the
parameter $\epsilon$. We also perform numerical continuation in $\epsilon$.
The results reported below extend those in [22] where the asymptotic result
was compared with numerical continuation at only one value of $\epsilon$,
$|\epsilon|=0.01$. For this purpose, we consider a one-dimensional periodic
domain of the same length $L=40\pi$ as in the numerical continuation and solve
Eq. (1) using a semi-implicit pseudo-spectral numerical scheme for spatial
derivatives [23] and a fourth-order Runge-Kutta time-stepping scheme. All DNS
results presented in this paper were obtained on a finely resolved uniform
spatial grid of $8192$ grid points (corresponding to approximately $100$ grid
points per pattern wavelength), except when specified otherwise. We use a time
step of $dt=0.01$. Larger values of $dt$ led to incorrect drift speeds.

To verify Eq. (10), we consider a state on the second lowest Z branch at
$\epsilon=0.13$, $r=-0.66$, shown in the inset in Fig. 6. The state is
asymmetric: a small change in peak/trough amplitudes with increasing
$\epsilon$ can be discerned, a consequence of symmetry breaking when
$\epsilon>0$.

Figure 6 shows the drift velocity $c$ as a function of $\epsilon$, with
excellent quantitative agreement between DNSs, numerical continuation and
theory for $\epsilon\lesssim 0.1$. For $\epsilon\gtrsim 0.1$, the numerically
obtained drift velocities deviate from the asymptotics, with the asymptotic
prediction overestimating the numerical value by less than $10\%$. However,
the DNSs and numerical continuation continue to show excellent agreement. At
$\epsilon\approx 0.13$ the state under consideration ceases to exist, since
beyond this point the bifurcation structure of the system is qualitatively
altered (cf. [18]). Additional runs with the smaller time step $dt=0.001$ were
also performed, and gave the same results as those for $dt=0.01$, indicating
that the DNSs are well resolved in time. We also repeated the DNSs on a
coarser grid with $2048$ grid points and obtained drift velocities $c$ that
are indistinguishable from those shown in Fig. 6, indicating that the
simulations are also well resolved in space.

Figure 7: Space-time plots of collision scenarios A-D for five different
values of $r$ indicated by gray dashed vertical lines in Fig. 12, shown with
time along the vertical axis and space along the horizontal axis. The
different collision scenarios involve different time scales, therefore
distinct time axes are specified for each case. A rich zoology of different
collision outcomes is found: one or two extrema may be deleted (for $r$ near
$-0.7)$, or bound states can form without the creation or deletion of any
extrema (in all scenarios except scenario D), or one, three, four or five
extrema may be added in the collision, depending on the value of $r$ and on
the scenario. Some of the resulting states travel while others are symmetric
and hence stationary.  Figure 8: The $L_{2}$ norm versus time for a collision
between two chasing asymmetric structures from scenario A ($r=-0.635$)
consisting of two stages. In the first stage, the smaller asymmetric LS is
rendered symmetric and hence stationary by the addition of an extremum. In the
second stage, the larger asymmetric LS suffers the same fate resulting in a
stationary bound state. Figure 9: Annihilation of a single-wavelength
asymmetric structure in collision with a symmetric LS at $r=-0.65$. This is a
special case of scenario B with an asymmetric structure consisting of only two
extrema, located at the leftmost edge of the range of stability of the lowest
Z branch in Fig. 3. Note the coarser colorbar compared to Fig. 7.

##  IV Collisions of localized structures

As shown in Fig. 3, multiple stable LSs can coexist in the domain, and these
may undergo collisions. In this section we present the results of extensive
DNSs of such collisions between different types of LSs performed using the
numerical solver described in the previous section. For all results described
below, we take $\epsilon=0.03$, unless stated otherwise.

###  IV.1 Overview of DNS results

Here we describe the rich collision phenomenology that is observed when
different types of LSs collide at various values of $r$, as illustrated in
Fig. 7. We distinguish the following four collision scenarios:

  * •

Scenario A: collision between two asymmetric states which differ in the number
of significant wavelengths, and travel in the same direction, with the
shorter, faster LS catching up to the longer, slower LS, as predicted by Eq.
(10). The two colliding extrema are of opposite sign. We focus specifically on
collisions between two LSs of length two and three wavelengths. Four different
possible outcomes are observed: deletion of one extremum, formation of a
drifting bound state without a change in the number of extrema, and the
creation of one or four new extrema. Note that the bound states at $r=-0.68$
and $r=-0.645$ shown in Fig. 7 differ in their separation. The cases $r=-0.7$
and $r=-0.635$ in fact involve two separate consecutive collisions due to the
periodicity of the domain, of which only the first is shown in Fig. 7. In the
former case, the second collision (not shown) results in a drifting bound
state, while in the latter case (shown in Fig. 8), the asymmetric LS is
rendered symmetric and stationary owing to the nucleation of an additional
extremum, resulting in a larger but stationary bound state.

  * •

Scenario B: a drifting asymmetric state collides with a stationary symmetric
state with a maximum at its center. The two colliding extrema are of opposite
sign. We focus on a two-wavelength asymmetric LS and a symmetric LS with three
positive and two negative large-amplitude extrema. Four different possible
collision outcomes are observed: deletion of one extremum, formation of a
drifting bound state without change in the number of extrema, and the creation
of one or four new extrema.

  * •

Scenario C: same as scenario B but with the symmetric LS flipped by $u\to-u$
so that the two colliding extrema are of the same sign. We focus on a two-
wavelength asymmetric LS and a symmetric structure with three negative and two
positive large-amplitude extrema. Five different possible collision outcomes
are observed: deletion of one extremum, formation of a drifting bound state
without change in the number of extrema, and the creation of one, three or
five new extrema.

  * •

Scenario D: head-on collision between two asymmetric states. The two colliding
extrema are of the same sign. We focus on a pair of identical two-wavelength
patterns (collisions between asymmetric LSs of distinct sizes were also
considered, and gave qualitatively similar results). Four different possible
collision outcomes were observed: deletion of two extrema and the creation of
three or five new extrema. No pure bound states were observed in this case.

Figure 7 shows a summary of the space-time plots depicting the different
possible collision outcomes. All collisions were simulated using DNS on a
uniform grid of $2048$ points to facilitate longer simulation times. Some
cases were also repeated with $8192$ grid points and no change in the
collision behavior was observed.

Three qualitatively distinct collision outcomes are observed across all
collision scenarios: deletion of extrema, formation of a bound state, and the
creation of new extrema. Some of these states travel while others are
symmetric and hence stationary. One notices that scenarios A and B are quite
similar to one another in terms of the outcome realized at a given $r$ (except
near $r=-0.645$, see also Fig. 12). Scenarios C and D are similar in the same
sense. Since the closest extrema in the collisions in scenarios A and B are
both of opposite sign but are of the same (negative) sign in scenarios C and
D, one might surmise that the relative sign of interacting extrema is
important in determining the qualitative collision outcome, cf. [18]. Our
detailed results broadly support this suggestion but indicate that this sign
is not the sole factor determining the collision outcome since differences
between scenarios persist.

Figure 7 reveals that the speed $c$ of the traveling state plays a significant
role in determining the collision outcome. This speed is controlled by the
value of the parameter $r$ since $r$ controls the degree of asymmetry of each
traveling state, but it also depends on the length of the structure, cf. Fig.
4. For example, in scenario A at $r=-0.7$ a narrow LS catches up to a wider
LS, with the resulting interaction stopping the former and leaving the latter
unaffected. Because of periodic boundary conditions the wider, slower moving
LS collides with the stationary state in a subsequent collision, creating a
drifting bound state (not shown, but see Fig. 8 for a similar multiple
collision at $r=-0.635$). Such drifting bound states are not observed in any
of the other scenarios at this value of $r$ where a narrower and so faster LS
interacts with a broader LS at rest. Figure 9 shows an extreme case of
scenario B, where a single-wavelength asymmetric LS collides with a stationary
symmetric structure and is annihilated. This is because this asymmetric LS is
located on the lowest Z branch in Fig. 3 and so is minimal in the sense that
no stable LS shorter than one significant wavelength exists. Hence, the
deletion of an extremum implies annihilation of this LS.

###  IV.2 Change in $L_{2}$ norm before and after collision

To quantify the change in pattern size during a collision, we measure the
$L_{2}$ norm defined in Eq. (3). Figure 10 shows a sample time series of
$\|u\|_{2}$ from a chasing collision (scenario A) at $r=-0.65$, where the
collision first produces a metastable bound state, but three new extrema are
generated at late times by nonlinear interactions. Figure 10 suggests the
simple metric

$\Delta\|u\|_{2}=\|u_{\textrm{final}}\|_{2}-\|u_{\textrm{initial}}\|_{2},$
(11)

where $u_{\textrm{initial}}$ is the state before the collision, i.e., when the
distance between the two structures is still large, and $u_{\textrm{final}}$
is the resulting state a long time after the collision has occurred.
Importantly, the collision dynamics are independent of the initial distance
between the colliding LSs, since there is no inertia in the system. This is
illustrated in Fig. 11, where the initial distance in scenario B (symmetric-
asymmetric collision) at $r=-0.65$ is varied by fractions of the wavelength
associated with the spatial eigenvalue, $\lambda=2\pi/\beta$, with $\beta$
defined in Eq. (12). Figure 11 shows that the collision dynamics are self-
similar under shifts in the initial distance: simply accounting for the time
required to propagate over the additional distance leads to data collapse.
This has also been verified explicitly for other collisions, including from
scenario D. We deduce from these observations that the only variables
affecting collision outcomes are indeed the chosen pair of colliding
structures, and the control parameter $r$ (at fixed $\epsilon$).

The phenomenology illustrated in Fig. 10 is reminiscent of collision dynamics
described in terms of scattors [24, 25]: unstable stationary or time-periodic
patterns which direct the evolution in state space during the collision
process along their stable and unstable manifolds. While these ideas were
proposed in the context of models other than SH35, they may be applicable for
some of the collision events observed here.

Figure 10: The $L_{2}$ norm versus time for a collision of two chasing
asymmetric LSs (scenario A) at $r=-0.65$. The collision of the two chasing
structures occurs at $t\approx 20~{}000$, leading to the creation of a
traveling metastable bound state, whose structure changes slowly over time, as
visible from the slope in $\|u\|_{2}$. At $t\approx 40~{}000$, three
additional extrema are nucleated by nonlinear interactions, behavior resulting
in the jump in $\|u\|_{2}$. Figure 11: The $L_{2}$ norm versus time for
symmetric-asymmetric collisions from scenario B (at $r=-0.65)$ for initial
conditions shifted by different fractions $\Delta x$ of the wavelength
$\lambda=2\pi/\beta$, with $\beta$ the imaginary part of the spatial
eigenvalue defined in Eq. (12). Inset shows the same data with time shifted by
$\Delta t=\Delta x/c$, where $c$ is the drift speed of the asymmetric LS: all
curves collapse exactly, indicating that the collision dynamics are
independent of initial separation.

Figure 12 shows $\Delta\|u\|_{2}$ as a function of $r$ for all four scenarios.
A value of $\Delta\|u\|_{2}\approx 0$ corresponds to the formation of a bound
state without a change in the overall pattern size, a positive value indicates
the creation of one or more extrema, while a negative value indicates that one
or more extrema were deleted.

The range of $r$ shown in Fig. 12 corresponds to the interval with stable
propagating solutions when $\epsilon=0.03$. Vertical gray dashed lines
indicate the cases depicted in Fig. 7. Near the leftmost edge of the existence
interval, at $r=-0.7$, one (A-C) or two extrema (D) are deleted in the
collision. Near the rightmost edge, at $r=-0.635$, new extrema are created:
either one (A, B) or five extrema (C, D). Away from these boundary cases, the
scenarios differ more substantially. At $r=-0.68$, one observes either the
formation of bound states (A, B), or the creation of three additional extrema
(C, D). At $r=-0.66$, new extrema are added in the collision in all cases:
either four additional extrema (A, B) or three (C, D). At $r=-0.645$ and
$r=-0.65$ (not shown), scenarios A and C revert to $\Delta\|u\|_{2}\approx 0$,
i.e. the formation of a bound state, while new extrema continue to be added in
scenarios B (four extrema) and D (three extrema).

As already mentioned in the discussion of Fig. 7, scenarios A and B are
similar to one another, as are C and D, in the sense that for most values of
$r$, they display the same number of extrema added or deleted. However, Fig.
12 reveals deviations between these scenarios in the interval $-0.65\lesssim
r\lesssim-0.64$. To ensure that these are not numerical artefacts, we repeated
the runs in scenarios A and B in this range of $r$ at higher spatial
resolution, using $8192$ instead of $2048$ grid points, and continued the
simulation until very late times ($t\approx 150~{}000$). We also repeated the
runs with significantly smaller time steps, $dt=0.005$ and $dt=0.001$, with
the same collision outcome, confirming that the nonmonotonic dependence of
$\Delta\|u\|_{2}$ on $r$ is a robust result that we discuss further below.

The black dash-dotted vertical line in Fig. 12 indicates the location of the
Maxwell point $r=r_{M}$ for the variational case $\epsilon=0$. The figure
shows that the values of $r$ where $\Delta\|u\|_{2}$ changes differ between
scenarios and that, in addition, these locations lie far from $r=r_{M}$. We
conclude that the Maxwell point of the variational case has little, if any,
relevance in determining the collision outcome in the nonvariational case,
even in the weak symmetry-breaking case considered here.

The overall trend of increasing $\Delta\|u\|_{2}$ with increasing $r$ can be
viewed as a reflection, in the nonvariational regime, of similar behavior in
the variational case: when $\epsilon=0$, the free energy $\mathcal{F}$ of the
stable pattern state decreases with increasing $r$, as shown in Fig. 2, and
the periodic pattern becomes increasingly energetically favored. Apart from
some exceptions at $r\geq-0.65$ (to be discussed further below), the
increasing trend in the number of extrema added in a collision is compatible
with this intuition. Finally, even for a fixed number of extrema added or
lost, Fig. 12 reveals a slow decrease in $\Delta\|u\|_{2}$ with increasing
$r$. This is associated with small changes in the profile of the extrema as
$r$ changes.

Figure 12: Change in the $L_{2}$ norm before and after a collision,
$\Delta\|u\|_{2}$ defined in Eq. (11), for different collision scenarios,
shown as a function of $r$. Labels indicate the number of extrema gained or
lost in the collision. Vertical dashed gray lines correspond to the values of
$r$ shown in Fig. 7. The vertical black dash-dotted line shows the location of
the Maxwell point $r_{M}$ for $\epsilon=0$ (all remaining results shown are
for $\epsilon=0.03$).

###  IV.3 Bound states: isolas and stability

Figure-eight isola structures extending over the entire width of the snaking
region have been previously described in the quadratic-cubic Swift-Hohenberg
equation without nonvariational terms, arising from multi-pulse solutions
consisting of two LSs bound together in close proximity [26, 27]. Here, we
find that for bound states formed from collisions at $\epsilon>0$, the
properties of the isolas to which these states belong depend strongly on the
collision scenario.

Figure 13: Traveling bound states obtained from numerical continuation of
collision results in scenarios A and B at $\epsilon=0.03$. Thin pink line:
traveling state obtained from a chasing collision (scenario A) at $r=-0.68$.
Thin light blue line: similar traveling state obtained from a symmetric-
asymmetric collision (scenario B) at $r=-0.69$. Starting points for
continuation in $r$ are marked with circles. Continuation in $\epsilon$ to the
variational case $\epsilon=0$ yields isolas of stationary two-pulse states.
Thick red line: scenario A at $\epsilon=0$. Thick dark blue line: scenario B
at $\epsilon=0$. Figure 14: Isola structures (colored) corresponding to the
traveling two-pulse bound states generated in different collision scenarios
superposed on the bifurcation diagram from Fig. 3 (gray). All results are for
$\epsilon=0.03$.

Figure 13 shows two isolas at $\epsilon=0.03$ (thin lines) obtained from a
continuation in $r$, starting from two bound states generated by a collision
in scenarios A and B at $r=-0.68$ and $r=-0.69$, respectively (points
highlighted in Fig. 13), together with the results of continuing these isolas
in $\epsilon$ to the variational case $\epsilon=0$ (thick lines). The figure
shows that in the former case (red curves) the isola retains its shape as it
transforms into the corresponding isola of the variational problem. In the
latter and more typical case, the isola at $\epsilon=0.03$ takes the form of a
spaghetti-like tangle (light blue), requiring substantial simplification prior
to reaching the corresponding isola of the variational problem (thick blue).
Figure 14 shows several additional isolas of traveling bound states at
$\epsilon=0.03$ obtained by colliding longer initial LSs and superposed on the
full bifurcation diagram (colors match between Figs. 13 and 14). All figure-
eight isolas seen in Fig. 14 are from scenario A, while tangled isolas are
from scenario B.

It is important to note that the rightmost edge of the figure-eight isola in
scenario A shown in Fig. 13 (thin red line) is located near $r\approx-0.66$,
and not at $r\approx-0.63$, the right boundary of the snaking region for
$\epsilon=0.03$. This is a consequence of the fact that two-pulse states with
a small separation between the LSs lie on isolas that are significantly
smaller than those for bound states with a larger separation. This is so not
only when $\epsilon=0$ [26] but also when $\epsilon=0.03$. Figure 15 shows an
example of a narrow figure-eight isola at $\epsilon=0$ obtained from a bound
state at $\epsilon=0.03$ and $r=-0.68$ together with a wider figure-eight
isola obtained by continuing a bound state, also obtained from a collision in
scenario A with $\epsilon=0.03$, but at $r=-0.645$. Despite their similarity
(see profiles in Fig. 15, top panel) these profiles are part of separate
isolas: the solution profiles of the narrow isola feature an interaction
region between the two LSs comprising the bound state which is one wavelength
shorter than that of the larger isola solution profiles, a property that is
preserved upon numerical continuation. As the separation between the two LSs
decreases and their interaction becomes stronger, the figure-eight isolas
shrink, with the rightmost edge moving farther and farther to the left.

This decrease in size plays a key role in determining collision outcomes: the
value $r\approx-0.66$ (Fig. 13, thin read line) coincides with the parameter
value in Fig. 12 where new extrema are first created in scenario A. This
observation suggests the following way of understanding the observed
phenomenology: when two LSs approach one another in a collision, there are two
possibilities. Either (a) a stable multipulse bound state exists at the value
of $r$ in question for the given pair of LSs, or (b) no such state exists or
at least it is not stable. In case (a), the collision will lead to the
formation of a stable bound state. In case (b), on the other hand, no such
state is available, and the system creates or deletes extrema to reach another
stable solution.

With this hypothesis, the nonmonotonic behavior in the number of extrema
created summarized in Fig. 12 can be explained as follows: at parameter values
$-0.65\lesssim r\lesssim-0.64$, there is an island of stability where stable
bound states exist. When no stable bound state exists, a strongly nonlinear
interaction must occur, resulting in LSs with a different number of extrema;
this number increases with increasing $r$, as discussed earlier.

To test the hypothesis that the existence and stability of a bound state is
the determining factor in setting the collision outcome, we perform stability
experiments. Specifically, to verify that the bound state either does not
exist or is unstable in scenarios A, B with $-0.66\lesssim r\lesssim 0.65$
(where new extrema form), we performed DNS at these values of $r$,
initializing with the post-collision final states from $r=-0.68$, $r=-0.67$,
$r=-0.65$ and $r=-0.645$. In each of these cases we observed the eventual
generation of additional extrema. Only one exception was observed, at
$r=-0.66$ in scenario B, when a bound state from $r=-0.68$ was used as the
initial condition and found to remain a pure bound state indefinitely, without
any change in the number of extrema. Thus collisions may trigger the creation
of new extrema despite the existence of a stable bound state, provided the
bound state has only a small basin of attraction.

We repeated the above study for scenario C: stable bound states formed from
collisions at $r=-0.645$ and $r=-0.65$ were used as initial conditions at
$r=-0.635$, $r=-0.64$, $r=-0.655$, $r=-0.66$. In all cases but one, the same
number of extrema was spontaneously created as observed in the collision
experiments. The exception, where a bound state remains indefinitely stable,
was again found near the transition from a bound state to the creation of new
extrema in the collision, namely at $r=-0.64$. In summary, these stability
experiments largely confirm the proposed hypothesis for explaining Fig. 12:
collisions yield bound states when these exist as stable states, and otherwise
lead to a change in the number of extrema (typically creation of new extrema).
The only deviations from this paradigm are observed when stable bound states
exist but have a small basin of attraction. Then the finite perturbation
provided by a collision can trigger a change in the number of extrema.

Figure 15: The $L_{2}$ norm vs. control parameter $r$ for two isolas at
$\epsilon=0$, obtained by numerical continuation from traveling bound states
in scenario A with $\epsilon=0.03$ resulting from collisions at $r=-0.645$
(large isola) and $r=-0.68$ (small isola), respectively. Color-matched
solution profiles are shown in the top panel, where the red solutions (right),
corresponding to the small isola, feature a separation between the two single-
pulse LSs shorter by one wavelength as compared with the blue profiles (left),
which correspond to the large isola.

The isolas corresponding to multi-pulse bound states at $\epsilon>0$ can take
other forms as well. Figures 16 and 17 show two such cases, in which a complex
tangled branch terminates at either end in bifurcation points located at the
center of a spiral structure. These figures correspond to the continuation of
bound states resulting from scenario A collisions at $r=-0.645$ and $r=-0.7$,
respectively. As mentioned in Sec. IV.1, the latter is in fact the result of
two separate collisions. The first collision occurs between two moving LSs and
leaves one stationary and symmetric LS and one moving LS as shown in Fig. 7.
The second collision (not shown) occurs after this, when the moving LS
collides with the stationary LS from the other side due to periodic boundary
conditions, forming a traveling bound state.

Figure 16: Scenario A bound state isola continued from $r=-0.645$ and ending
in T-points at either end. Solution profiles and parts of the spiral branch
are shown at increasing levels of magnification in the upper panels. Figure
17: Scenario A bound state isola continued from $r=-0.7$ and ending in
T-points at either end. Solution profiles and parts of the spiral branch are
shown at increasing levels of magnification in the upper panels.

The bifurcations identified in Figs. 16 and 17 are of codimension two and
occur as the separation between the two LSs in the bound state grows without
limit, while the structures independently change their amplitude and degree of
asymmetry such that their velocities remain the same. As the branch spirals
towards the center, the separation between the constituent LSs grows. These
bifurcations are highly reminiscent of T-point bifurcations, also known as
Bykov points [28, 29], which have been described as corresponding to an
equilibrium-to-equilibrium heteroclinic cycle, with two branches of primary
homoclinic orbits bifurcating from it. Bifurcations of this type have been
observed in a variety of systems, including the Lorenz63 system [30], a model
of calcium pulses in pancreatic cells [31], and in type-I excitable media
[32]. In contrast, here the two constituent LSs individually approximate a
homoclinic orbit to the trivial state $u=0$ and at the T-point the state of
the system corresponds to a double homoclinic cycle. Of course, in a finite
domain we cannot reach this T-point, and our numerical continuation results
therefore only produce a periodic structure that approximates this double
homoclinic cycle. To the best of our knowledge, a T-point bifurcation from
such a double homoclinic cycle described here has not been observed
previously. Note that such T-points cannot occur for stationary LSs.

##  V Reduced model of interacting localized structures

An important concept in the analysis of LSs is the notion of spatial dynamics
that applies equally well in the comoving frame. In the simplest case one
linearizes the equation of motion (1) about the trivial state, and considers
solutions of the form $u\propto e^{\lambda x}$, cf. [12, 9]. The roots
$\lambda$ of the resulting characteristic polynomial, known as spatial
eigenvalues, govern the behavior of the system in the vicinity of the trivial
solution $u\equiv 0$. They are given by

$\lambda=\pm\sqrt{\pm\sqrt{r}-1}.$ (12)

In the case of interest, $r<0$ so that LSs exist, $\lambda$ is complex:
$\lambda=\pm\alpha\pm i\beta$ with $\alpha,\beta>0$. The (positive) real part
of $\lambda$ describes the exponential growth of $u$ away from $u=0$ and
towards LS on a $1/\alpha$ spatial scale, while the nonzero imaginary part of
$\lambda$ implies that this growth is not monotonic, but oscillatory, with
wavelength $2\pi/\beta$. For stationary LS the profile decays towards $u=0$ in
the same manner.

Interactions between LSs will be mediated by these exponentially
growing/decaying oscillatory tails, unless the proximate (i.e. most strongly
interacting) large-amplitude extrema of the two respective structures come
close enough for nonlinear effects to become relevant. Such interactions via
oscillating tails are found in a wide range of problems and have been studied
extensively [33, 34, 35, 36, 37]. Here, in a spirit similar to these works, in
particular [33], we propose a simple reduced model to quantitatively describe
these interactions.

Consider two LSs, with the closest order one amplitude extrema of either
structure located at $x_{i}$, $i=1,2$. Each LS is characterized by a drift
velocity $c_{i}$, $i=1,2$, known from the analysis in section III. Note that
$c_{i}$ is zero for symmetric LS, and nonzero for asymmetric LS. The proposed
reduced model equations, reminiscent of overdamped particle dynamics, read

$\displaystyle\frac{dx_{1}}{dt}$
$\displaystyle=c_{1}+g_{1}\cos\left(\beta|x_{1}-x_{2}|-\phi\right)e^{-\alpha|x_{1}-x_{2}|}$
$\displaystyle\frac{dx_{2}}{dt}$
$\displaystyle=c_{2}+g_{2}\cos\left(\beta|x_{1}-x_{2}|-\phi\right)e^{-\alpha|x_{1}-x_{2}|}.$
(13)

Equations (13) are integrated using a fourth-order Runge-Kutta method with
initial conditions corresponding to DNS data for any given collision. The
model involves three unknown parameters describing the interaction: two
amplitudes $g_{1},g_{2}$ and a phase $\phi$. The values of these parameters
will depend on $r$ and on the colliding structures in question.

Below, we compare the predictions of the reduced model given by Eq. (13) with
high-resolution DNS results using a gradient descent optimization to determine
$g_{1}$, $g_{2}$ and $\phi$, as described in appendix A. To conveniently
assess the quantitative agreement between the particle-like trajectories, we
consider the deviation of the relative distance between proximate extrema from
free propagation. For $x_{2}>x_{1}$ (without loss of generality), we define

$\chi(t)\equiv x_{2}(t)-x_{1}(t)-(c_{2}-c_{1})t-[x_{2}(0)-x_{1}(0)].$ (14)

###  V.1 Comparison between reduced model and DNS

####  V.1.1 Chasing: scenario A

Figure 18 shows a comparison between the DNS results (contour plot) and the
model predictions (dashed yellow lines) for a chasing collision of two
asymmetric extrema (scenario A). There is good agreement between the extrema
trajectories in both cases, including at a quantitative level, as can be seen
from Fig. 19.

Figure 18: Overlay of the reduced model trajectories on top of DNS data for
scenario A. Colored contour plot shows the DNS results at $\epsilon=0.03$,
$r=-0.68$ in scenario $A$ (collision between two asymmetric structures that
are two and three wavelengths long, respectively). Yellow dashes indicate
trajectories predicted by the reduced model (13) with parameters
$g_{1}=0.53489982$, $g_{2}=-0.30631508$ and $\phi=-1.08804942$, as well as
$c_{1}=0.0020$ and $c_{2}=0.00128276$. See also Fig. 19.  Figure 19: Residual
distance $\chi$ between proximate extrema (defined in Eq. (13)) versus time
from the DNS and the reduced model in scenario A, corresponding to Fig. 18. A
good quantiative agreement is observed between DNS and model. The interaction
is effectively repulsive, cf. section V.2.

Figure 20 shows the model parameters obtained by the gradient descent method
described in appendix A at different values of $r$ for the chasing collision
shown above. In all cases shown here, a pure bound state is formed in the
collision, without the addition or deletion of any extrema. The parameter
values are seen to depend only weakly on $r$, reflecting small changes in
spatial structure as $r$ is varied for a given type of LS.

Figure 20: Parameters of the educed model obtained for scenario A and
different values of the control parameter $r$. Red diamonds: $g_{1}$, blue
circles: $g_{2}$, green triangles $\phi$. The parameters depend only weakly on
$r$, reflecting minor changes in profile with $r$. All cases shown correspond
to collisions generating pure bound states, except $r=-0.65$, where a
metastable bound state forms (Fig. 10), leading to the appearance of
additional extrema at late times. In that special case $r=-0.65$, the fit was
performed with a cut-off time $t_{f}$ (cf. appendix A) after the initial
collision but before the generation of the additional extrema. In the range
$-0.67\lesssim r\lesssim-0.65$, and at $r\gtrsim-0.64$, the scenario A
collisions do not result in bound states, but instead generate new peaks (cf.
Fig. 12).

####  V.1.2 Symmetric-asymmetric collision: scenario B

Figure 21 shows an overlay of the reduced model trajectories on top of the DNS
results, similar to Fig. 18, but for a collision between a symmetric LS with a
maximum at its center and a two-wavelength asymmetric LS (scenario B). There
is again good agreement between the trajectories of the closest extrema in
both cases, including at a quantitative level, as one can see in terms of the
quantity $\chi$ from Fig. 22.

Figure 21: Colored contour plot of the DNS of Eq. (1) at $\epsilon=0.03$,
$r=-0.69$ in scenario $B$ (collision between symmetric and asymmetric LSs,
with extrema of opposite sign facing each other). Yellow dashes indicate the
trajectories predicted by the reduced model (13) with parameters
$g_{1}=0.3539801$, $g_{2}=-0.51798143$, and $\phi=-1.038100556$, as well as
$c_{1}=0$ and $c_{2}=-0.0019615$. See also Fig. 22.  Figure 22: Residual
distance $\chi$ between proximate extrema (defined in Eq. (13)) versus time
from the DNS and the reduced model for scenario B, corresponding to Fig. 21.
The zoom shows the equilibrium separation. The interaction is effectively
attractive, cf. section V.2.

####  V.1.3 Flipped symmetric-asymmetric collision: scenario C

Figure 23 shows an overlay of the reduced model trajectories on top of the DNS
results, similar to Fig. 18 for a collision from scenario C. There is good
agreement between the trajectories of the proximate extrema in both cases,
including at a quantitative level, as one can see in terms of the quantity
$\chi$ from Fig. 24.

Figure 23: Colored contour plot of the DNS of Eq. (1) at $\epsilon=0.03$,
$r=-0.645$ in scenario $C$ (collision between symmetric and asymmetric LSs,
with extrema of opposite sign facing each other). Yellow dashes indicate the
trajectories predicted by the reduced model (13) with parameters
$g_{1}=-0.43756$, $g_{2}=0.55907$, and $\phi=-1.013089$, as well as $c_{1}=0$
and $c_{2}=-0.001895$. See also Fig. 24.  Figure 24: Residual distance $\chi$
between proximate extrema (defined in Eq. (13)) versus time from the DNS and
the reduced model for scenario C, corresponding to Fig. 23. The interaction is
effectively attractive, cf. section V.2.

####  V.1.4 Head-on collision: scenario D

Figure 25 shows an overlay of the reduced model trajectories on top of the DNS
results for a collision from scenario D. While the trajectories in Fig. 25
look qualitatively consistent, the quantitative perspective provided by Fig.
26 reveals that the reduced model with optimal parameters clearly deviates
from the DNS results in this case. The gradient descent method was applied
with $t_{f}=5~{}500$ (see appendix A) to find the set of parameters used in
Fig. 26, but no better agreement was found for other choices of $t_{f}$, which
were also tested.

The deviation between the reduced model trajectories and the DNS result is in
agreement with expectations. Since the reduced model is exclusively built on
the linear structure of SH35, while the creation of new extrema is a highly
nonlinear process, the reduced model is not expected to capture the dynamics
correctly once nonlinear interactions become dominant. At earlier times, the
values of $\chi$ produced by the reduced model remain close to the DNS data,
and only deviate near the time of extremum creation. This highlights the
limitations of the otherwise very successful reduced model analysed here.

Figure 25: Colored contour plot of DNS of Eq. (1) at $\epsilon=0.03$,
$r=-0.65$ in scenario $D$ (head-on collision of two identical structures).
Yellow dashes indicate trajectories predicted by the reduced model (13), with
parameters $g_{1}=-g_{2}=-0.551$, $\phi=-1.90324$, as well as
$c_{1}=-c_{2}=0.0019291$. See also Fig. 26.  Figure 26: Residual distance
$\chi$ between proximate extrema (defined in Eq. (13)) versus time from the
DNS and the reduced model, corresponding to Fig. 25. Inset shows zoom on early
times. The vertical dashed line indicates the time $t\approx 5500$ when new
extrema appear. The best-fit model prediction deviates from the DNS slightly
before this, as a consequence of nonlinear effects not captured by the model.

###  V.2 Sign of the interaction: attractive or repulsive?

It is interesting to note that in all the cases described above, when the
interacting extrema were of opposite signs (scenarios A, B), we found
$g_{1}>0$ and $g_{2}<0$ (for $-\pi<\phi<0$). Conversely, when the two
interacting extrema were of the same sign (scenarios C, D), we systematically
found $g_{1}<0$ and $g_{2}>0$ (for $-\pi<\phi<0$). However, it is nontrivial
to deduce what this implies for the sign of the effective interaction, due to
its oscillatory nature. We can define an effective sign of the interaction as
follows. First we introduce the equilibrium distance $d_{eq}$ as the value of
$|x_{2}-x_{1}|$ in the bound state after the collision has occurred. From Eq.
(13), we see that equilibrium implies

$\cos\left(\beta d+\phi\right)\exp(-\alpha
d)=\frac{c_{2}-c_{1}}{g_{1}-g_{2}}.$ (15)

In general, Eq. (15) has multiple solutions, as illustrated in Fig. 27. One
may surmise that these multiple solutions are related to the overlapping
isolas shown in Fig. 15, which differ in their equilibrium distance by one
wavelength. When $c_{2}=c_{1}$, Eq. (15) has infinitely many solutions,
similar to the family of bound two-pulse states described in [26] that differ
only in their equilibrium separation. By contrast, for $c_{2}-c_{1}\neq 0$,
the number of solutions is finite. In a generic collision at $\epsilon>0$,
since the colliding LSs approach from a large distance, we expect that the
physical value of $d_{eq}$ is given by the largest positive value of $d$ that
solves Eq. (15).

Figure 27: Illustration of how the equilibrium distance $d_{eq}$ is determined
from Eq. (15).

We assume that $d_{eq}$ is known for a given collision of two LSs starting out
at a distance $|x_{2}-x_{1}|=d_{0}$, with $c_{1}>c_{2}$ (which applies to all
examples shown above). In the absence of interactions, the free propagation of
the two LSs will reduce the distance from $d_{0}$ to $d_{eq}$ in a time

$t_{free}=\frac{d_{0}-d_{eq}}{c_{1}-c_{2}}.$ (16)

In the presence of interactions, this time will change to $t_{eq}$, a time
that may be larger or smaller than $t_{free}$ due to the oscillatory nature of
the interaction.

We propose the following terminology. If the equilibrium distance in the
collision is reached earlier than for free propagation, i.e. if
$t_{eq}<t_{free}$, then we say that the interaction is effectively attractive.
By contrast, if the equilibrium distance is reached later, i.e.
$t_{eq}>t_{free}$, then we say that the interaction is effectively repulsive.
We note that the effective sign of the interaction can be determined
graphically from the sign of $\chi(t_{eq})$, since Eq. (14) implies

$t_{eq}=t_{free}+\frac{\chi(t_{eq})}{|c_{1}-c_{2}|}.$ (17)

The time $t=t_{eq}$, where the distance $|x_{2}-x_{1}|=d_{eq}$ is highlighted
by an arrow in Figs. 19, 22, 24 and 26; since $x_{2}-x_{1}$ is constant at
$t>t_{eq}$, a linear increase in $\chi=const+(c_{1}-c_{2})t$ sets in. The sign
of $\chi(t_{eq})$, i.e. the onset of this linear increase, determines the
effective sign: if $\chi(t_{eq})>0$, then the interaction is effectively
repulsive but if $\chi(t_{eq})<0$, then the interaction is effectively
attractive.

We highlight that even though the collisions from scenarios A and B shown in
Figs. 18 and 21 both feature interactions between extrema of opposite signs,
the relative signs of their interaction differ: the interaction is effectively
repulsive in scenario A (Fig. 19), while it is effectively attractive in
scenario B (Fig. 22). This indicates that while the relative sign of
interacting extrema does appear to be important, it is not the only factor
determining collision outcomes. This is likely related to the observed
deviation between scenarios A, B and C, D, respectively, in the parameter
range $-0.65\lesssim r\lesssim-0.64$.

##  VI Conclusions

In this paper, we have provided an in-depth analysis of LSs in the
nonvariational SH35, Eq. (1), using numerical continuation, DNSs, asymptotics
and reduced-order modeling to shed new light on the propagation and
interaction of LSs in a canonical driven dissipative system. These
interactions are highly inelastic, in contrast to systems described by
integrable partial differential equations, and lead to both stationary and
drifting structures. Moreover, the interactions can be attracting or repelling
depending on the nature of the interacting LSs and the parameters.

The asymptotic theory predicts a linear dependence of the drift speed of LSs
on $\epsilon$, with excellent quantitative agreement with DNSs and numerical
continuation at $\epsilon\lesssim 0.1$, but significant deviation for
$\epsilon\gtrsim 0.1$. The collisions resulting from this drift are shown to
display rich phenomenology: different numbers of extrema can be added or
deleted in a collision, depending on the types of structures interacting and
the value of $r$. Alternatively, a pure bound state can be formed, which
preserves the number of extrema of the initial structures. We have found that
the stability properties of these bound states play a key role in determining
whether a collision changes the number of extrema or not: if the bound state
exists stably, and has a sufficiently large basin of attraction, then the
perturbation resulting from a collision does not change the number of extrema.
However, if this is not the case, then extrema will be added or deleted in the
collision. When extrema are deleted, this can lead to smaller LSs, or to
annihilation in the case of the minimal, single-wavelength asymmetric
structure. When extrema are created, metastable states may arise from
collisions, which undergo a further nonlinear interaction at late times. This
phenomenology is reminiscent of what has been attributed to a scattor [24,
25]: unstable stationary or time-periodic patterns which direct orbits during
the collision process in the infinite-dimensional state space along their
stable and unstable manifolds. While these ideas were proposed in the context
of models other than SH35, they may be applicable to some of the collision
events observed here.

Bound states arising from collisions were shown to lie on isolas which are
embedded within the snakes and ladders bifurcation structure. In contrast with
the variational case $\epsilon=0$, the isolas at $\epsilon>0$ are typically
not of figure-eight form, but rather form a complex, spaghetti-like tangled
set. Only for the specific subclass of collisions between chasing asymmetric
LSs, are the resulting multi-pulse bound states found to lie on isolas whose
figure-eight form is preserved at $\epsilon>0$. In addition, we have also
described a novel example of an isola at $\epsilon>0$ which is not closed, but
features T-point bifurcations at either end, each of which involves the
gradual separation of a bound state into two separate asymmetric LSs, each
separately corresponding to a homoclinic connection to $u=0$ in the comoving
frame. States of this type require drift and so can only be found in
nongradient systems of the type studied here.

Finally, a reduced model was also proposed, consisting of two coupled ordinary
differential equations, describing the interaction of LSs via their
oscillatory exponential tails. The reduced model was shown to reproduce a wide
range of collisions, with quantitative accuracy, provided no large-amplitude
extrema are created or destroyed in the collision. If the number of large-
amplitude extrema is altered in the collision, the model still describes the
trajectories of the interacting proximate extrema up until shortly before the
time at which this occurs. After this time, the model fails to describe the
DNS results since it does not include the nonlinear structure of Eq. (1). The
fact that the trajectories of LSs can be accurately reproduced by the reduced
model for a sizable fraction of all observed cases is indicative of the fact
that the collisions between localized patterns are to a large extent particle-
like, with a nontrivial interaction potential corresponding to the linear
structure of Eq. (1). This model led to considerable insight into the nature
of the interactions between LSs in the SH35 model and allowed a relatively
simple determination of the conditions under which the interaction is
effectively attracting or repelling.

Collisions of convectons in binary fluid convection, such as those described
in [19], feature dynamics beyond the scope of the order parameter description
provided by SH35, in particular due to the possibility of complex temporal
behavior. Nonetheless, it remains to be understood to what extent the SH35
collision phenomenology described here can be observed between convectons or
similar structures in other continuum systems. Specifically, [19] do not
observe bound states, while these play an important role in the selection of
collision outcomes described here. It remains an open question whether such
bound states exist between convectons, given that binary fluid convectons
interact nonlocally via a large-scale solute field generated by the pumping
mechanism described in [4].

A simple possibility for obtaining more complex temporal behavior in SH35
would be the inclusion of higher-order time-derivatives, e.g. a second-order
term that allows for inertia and wave-like solutions. In this case, the
interactions described here, which rely exclusively on the exponential tails
of LSs, would be enriched due to the possibility of wave radiation. The study
of this modified SH35 with higher-order time derivatives is left for a future
study.

###### Acknowledgements.

We acknowledge support from the National Science Foundation under grants
DMS-1908891 and DMS-2009563. We also thank Nicolás Verschueren van Rees for
fruitful discussions and for providing us with a high performance C++ solver
which we adapted for our DNS studies.

##  Appendix A Reduced model parameter estimation by gradient descent

Figure 28: Illustration of the robustness of the gradient descent method:
starting from various different initial guesses for the model parameters the
method converges to a unique set of parameters $g_{1}$, $g_{2}$, $\phi$.
Example shown corresponds to scenario $B$ at $r=-0.69$.

In this appendix, we describe the details of the gradient descent method used
to systematically determine the reduced model parameters $g_{1}$, $g_{2}$ and
$\phi$ corresponding to SH35 solutions for a given $r$ and a given collision
scenario, based on high-resolution DNS data.

Given $u(x,t)$ from DNS on a uniform grid with $8192$ points, we track the
position of the two extrema that are located closest to one another in the
collision, as a function of time $t$ (note that a high spatial resolution is
required for well-resolved tracking). We denote the positions of the closest
extrema by $x_{1}^{DNS}(t)$, $x_{2}^{DNS}(t)$, corresponding to $x_{1}(t)$,
$x_{2}(t)$ in Eq. (13). To determine which set of reduced model parameters
$g_{1}$, $g_{2}$, $\phi$ most accurately represents the DNS results, we define
a cost function

$C_{g_{1},g_{2},\phi}[x_{1},x_{2}]\equiv\sqrt{\frac{1}{t_{f}}\int_{0}^{t_{f}}\sum_{i=1}^{2}\left(x_{i}(t)-x_{i}^{DNS}(t)\right)^{2}dt},$
(18)

where at $t=0$ the colliding LSs are separated by many wavelengths, and
$t_{f}$ can be chosen to be long after the collision, if no large-amplitude
extrema are created or destroyed (e.g. Fig. 18), or instead may be chosen just
before nonlinear effects set in (e.g. Fig. 26). The global minimum of the
functional $C_{g_{1},g_{2},\phi}[x_{1},x_{2}]$ corresponds to the best fit
between reduced model and DNS.

To perform the gradient descent optimization, we start with an initial guess
for $g_{1}$, $g_{2}$, $\phi$, and compute the gradient $\nabla C=(\partial
C/\partial g_{1},\partial C/\partial g_{2},\partial C/\partial\phi)$. Then, in
each step, we update parameters by the rule

$(g_{1},g_{2},\phi)\to(g_{1},g_{2},\phi)-\lambda\nabla C,$ (19)

with a parameter $\lambda$ measuring the step size in parameter space.

Figure 28 shows sample trajectories in parameter space (top panel: projection
onto the $g_{1},\phi$ plane, bottom panel: projection onto the $g_{2},\phi$
plane) obtained by this method for a collision from scenario $B$ (asymmetric
LS colliding with a stationary symmetric LS) at $r=-0.69$ (collision results
in pure bound state). For a wide range of initial conditions, the method
converges to a well-defined set of parameters. The resulting model
trajectories accurately reproduce those observed in the DNS, as can be seen in
Figs. 18 and 19. The significance of Fig. 28 is that it indicates that the
method is robust with respect to the precise initial guess for the model
parameters $g_{1},g_{2},\phi$. However, several complicating factors need to
be taken into account to obtain the correct optimal model trajectories with
this procedure.

  1. 1.

Due to the very long simulation times of $O(10^{4}-10^{5})$ time units, the
propagation velocities $c_{1},c_{2}$ have to be specified with high precision
in order to accurately describe the free propagation of the LSs. Only after
carefully determining $c_{1},c_{2}$ with a high precision does the gradient
descent algorithm converge to an accurate fit between model and DNS.

  2. 2.

Equations (13) are invariant under $\phi\to\phi+2\pi$ and $\phi\to\phi+\pi$,
$(g_{1},g_{2})\to(-g_{1},-g_{2})$. The global minimum of $C$ is only unique if
the phase is limited to an interval of size $\pi$. Here, we choose
$\phi\in[-\pi,0]$.

  3. 3.

The cut-off time $t_{f}$ can affect the quality of the fit: if it is too
large, the final bound state will be weighted excessively compared to the
actual collision dynamics of interest. If $t_{f}$ is too short, it may fail to
correctly capture the later stages of the collision. We typically choose
$t_{f}$ to be several thousand time units after the equilibrium distance is
reached, for pure bound states, and before the nonlinear creation of new
extrema, if any.

## References

  * [1] Diederik Johannes Korteweg and Gustav De Vries.  XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 39(240):422–443, 1895.
  * [2] Carlo R Laing, William C Troy, Boris Gutkin, and G Bard Ermentrout.  Multiple bumps in a neuronal model of working memory.  SIAM J. Appl. Math., 63(1):62–97, 2002.
  * [3] Kassem Ghorayeb and Abdelkader Mojtabi.  Double diffusive convection in a vertical rectangular cavity.  Phys. Fluids, 9(8):2339–2348, 1997.
  * [4] Oriol Batiste, Edgar Knobloch, Arantxa Alonso, and Isabel Mercader.  Spatially localized binary-fluid convection.  J. Fluid Mech., 560:149–158, 2006.
  * [5] Mohamed Gad-El-Hak, Ron F Blackwelderf, and James J Riley.  On the growth of turbulent regions in laminar boundary layers.  J. Fluid Mech., 110:73–95, 1981.
  * [6] Tobias M Schneider, Daniel Marinc, and Bruno Eckhardt.  Localized edge states nucleate turbulence in extended plane Couette cells.  J. Fluid Mech., 646:441–451, 2010.
  * [7] Pierre Coullet, C Riera, and Charles Tresser.  Stable static localized structures in one dimension.  Phys. Rev. Lett., 84(14):3069, 2000.
  * [8] Arik Yochelis, Yin Tintut, Linda L Demer, and Alan Garfinkel.  The formation of labyrinths, spots and stripe patterns in a biochemical approach to cardiovascular calcification.  New J. Phys., 10(5):055002, 2008.
  * [9] Edgar Knobloch.  Spatial localization in dissipative systems.  Annu. Rev. Condens. Matter Phys., 6(1):325–359, 2015.
  * [10] Pierre C Hohenberg and Jack B Swift.  Effects of additive noise at the onset of Rayleigh-Bénard convection.  Phys. Rev. A, 46(8):4773, 1992.
  * [11] Yi-Ping Ma and Edward A Spiegel.  A diagrammatic derivation of (convective) pattern equations.  Physica D: Nonlinear Phenomena, 240(2):150–165, 2011.
  * [12] John Burke and Edgar Knobloch.  Localized states in the generalized Swift-Hohenberg equation.  Phys. Rev. E, 73(5):056211, 2006.
  * [13] John Burke and Edgar Knobloch.  Snakes and ladders: localized states in the Swift–Hohenberg equation.  Phys. Lett. A, 360(6):681–688, 2007.
  * [14] John Burke and Edgar Knobloch.  Homoclinic snaking: structure and stability.  Chaos, 17(3):037102, 2007.
  * [15] Gregory Kozyreff and Mustapha Tlidi.  Nonvariational real Swift-Hohenberg equation for biological, chemical, and optical systems.  Chaos, 17(3):037103, 2007.
  * [16] John Burke and Jonathan HP Dawes.  Localized states in an extended Swift–Hohenberg equation.  SIAM J. Appl. Dyn. Syst., 11(1):261–284, 2012.
  * [17] John Burke, S M Houghton, and Edgar Knobloch.  Swift-Hohenberg equation with broken reflection symmetry.  Phys. Rev. E, 80(3):036202, 2009.
  * [18] SM Houghton and Edgar Knobloch.  Swift-Hohenberg equation with broken cubic-quintic nonlinearity.  Phys. Rev. E, 84(1):016204, 2011.
  * [19] Isabel Mercader, Oriol Batiste, Arantxa Alonso, and Edgar Knobloch.  Travelling convectons in binary fluid convection.  J. Fluid Mech., 722:240–266, 2013.
  * [20] Eusebius Doedel, Alan R Champneys, Thomas F Fairgrieve, Yuri A Kuznetsov, Björn Oldeman, R Paffenroth, B Sandstede, Xianjun Wang, and C Zhang.  AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations, 2012.
  * [21] Hannes Uecker.  Numerical Continuation and Bifurcation in Nonlinear PDEs.  SIAM, Philadelphia, PA, 2021.
  * [22] Elizabeth Makrides and Björn Sandstede.  Predicting the bifurcation structure of localized snaking patterns.  Physica D, 268:59–78, 2014.
  * [23] Hannes Uecker.  A short ad hoc introduction to spectral methods for parabolic PDE and the Navier-Stokes equations.  Summer School Modern Computational Science, pages 169–209, 2009\.
  * [24] Yasumasa Nishiura, Takashi Teramoto, and Kei-Ichi Ueda.  Dynamic transitions through scattors in dissipative systems.  Chaos, 13(3):962–972, 2003.
  * [25] Yasumasa Nishiura, Takashi Teramoto, and Kei-Ichi Ueda.  Scattering of traveling spots in dissipative systems.  Chaos, 15(4):047509, 2005.
  * [26] John Burke and Edgar Knobloch.  Multipulse states in the Swift-Hohenberg equation.  In Conference Publications, volume 2009, page 109. American Institute of Mathematical Sciences, 2009.
  * [27] Margaret Beck, Jürgen Knobloch, David JB Lloyd, Björn Sandstede, and Thomas Wagenknecht.  Snakes, ladders, and isolas of localized patterns.  SIAM J. Math. Anal., 41(3):936–972, 2009.
  * [28] VV Bykov.  Bifurcations of dynamical systems close to systems with a separatrix contour containing a saddle-focus.  Methods of the Qualitative Theory of Differential Equations, pages 44–72, 1980.
  * [29] VV Bykov.  The bifurcations of separatrix contours and chaos.  Physica D, 62(1-4):290–299, 1993.
  * [30] Paul Glendinning and Colin Sparrow.  T-points: a codimension two heteroclinic bifurcation.  J. Stat. Phys., 43:479–488, 1986.
  * [31] Mónica M. Romeo and Christopher K. R. T. Jones.  The stability of traveling calcium pulses in a pancreatic acinar cell.  Physica D, 177(1):242–258, 2003.
  * [32] Pablo Moreno-Spiegelberg, Andreu Arinyo-i Prats, Daniel Ruiz-Reynés, Manuel A Matias, and Damià Gomila.  Bifurcation structure of traveling pulses in type-I excitable media.  Phys. Rev. E, 106(3):034206, 2022.
  * [33] P Coullet, C Elphick, and D Repaux.  Nature of spatial chaos.  Phys. Rev. Lett., 58(5):431, 1987.
  * [34] Christian Elphick, Glenn R Ierley, Oded Regev, and Edward A Spiegel.  Interacting localized structures with galilean invariance.  Phys. Rev. A, 44(2):1110–1122, 1991.
  * [35] Shin-Ichiro Ei and Takao Ohta.  Equation of motion for interacting pulses.  Phys. Rev. E, 50(6):4672–4678, 1994.
  * [36] Neil J Balmforth, Glenn R Ierley, and Edward A Spiegel.  Chaotic pulse trains.  SIAM J. Appl. Math., 54(5):1291–1334, 1994.
  * [37] Yasumasa Nishiura and Takeshi Watanabe.  Traveling pulses with oscillatory tails, figure-eight-like stack of isolas, and dynamics in heterogeneous media.  Physica D, 440:133448, 2022.