). Discrete snaking: Multiple Cavity Solitons in Saturable Media.

. This paper continues an investigation into a one-dimensional lattice equation that models the light ﬁeld in a system comprised of a periodic array of pumped optical cavities with saturable nonlinearity. The additional eﬀects of a spatial gradient of the phase of the pump ﬁeld are studied, which in the presence of loss terms is shown to break the spatial reversibility of the steady problem. Unlike for continuum systems, small symmetry-breaking is argued to not lead directly to moving solitons, but there remains a pinning region in which there are inﬁnitely many distinct stable stationary solitons of arbitrarily large width. These solitons are no-longer arranged in a homoclinic snaking bifurcation diagrams, but instead break up into discrete isolas. For large enough symmetry-breaking, the fold bifurcations of the lowest intensity solitons no longer overlap, which is argued to be the trigger point of moving localised structures. Due to the dissipative nature of the problem, any radiation shed by these structures is damped and so they appear to be true attractors. Careful direct numerical simulations reveal that branches of the moving soli-tons undergo unsual hysteresis with respect to the pump, for suﬃciently large symmetry breaking.

1. Introduction.Recently, the present authors [33,34] studied the rich variety of stable localised structures that can occur in a spatially discrete model for an optical cavity with an imposed periodic structure and saturable nonlinearity.In particular, it was shown how infinitely many distinct stable localised structures can occur for a wide range of parameter values due to the so-called homoclinic snaking mechanism that has received a lot of attention in continuum models of pattern formation [1,4,5,6,8,14,20,22,35].This paper explores the consequences of introducing a spatial gradient to the optical pumping field (the forcing term in the model), which breaks the spatial reversibility of the steady system.
The model can be written in dimensionless form as the one-dimensional lattice equation for a complex field A n ∈ C, n ∈ Z: The variable A n represents the amplitude of the field inside the n-th identical optical resonator in a one-dimensional array; see [33], which generalises earlier work of Egorov et al. [9,11,29] for related models.The coefficient c ≥ 0 represents the strength of the nearest-neighbour coupling between oscillators due to their evanescent field.The limit c → ∞ corresponds to the continuum limit, whereas c → 0 represents the so-called anti-continuum limit in which the oscillators are independent of each other.The real part of δ represents detuning of the pump frequency from the resonant frequency of the oscillators.The parameter α represents the strength of the Kerr effect of intensity-dependent refractive index, where the physically realistic effect of saturation for large intensity has been included.The field strength has been scaled so that the saturation constant (the coefficient multiplying |A| 2 in the denominator of the nonlinear term of (1)) has been set to unity.The parameter P measures the amplitude of an applied optical pump field whose phase we suppose to be a simple periodic function of space, with period 2π/q.The case q = 0, where there is no phase gradient of the pump, reduces (1) to the model studied in [33,34].
For this work we restrict attention to the case of a focusing nonlinearity α > 0, but in general allow δ to be complex, with Im(δ) representing a linear loss term.We consider only passive media which implies Im(δ) ≥ 0.
In order to seek stationary solitons (in what follows the term "soliton" will be used to mean simply a localised solution), it is helpful to view the equilibrium equation for stationary solutions of (1) as a spatial dynamical system in the discrete "evolution" variable n ∈ Z: Soliton solutions to (1) thus correspond to homoclinic orbits of the 2D map (2), (3) for (A n , B n ) ∈ C 2 , where the additional variable B n = A n − A n−1 measures the slope of the soliton at a point between lattice sites n and n − 1.
In [33] for the case q = 0 it was argued that one can expect two possible kinds of homoclinic solution to (2), (3) specifically, those that are invariant under the two separate reversibilities and reversal of the direction of iteration.Introduction of non-zero q breaks the invariance of ( 2), (3) under both of these reversibilities.Note however in the absence of linear loss terms, Im(δ) = 0, there are two different reversibilities in which * represents complex conjugation, under which the map is invariant with reversal of the direction of iteration.
If both q and Im(δ) are nonzero then the spatial dynamical system is no longer reversible and any homoclinic solutions will be asymmetric in general.Conversely if q = Im(δ) = 0 then the two sets of reversibilities Indeed, in that case the equations admit the invariant reduction in which A n (and hence B n ) can be assumed to be real.Such localised solutions are thus homoclinic solutions of a two-dimensional real map.The main result in [33] was a systematic investigation of parameter regions for α, c, P and δ in which bistability occurs between steady spatially homogeneous steady states and for which there is a Maxwell point between two of these states.(In this context, where there can be temporal dissipation, we define a Maxwell point to be where a heteroclinic connection exists in the spatial dynamics).It was then shown how the pinning regions [30] near such Maxwell points can lead to either bright solitons (with intensity of the core state being greater than the tail) or grey solitons (with tail intensity greater than that of the core) which have arbitrarily wide spatial extent.These solutions are organised in a homoclinic snaking bifurcation diagram as a single parameter (typically the pump) is varied through the pinning region.
In this work we shall explore the effect of adding the phase gradient q to these localised structures.Rather than a complete exploration of parameter space, we shall instead restrict attention to an explanation of the new effects that this symmetrybreaking term brings.Therefore, in what follows, we shall fix the values c = 1.0, α = 10, and Re(δ) = −9.2 and consider only variation of the symmetry-breaking parameters q and Im(δ) in addition to the primary bifurcation parameter P .Also, we shall restrict attention exclusively to bright solitons.The next section considers stationary localised structures, which are computed using the same methods as in [33], that is numerical continuation of localised solutions to (2), (3) along with spectral computations and direct simulations of (1) to analyse stability.Section 3 then goes on to study possible moving solitons, using direct numerical simulation only.Finally, Section 4 draws conclusions and discusses the implications of the results for other work.
2. Stationary solitons.The special case q = Im(δ) = 0 corresponds to a conservative system.Figure 1(a) shows the bifurcation diagram for on-site symmetric, off-site symmetric and asymmetric solitons in this case, together with the spatially homogeneous states.All solutions depicted have the same single coherent phase as the pump, that is A n is real for all n.Note that stable solitons exist throughout the interval of P -values for which the low-intensity homogeneous state is stable.For low pump, these stable solitons are on-site symmetric, whereas for higher P it is the off-site solitons that are stable.There is also an interval of bistability between the two kinds of soliton.This interval of bistability is bridged by a branch of asymmetric solitons that are neither site centred nor off-site centred.In fact, there are two such branches overlaid in the bifurcation diagram, which are left-right mirror images of each other.These solutions bifurcate from the symmetric solitons in pitchfork bifurcations at precisely the same parameter values at which stability is lost.
Panels (b), (c) and (d) of Fig. 1 illustrate by direct numerical simulation the different mechanisms of instability of the various unstable soliton branches.Note that small perturbations (due to numerical approximation) in all cases lead to an initial small translation of the soliton.For the unstable on-site solitons this causes an oscillation in both the position and the amplitude of the soliton core.Eventually this oscillation becomes unstable and the entire domain becomes is attracted to oscillatory states that are close to the high-intensity homogeneous state.The offsite soliton on the other hand seems to decay as it oscillates and eventually the low-intensity background state fills the whole domain.The asymmetric soliton also spread out to a high-intensity state, but the initial meander of the soliton core has a definite one-sided bias to it.Note that we have found no evidence of there being fundamentally different outcomes if definite perturbations are made that either add or subtract energy from the solitons.Figure 2 shows the effect on the soliton branches in Fig. 1 of introducing non-zero q and Im(δ), both separately and together.Panel (a) shows the effect of introducing small nonzero q.Note that now, strictly speaking, there are now no longer any spatially homogeneous states.However there are stationary states that have constant |A| and a spatially periodic phase that matches that of the pump field.For simplicity we shall still refer to these as spatially homogeneous states.Since the steady equations ( 2),(3) still admit two reversibilities, R * 1 and R * 2 , we find that both on-site and off-site symmetric solitons still occur in this case, although now A n is in general complex for all n.The bifurcation diagram of these solitons remains qualitatively similar to that in Fig. 1 in that there are symmetry-breaking pitchfork bifurcations of both the on-site and the off-site branches.These bifurcations give birth to a branch of asymmetric solitons that connects the two different types of symmetric soliton.The one new feature is that the introduction of the phase gradient creates a new window of instability on the branch of off-site solitons.This Figure 2(c) shows the effect of introducing a small linear loss term Im(δ) = 0.05 to the bifurcation diagram in Fig. 1(a), while keeping q = 0.In this case the main discernible difference with the conservative case is that neither the excited homogeneous state nor the solitons reach all the way to zero pump values.Thus there is a lower threshold on P as well as an upper one for which soliton solutions can be observed.As in the case with phase gradient alone, A n is once again complex for all n for all of the states, because of the non-zero imaginary part of the steady equations.Now, if we introduce both loss and phase gradient terms, then as discussed in the Introduction, all reversibilities of the spatial dynamics are destroyed, and so we expect all soliton branches to be asymmetric.Figure 2(d) illustrates what happens for small loss and gradient terms.The pitchfork bifurcations are immediately destroyed and are unfolded in the usual way into a single fold bifurcation and continuously connected branch.The two branches that were spawned at the pitchforks are no longer exactly overlaid.Also, the fold bifurcation of the low-intensity homogeneous solution (at parameter value P = P + = 4.255195 for the parameter values used here) can no longer be described by the normal form of a reversible saddle-node bifurcation [21] and so there is no small amplitude bifurcation of the pair of off-site and on-site solitons.Instead, the small-amplitude portions of these two branches (which are both unstable) now meet each other in a fold.Thus the bifurcation structure has split into two separate pieces.
Figure 3 shows what happens to the higher-intensity disconnected piece of this bifurcation diagram as either the loss or the gradient are further increased.As explained in more detail in [33], the parameter values (4) are close to those for which a codimension-two point two soliton branches in Fig. 1(a) connect in a further fold at precisely P = P + .In fact, this transition happens as either q > 0 or δ > 0 are increased, specifically for q ≈ 0.05 and δ ≈ 0.4; see Fig. 3. Beyond these transition values of either parameter, the the right-hand end right-hand end of the soliton Blow up of the small region of panel (b) in which the former thick solid and dashed curves are now depicted as thin black lines, with the addition of the asymmetric states that are spawned in pitchfork bifurcations from these on-site and off-site symmetric solitons.The solid and dashed lines in this panel represent the corresponding curves of solitons that occur after the inclusion of a small phase gradient q = 0.01.(d) Graphs of the stable asymmetric solitons at the points indicated in panel (c).branch closes up so that the entire curve becomes an isola.This isola continues to decrease in size as either the loss or the gradient is further increased, until a further critical value at which the isola disappears altogether by shrinking to a point.
In [33] it was shown for higher losses and no phase gradient that the bifurcation diagram in Fig. 1 turns into a so-called homoclinic snaking curve.See Fig. 4(a),(b) for the case Im(δ) = 1.This occurs because as Im(δ) is increased through a critical value ≈ 0.75, a heteroclinic connection occurs between the high-intensity state and the low intensity one for precisely P = P + .As Im(δ) is further increased, this heteroclinic connection moves into the P -interval of bistability between the upper and lower homogeneous states and this causes the curves of solitons to snake around this Maxwell point.Figure 4(c) shows what happens to this snaking curve when the reversibility is completely broken by adding a small gradient term q > 0. The thin lines in the diagram shows a detail of the zero-q snake, including how the (pairs of) asymmetric solitons form 'rungs' of the 'snakes and ladders' bifurcation structure described by Burke and E.Knobloch [5,4] for the continuum Swift-Hohenberg equation.The asymmetric branches connect the on-site and offsite solitons.This is also in accord with the general theory of homoclinic snakes derived by Beck et al [1].The rungs are born in pitchfork bifurcations that occur very close to the parameter values of the turning points (folds) of the symmetric orbits.
That the bifurcation diagram must split like this when symmetry is broken can be explained approximately by considering how the stable and unstable manifolds unfold in the neighbourhood of a turning point of a homoclinic snaking curve.The sketch in Fig. 5 represents the behaviour of intersections between stable and unstable manifolds in a two-dimensional reversible map. Figure 6  picture when reversibility is broken.Note that the real and imaginary parts of (2), (3) effectively represent a four-dimensional real map for the parameter values in question.But, since the stable and unstable manifolds are both two-dimensional, one can depict the dynamics by taking a cross-section that ignores one dimension in each of the stable and unstable directions (in essence appealing to the "homoclinic centre manifold theorem" [31]).Hence a two-dimensional diagram can be used to capture the necessary topology.For slightly higher q-values than those depicted in Fig. 4(c), the isolas each shrink to a point and the infinite multi-stability of stationary localised structures is lost.Specifically Fig. 7 shows what happens to the highest mass isola in Fig. 4(c) as the phase gradient is further increased.Note that the formation and destruction of isolas as one breaks the symmetry is qualitatively similar to that observed in continuum systems when the reversibility is destroyed by the introduction of a third-order dispersion term in the Swift-Hohenberg equation [2].The difference with that case, though is that there the isolas represent travelling solitary structures, whereas here they are stationary.Other relevant results for continuum systems is the analysis in [19,3] in which two-pulse stationary solitons near a homoclinic snake are shown to generically lie on isolas with a similar topology to those encountered here, even in the absence of symmetry-breaking.
3. Moving solitons.It is well-known that solitons of 1D discrete-NLS type lattices, unlike their continuum counterparts, cannot move with arbitrarily small wave speed due to the so-called Peierls-Nabarro energy barrier.For such models in the absence of dissipation, the on-site and off-site solitons typically have different values of the Hamiltonian energy of the system.Thus if an infinitesimal motion is given Figure 7. Zoom of the upper isola in Fig. 4(c) and its transformation as q is further increased.
to the soliton, it must repeatedly climb and descend energy barriers as it passes through the location of the stationary on-site and off-site solitons.This process causes the shedding of radiation and the ultimate destruction of the soliton (see, e.g.[17] and references therein).Nevertheless, in recent studies it has been shown that special nonlinearities can lead to so-called transparent points at which the on-site and off-site solitons have the same energy, see [28] and references therein.Thus, there can be values of parameters at which exact travelling solitons that do not shed radiation can be shown to occur [28,25,26].Albeit, theses solitons cannot occur for arbitrarily small wave speeds, due to additional resonances with phonons [24,23].
For models with forcing and damping as in (1), there is no analogue of the Hamiltonian that will create a Peierls-Nabarro barrier and so one might expect solitons to move as soon as a phase gradient is added.However, for a related model Egorov et al [10] proposed a different argument to show why one should not expect immediate onset of moving localised structures.Let us explain this argument in the context of the present model.At the parameter values of interest, e.g. in Fig 2(c) for 1 < P < 3, there is at least one stable stationary soliton.So the introduction of infinitesimal motion to a soliton is likely to make it become attracted to this stationary state.To overcome this attraction, a sufficiently large phase gradient q > q p needs to be introduced before one might expect to see moving structures.In the asymptotic limit of large c, Egorov et al [10] provide an asymptotic description of this 'pinning region' q < q p for which solitons will remain stationary.
Here we should like to propose another criterion for determining when solitons should start to move.Intuitively, there would be no barrier if at a critical value of q, the translational mode of the two stationary solitons each had a zero eigenvalue.Actually we can find such a condition precisely at a point in the (P, q) parameter plane at which both branches of stationary solitons simultaneous undergo a fold with respect to P .This we should expect to be the onset point of a branch of moving solitons.The numerical results in Fig. 8(a),(b) show that this criterion is indeed true.For q-values just higher than this 'fold-disconnect' value, numerical simulations show what appears to be a branch of travelling coherent localised structures whose wave speeds tend to zero at the P -value of either of the folds.For higher loss values, as the two folds occur for more widely spread P -values, as in Fig. 8(c),(d),(e), we find a wider branch of these 'moving solitons', which have higher approximate wave speeds.
Note from Figure 8(a),(c) that the maximum of the moving soliton amplitude has a discontinuity at the point at which the wave speed reaches zero.This is because the moving solitons appear actually to be periodic solutions, in which the maximum value of the soliton oscillates.Hence the peak of the oscillation has a value that is greater than the maximum amplitude of the field of a stationary soliton at the same parameter values.
The dynamics of the moving solitons is much richer if the gradient is higher and/or the loss is lower, as illustrated in Fig. 9 which has double the value of q from those used in Fig. 8(c),(d),(e) and a quarter of the value of Im(δ).Panel (a) shows the bifurcation diagram.Note that in the vicinity of the folding points the stationary solitons are now unstable.The instability leads to the formation of a stable moving soliton.The development of the instability is shown in panel (b) for parameter values near the fold.The only perturbation added is due to numerical error, so the dynamics remains close to the stationary state for a long time before eventually a weak exponential instability sets in leading to the emission of radiation and the one-sdied drift of the core of the soliton.The middle panel is plot on a logarithmic scale to make visible radiative tail that appears behind the moving soliton.At these parameter values the radiation is rather weak and it decays rapidly because the loss is relatively high.Note that the core of the soliton quickly establishes a definite steady velocity (the measurement of which is how we calculate v in Figs. 8 and 9).As it moves, the soliton continues to emit raditation, but in a regular manner that is compensated by the pumped internal dynamics of the soliton, so that the entire structure is approximately a coherent travelling wave.Figure 9(c) shows the dependence of the soliton velocity on the pump.Note the parameter regions of hysteresis where both stable stationary and moving solitons occur.
To explain the hysteresis further let us describe what happens when one decreases P quasi-statically starting from a stable soliton belonging to the right-hand branch of stationary solitons.The soliton will initially follow the bifurcation curve of the stable stationary solitons until it reaches the point where the soliton becomes unstable.At this point there is a large dynamic jump in the amplitude and velocity of the observed long-time dynamics, to a moving soliton.As P is further decreased, the soliton continues to move with finite velocity, even when P is reduced to the point at which the stationary solitons become stable.
For P ≈ 1 several sudden jumps are observed as one continues to decreases the pump.A blow-up of this situation is depicted in Figure 10(a).The first jump is to a soliton moving with different lower velocity and having different intensity.There is at least one further jump to a soliton moving with yet lower speed, and the final jump is to the branch of stationary solitons.The dynamics of the three distinct moving solitons are shown indicated in Fig. 10(b)-(d).Note there is again hysteresis, in that if start from the stationary soliton at P = 1 and increase the pump, then the soliton only starts to move at the fold point where the stationary soliton becomes unstable.There also seems to be a similar 'staircase' structure to the hysteresis loop at the high-intensity end of the bifurcation diagram in Fig. 9(c) although we have not investigated this in detail.

Conclusion.
This paper has thrown up a number of interesting features that occur due to a combination of the discreteness of the model studied, and the addition of terms that break the symmetry.First we showed that the homoclinic snake is only broken by the inclusion of both loss terms and a phase gradient.One effect alone does not break the spatial reversibility of the underlying discrete map, and so the reversible theory of homoclinic snaking [1,6,35] still applies.Breaking the symmetry in the discrete system studied here has similarities with what happens in continuum systems that possess homoclinic snakes, see [2].In that case, one also sees isolas that are rather like those in Figure 4.However, in the continuum case, the states that form the isolas are moving solitons.Here, instead we find that the isolas represent stationary states.This can be understood because of the fundamental difference between discrete-time and continuous-time dynamical systems in the absence of reversibility or Hamiltonian structure.As is well known, for continuous-time systems the existence of homoclinic orbit to a saddle point equilibrium is of codimension one.That is, in order to see a curve of homoclinic orbits in a bifurcation diagram, one has generically to vary two parameters, in this case P and the soliton's velocity v.In contrast, for a discrete-time system, the existence of a homoclinic orbit to a saddle fixed-point is a persistent phenomenon.That is, if stable and unstable manifolds intersect transversally for one parameter, then they will continue to do so for nearby parameters, and one sees a curve of homoclinic solutions as one varies a single parameter, P say.It is only in the presence of reversibility or Hamiltonian structure that continuum homoclinic orbits are also persistent phenomena.The arguments justifying the shape of the isolas on breaking the symmetry were only sketched in Figs. 5, 6. Presumably, more rigorous versions of these results may be obtained using the methodology in [1].Indeed, [19] have recently used such methods to prove that multi-pulse versions of the homoclinic orbits that snake in continuum systems lie on similar shaped isolas.
With regard to moving solitons, we have found similarities with discrete conservative systems, in agreement with the results in [10] for a related system in that there sufficiently large symmetry-breaking must be introduced before solitons can move.One difference with discrete NLS-type systems though is that these fundamentally dissipative structures can be attractors, since any radiation that is shed will be damped.A particularly novel aspect of the present work is a geometric interpretation of the threshold for solitons to move in terms of the shape of the bifurcation diagram of stationary solitons.Namely, a branch of moving solitons will be born when two independent stationary solutions undergo a fold for the same parameter value.
Further work is necessary to understand the precise nature of the moving solitons.The evidence here suggests that these moving solitons are periodic solutions rather than structures that are stationary in a moving frame.It would be interesting to understand the origins of the unusual 'staircase hysteresis' apparent in Figs. 9 and 10.Could the lattice discreteness have something to do with these apparent sudden jumps in velocity?In truth it could be that these jumps are caused by fold bifurcations, in order to find which we would need to be able to continue unstable versions of these travelling localised periodic solutions.Perhaps dedicated continuation algorithms for solutions of this type, such as those used in [23], can be used to probe more details.
It is worth commenting on the physical applicability of the results contained in this paper.Optical localised structures have attracted the attention of physicists due to the possibility of their use in optical memory, where the solitons serve as bits of information [27].Another possible application of cavity solitons is in the manipulation of Bose-Einstein condensates of excitons in systems with pumps, where bright and dark cavity solitons can form [36,12].In either case, the positions of the solitons can be controlled by a holding beam, which gives the possibility to process the information coded by the soliton [13].Microstructuring of the cavity allows for better control and localisation of the cavity solitons serving as bits of information.In the case of strong localisation, the dynamics of such cavities with microstructure systems can be described by discrete models like that considered in this paper.Also we believe that the results on both moving discrete cavity solitons reported here can be of practical importance in understanding how such bits of information be manipulated with minimal or no control input.
Finally, let us point to another avenue of research that is worthy of exploration, namely the link between the results obtained here and related results for continuum models with (spatially periodic) microstructure.If one doesn't apply the strong localisation limit, then models of optical cavities with microstructure can also posed in terms of continuum systems with periodic potentials [18].For such systems, when viewed as a spatial dynamical system, the fundamental states corresponding to the homogeneous states studied here, are perioidc orbits rather than equilibria.Hence, upon taking a Poincaré map approach, the arguments presented here about the break-up of snakes upon breaking the symmetry are also likely to apply.It would therefore be interesting to study in detail whether similar phenomena to those reported here are found in model systems of pumped optical cavities with periodic potentials.More generally, issues of symmetry-breaking due to inhomogeneities causing localised states to propogate has been observed in other physical systems, see for example [16,32].Explanations for such behaviour have been proposed by arguing the existence of certain global bifurcations in the spatial dynamics [32,7].Conversely, in fundamentally dissipative systems, it is known that inhomogenities can cause the failure of localised pulses to propogate [15] via a delayed pitchfork bifurcation.Clearly, there is scope for a more general understanding of the bifurcation events that can cause the onset of motion of localised structures in both discrete and inhomogeneous dissipative systems.

Figure 1 .
Figure 1.(a) Bifurcation diagram for stationary bright soliton solutions to (1) with parameter values (4) in the conservative case q = 0, Im(δ) = 0. Solid lines represent positions of the soliton branches that are temporally stable, whereas dashed lines represent unstable solitons.Thiner lines show the bifurcation diagram for homogeneous states, with dashed lines being stable and solid lines unstable.(b),(c),(d) Numerical simulation results in which the unstable solitons at the three labelled points in panel (a) were used as initial conditions.

Figure 2 .
Figure 2. (a) Similar to Fig. 1(a) but for q = 0.025.(b) Numerical simulation results in which the unstable solitons at the point labelled in panel (a) was used as initial conditions.(c) Similar to panel (a) but for q = 0 and Im(δ) = 0.05.(d) Similar but for q = 0.025, Im(δ) = 0.05.

Figure 3 .
Figure 3. (a) The transformation of the right-hand disconnected portion of the bifurcation diagram in Fig. 2(d) with increase in q.(b) The transformation of the same disconnected portion righthand disconnected portion of the bifurcation diagram with increase in Im(δ).

Figure 4 .
Figure 4. (a) Similar to Fig. 1(a) but for Im(δ) = 1.(b) The same results plot using the 'mass' norm M= n (|A n − A ∞ | 2 ).(c)Blow up of the small region of panel (b) in which the former thick solid and dashed curves are now depicted as thin black lines, with the addition of the asymmetric states that are spawned in pitchfork bifurcations from these on-site and off-site symmetric solitons.The solid and dashed lines in this panel represent the corresponding curves of solitons that occur after the inclusion of a small phase gradient q = 0.01.(d) Graphs of the stable asymmetric solitons at the points indicated in panel (c).

Figure 5 .
Figure 5. (a),(b),(c) Depicting the intersections between the stable and unstable manifolds of the fixed point of a two-dimensional reversible map near the fold point of a homoclinic snaking curve.In all three figures the unstable manifold is represented by the upper solid curve, the stable manifold by the lower dotted curve, and the fixed point set of the reversibility by a dashed horizontal line.Symmetric solitons occur when the stable and unstable manifolds intersect on the fixed point set, asymmetric solitons are represented by intersections away from this line.(d) The corresponding bifurcation diagram of homoclinic orbits, in which a solid line is used to represent a soliton that is temporally stable.The numbers next to each curve represent the number of unstable eigenvalues of the corresponding soliton.

Figure 8 .
Figure 8. Moving solitons obtained by direct numerical integration of (1) for parameter values (4) and Im(δ) = 0.4 and (a),(b) q = 0.025 and (c),(d),(e) q = 0.05.Panels (a) and (c) show stationary (gray) and moving (black) solitons as the maximum intensity against P .Panels (b) and (d) show the moving solitons as wave speed v against P .Panel (e) shows individual profiles of the moving solitons at specific instances of time.

inFigure 9 .
Figure 9. (a) Similar to Fig. 8(a) but for Im(δ) = 0.1 and q = 0.1.(b) Simulation results for the moving soliton starting from the unstable stationary soliton indicated in panel (a).The upper plot shows the intensity as a function of time and space.The middle plot shows the same information with the colour map plot on a logarithmic scale.The lower plot shows specific solution profiles at individual instants of time.(c) The same information as in (a) with the wave speed v plot against P .Arrows indicate jumps that occur under slow variation of parameters.

Figure 10 .
Figure 10.(a) A zoom of the indicated area in Fig. 9 near P = 1 upon decreasing P .Although not depicted here, a small region of hysteresis is observed close to each "jump" in v. (b),(c),(d) Field intensity profiles at isolated instances of time at the three parameter values indicated in panel (a).Note the different velocities of the travelling localised structures.