Generalized and multi-oscillation solitons in the nonlinear Schrödinger equation with quartic dispersion

Solitons of the GNLSE (1) are given by the $u 1$-component of homoclinic orbits to the origin $0=(0,0,0,0)$ of system (2), as a result of the traveling-wave ansatz.^8,10 Note that $0$ is an equilibrium in the zero-energy level $H≡0$ and, moreover, $0∈ Σ 1∩ Σ 2$⁠. Hence, when $0$ is a saddle, homoclinic orbits to it (approaching $0$ in both forward and backward time) may be $R 1$-symmetric, $R 2$-symmetric, $R ∗$-symmetric, or non-symmetric.⁸ In this Hamiltonian setting, homoclinic orbits lie in the same energy level as the equilibrium, and they are the limits of a family of periodic orbits (which are locally parametrized by the energy $H$⁠). Moreover, any homoclinic orbit is structurally stable; that is, it persists under changes of parameters—unless it undergoes a bifurcation. There has been extensive work on fourth-order, reversible, and Hamiltonian systems^{11–13,15–18} to understand the mechanism in which homoclinic solutions arise and disappear as parameters are varied.

System (2) is known^3,7,8 to have a pair of primary $R 1$-symmetric homoclinic orbits to $0$⁠, which are each other’s counterparts under the reversibility given by $R 2$⁠. The eigenvalues of the equilibrium $0$ become complex conjugate while the pair of homoclinic orbits exists, and this is known as a Belyakov–Devaney (BD) bifurcation.^11,12 Moreover, the primary homoclinic orbit disappears in a Hamiltonian Hopf (HH) bifurcation when $0$ ceases to be a saddle and becomes a center. Without loss of generality, we may consider the $( β 2,μ)$-plane for $β 4=−1$ and $γ=1$ to find the relevant bifurcations, which all occur along semi-parabolas.^8, Figure 1(a) shows how the curves BD and HH together form a single (analytically known) parabola that divides the $( β 2,μ)$-plane into regions I–III.^5,8 In regions I and II, the equilibrium $0$ is a real saddle and a saddle-focus, respectively. The primary $R 1$-symmetric homoclinic orbit exists in both these regions, and it is illustrated in panels (c) and (d) with temporal traces of $u 1$ that represent the corresponding soliton of the GNLSE (1). For simplicity, we will refer to such $u 1$-traces representing solitons as homoclinic solutions (as opposed to homoclinic orbits). Note that the soliton has non-oscillatory exponentially decaying tails in region I and oscillating decaying tails in region II.^3,7,8 In region III, the equilibrium $0$ is a center with purely imaginary eigenvalues and, hence, can no longer support homoclinic or other connecting orbits.

The Belyakov–Devaney bifurcation^11,12 implies that infinitely many homoclinic solutions of different symmetry types exist to the right of the curve BD in Fig. 1(a).^8,18 Figure 1(e) shows the primary $R 2$-symmetric homoclinic solution, which is special because it emerges from BD and exists throughout region II; that is, it disappears only at the Hamiltonian Hopf bifurcation curve HH. All other homoclinic orbits arise from BD as infinite families of different symmetry types, and they are associated with heteroclinic connections between the equilibrium $0$ and different periodic orbits, which we refer to as EtoP connections. The EtoP connections also persist under parameter variation; they come in pairs that disappear (as $β 2$ is increased) along curves, again parabolas, of fold bifurcations in region II, some of which are shown in Fig. 1(a). As we will discuss in more detail in Sec. III, the associated homoclinic orbits of a given family also disappear in pairs at fold curves that are close to the fold curve of the corresponding EtoP connection;⁸ the enlargement in Fig. 1(b) illustrates the ordering of these fold curves. We refer to the overall bifurcation structure as BD-truncated homoclinic snaking (as opposed to regular homoclinic snaking^19,20) and note that this phenomenon is also observed in a certain parameter regime of the Lugiato–Lefever equation.¹⁴ 

These earlier results form the starting point for the work presented in this paper. As the name suggests, the phenomenon of BD-truncated homoclinic snaking in the GNLSE (1) was investigated previously⁸ only in terms of homoclinic connections that come close to a single saddle periodic orbit, which are organized by pairs of EtoP connections with different symmetry properties. The key new insight we establish here is the existence of heteroclinic connections between saddle periodic orbits in the same energy surface, which are referred to as PtoP connections. These connecting orbits also persist with parameters and allow us to construct a plethora of new solitons of the GNLSE (1). These come in two distinct types.

1. PtoP connections are organizing centers that “generate” associated infinite families of homoclinic orbits to saddle periodic orbits with different symmetry properties. These solutions are not confined to the zero-energy surface, and they correspond to solitons with non-decaying oscillating tails of the GNLSE (1). This type of solution to a “periodic background” is referred to as a generalized soliton in the context of optics.²¹ 

2. PtoP connections can be found in the zero-energy level, which allows us to combine them with EtoP connections to obtain heteroclinic cycles that, in turn, give rise to additional infinite families of solitons with decaying tails. These solitons are characterized by the corresponding homoclinic orbits “visiting” not just one but several periodic orbits: they follow an EtoP connection from $0$ to a periodic orbit, perform a certain number of oscillations near it, then follow a PtoP connection to another periodic orbit, perform a certain number of oscillations near it, possible follow another PtoP connection, and so on, until they finally return to $0$ via another EtoP connection. In contrast to the solitons considered before,⁸ which visit a single periodic orbit, we refer to these solitons organized by PtoP connections as multi-oscillation solitons.

The overall picture is that there are infinitely more infinite families of solitons that are an integral part of the phenomenon of BD-truncated homoclinic snaking in system (2), in addition to those known from previous studies of (BD-truncated) homoclinic snaking.^8,14,19,20 All multi-oscillation solitons, which are homoclinic solutions to $0$ after all, are also created at BD, may have different symmetry properties, and disappear in fold bifurcations before HH is reached. Moreover, as we will show, the PtoP connections in the zero-energy surface also arise in pairs from the Belyakov–Devaney bifurcation BD and disappear in fold bifurcations before HH is reached. In turns, each such pair of PtoP connections gives rise to infinite families of homoclinic orbits of different symmetry properties to (each of) the two periodic orbits involved. In other words, we also find BD-truncated homoclinic snaking of homoclinic orbits to periodic orbits that are organized by the corresponding PtoP connections. This phenomenon for homoclinic orbits to periodic orbits, which has not been reported previously in the literature, effectively generalizes BD-truncated homoclinic snaking of homoclinic orbits to equilibria.^8,14 Since PtoP connections are not confined to the zero-energy level, they form surfaces in phase space. We compute and present such surfaces, and this provides a connection to theoretical results on the existence of surfaces of homoclinic solutions to periodic solutions in reversible systems.²² 

Clearly, PtoP connections require and go hand-in-hand with periodic orbits that they connect. Periodic orbits also form surfaces in phase space that are locally parametrized by the Hamiltonian function of the ODE.^11,12 A central role in the overall organization of connecting obits is played by the families of the periodic orbits that have as limits the pair of primary homoclinic orbits. Initially, for $β 2$ in region II and not too far from the curve BD in Fig. 1(a), there are three disjoint basic surfaces of periodic solutions. We compute these surfaces and show that they each intersect the zero-energy level infinitely many times. The overall picture of the soliton structure of the GNLSE (1) that, therefore, emerges is that the pair of primary homoclinic orbits generates infinitely many saddle periodic orbits in the zero-energy level. Sufficiently close to BD in region II, each such periodic orbit generates a pair of EtoP connections of $0$⁠. Moreover, any two periodic orbits in the zero-energy level can be connected by a pair of PtoP connections, which generate infinite families of homoclinic orbits to either of the periodic orbits. Concatenating any pair of EtoP connections with any number of PtoP connections creates infinitely many heteroclinic cycles from and to $0$⁠, each of which generates an entire menagerie of additional and arbitrarily complicated homoclinic orbits to the origin—and, hence, novel types of solitons of the GNLSE. In turn, each homoclinic orbit is the limit of families of periodic orbits, infinitely many of which lie in the zero-energy level so that they also generate pairs of EtoP and PtoP solutions and so on.

The work presented here makes heavy use of state-of-the-art computational methods that allow us to find different periodic, homoclinic, and heteroclinic orbits; to continue them in parameters; and to detect their bifurcations. In particular, we formulate connecting orbits between saddle objects as solutions of suitably defined two-point boundary value problems (2PBVPs).^8,23 In an approach called Lin’s method,²⁴ this involves finding intersections of computed parts of stable and unstable manifolds of equilibria and periodic orbits in a three-dimensional section. Since system (2) is Hamiltonian, it is necessary to ensure that solutions to the respective 2PBVPs are isolated, which is achieved by introducing the gradient of the conserved quantity $H$ multiplied with a perturbation parameter^8,25 to system (2). All these computations are performed with the continuation package Auto-07p²⁶ and its extension HomCont.²⁷ 

The paper is structured as follows. In Sec. II, we introduce the three basic surfaces $S 1 ±$ and $S ∗$ of periodic orbits that accumulate on the primary homoclinic orbit and its $R 2$-counterpart. In Sec. III, we then briefly summarize previous results⁸ regarding how EtoP connections to specific periodic orbits in these surfaces organize associated families of homoclinic orbits to $0$⁠; here, we also show a one-parameter bifurcation diagram that clarifies over which $β 2$-ranges they exist. This sets the stage for the discussion in Sec. IV of PtoP connections and how they generate generalized solitons. We first consider PtoP connections in the zero-energy, where they exist, and how they organize corresponding families of homoclinic orbits with different symmetry properties. We also calculate and present the surfaces of these PtoP connections, which are not restricted to have zero energy. We then provide in Sec. V a schematic diagram that shows how EtoP and PtoP connections can be concatenated into heteroclinic cycles. Specifically, we show how families of homoclinic orbits to $0$ with different symmetry properties are generated by combining different EtoP and PtoP connections. These additional types of homoclinic orbits are indeed multi-oscillation solitons of the GNLSE, and one-parameter bifurcation diagrams show their $β 2$-ranges of existence. We conclude in Sec. VI with a brief summary and an outlook on open questions and future research.

The primary homoclinic orbit from Fig. 1 is $R 1$-symmetric but not $R 2$-symmetric, meaning that it has an $R 2$-counterpart. As is illustrated in Fig. 2, for $β 2=0.4$⁠, this pair of primary homoclinic solutions gives rise to three basic surfaces; here and throughout, we fix, without loss of generality, $β 4=−1$ and $γ=1$⁠, and set $μ=1$ since all bifurcation curves in the $( β 2,μ)$-plane are parabolas; note here that Eq. (1) can be rescaled appropriately for different signs of $β 2$ and $β 4$⁠, as was shown previously.^8, Figure 2(a1) shows the pair of primary homoclinic orbits in terms of their $u 1$-profile, while panel (a2) shows them in the $( u 1, u 2)$-plane. They each come with a surface of periodic orbits that accumulates on the respective homoclinic orbit in the zero-energy level. These two surfaces $S 1 −$ and $S 1 +$ are shown in projection onto $( u 1, u 2,H)$-space in panels (b) and (c), respectively. The parametrization by the energy $H$ is explicit, and this allows us to easily identify the periodic orbits in the zero-energy level, which contains the origin $0$ and the pair of primary homoclinic orbits, which are also shown in Figs. 2(b) and 2(c). Note that $S 1 −$ and $S 1 +$ are both $R 1$-symmetric and $R 2$-counterparts of one another.⁸ Moreover, they each extend over only a bounded range of $H$⁠, where the largest value of $H$ occurs at the equilibria $E −$ and $E +$⁠. These are centers at $E ±= ( ± μ γ , 0 , 0 , 0 )$ (so only exists if $μγ>0$⁠) and each others $R 2$-counterparts with $H( E ±)=−6 μ 2/(γ β 4)$⁠. The union of the pair of primary homoclinic orbits is the limit of a third, $R ∗$-symmetric surface $S ∗$⁠, which is also shown in $( u 1, u 2,H)$-space in Fig. 2(d) together with the pair of primary homoclinic orbits. This surface has a maximal value of $H$ but extends to arbitrarily large negative values of $H$⁠; that is, it is not bounded from below.⁸ 

We refer to the surfaces $S 1 ±$ and $S ∗$ as the surfaces of basic periodic orbits, and they are shown together in Fig. 3 in projections onto $( u 1, u 2,H)$-space, namely, in a cutaway view that shows their parts for $u 2≤0$ up to the section $Σ={( u 1,0, u 3, u 4)}$ defined by $u 2=0$⁠. Notice that the respective other half, for positive $u 2$ on the other side of $Σ$⁠, is obtained from that shown by applying the transformation $R 1$⁠. In Fig. 3, this corresponds to a reflection in the $Σ$-plane; compare with Figs. 2(b)–2(d).

This cutaway representation allows us to compute and show the intersection curves of the surfaces $S 1 ±$ and $S ∗$ with the plane $Σ$⁠, which is very useful for understanding the geometry of these surfaces. The gray line in $Σ$ indicates the zero-energy level with $u 2=0$⁠, which contains the basic homoclinic orbits. As Fig. 3 illustrates, the basic surfaces $S 1 ±$ and $S ∗$ each spiral as they limit on the (union of) the two basic homoclinic orbits, which means that they intersect the zero-energy level infinitely many times. To be more specific, starting from the equilibria $E ±$ and decreasing the energy $H$⁠, one encounters in the zero-energy level the $R 1$-symmetric and $R 2$-related pair of periodic orbits $Γ 1 −$ on $S 1 −$ and $Γ 1 +$ on $S 1 +$⁠. Similarly, starting from large negative $H$ and increasing the energy, one encounters in the zero-energy level the $R ∗$-symmetric periodic orbit $Γ ∗$ on $S ∗$⁠. The three periodic orbits $Γ 1 ±$ and $Γ ∗$ are shown in full in Fig. 3. As can be seen, they are the “first” periodic orbits one encounters in the zero-energy level as the respective surface accumulates on the pair of homoclinic orbits. For this reason, we refer to $Γ 1 ±$ and $Γ ∗$ as the primary periodic orbits; note, in particular, that each of them is a single-loop orbit with exactly two intersection points with the section $Σ$⁠.

In this work, we demonstrate with the specific examples of the primary periodic orbits $Γ 1 ±$ and $Γ ∗$⁠, and for $β 2=0.4$⁠, how connecting orbits between them give rise to PtoP connections, which, in turn, generate generalized and multi-oscillation solitons. Indeed, the same is true for the further periodic orbits of $S 1 ±$ and $S ∗$ with $H=0$⁠, of which there are countably infinitely many. In contrast to the primary periodic orbits, further periodic orbits of $S 1 ±$ and $S ∗$ in the zero-energy level are multi-loop orbits with an increasing number of additional intersection points with the (three-dimensional) section $Σ$⁠. How they arise from the associated structure of $S 1 ±$ and $S ∗$⁠, and how these surfaces depend on the parameter $β 2$ are interesting subjects beyond the scope of this paper that will be discussed elsewhere.²⁸ 

As a starting point for our introduction of PtoP orbits, we briefly recall and discuss here what families of solitons of different symmetry types are generated in region II by EtoP connections from $0$ to the primary periodic orbits $Γ 1 ±$ and $Γ ∗$⁠,⁸ again for $β 2=0.4$ (and $μ=1$⁠). The dependence of the EtoP connections and associated solitons on $β 2$ is subsequently presented as a one-parameter bifurcation diagram, which introduces the different fold curves already shown in Fig. 1.

Figures 4 and 5 show how the EtoP connections to $Γ ∗$ and to $Γ 1 ±$⁠, together with their corresponding symmetric counterparts, generate heteroclinic cycles that organize homoclinic orbits to $0$ under mild conditions.^8,10,29,30 These homoclinic orbits closely follow an EtoP connection and make a number of loops around the respective periodic orbit before converging back to $0$ along a return EtoP connection. Homoclinic orbits with any number of loops can be found, and they form infinite families with increasing numbers of “humps” of the corresponding solitons. Since system (2) possesses two reversible symmetries, families of homoclinic orbits with different symmetry properties can be obtained through this mechanism.

Figure 4 shows the EtoP connections to the primary periodic orbit $Γ ∗$ and associated homoclinic solutions. Panels (a) and (b) show two EtoP connections from $0$ to $Γ ∗$ that exist simultaneously for $β 2=0.4$⁠. Here, panels (a1) and (a2) show the EtoP connections in $( u 1, u 2,H)$-space together with a cutaway view of the surface $S ∗$⁠; and panels (b1) and (b2) show their $u 1$-traces, which illustrates that the EtoP connections converge backward in time to $0$ and forward in time to $Γ 1 +$⁠. Note that these two EtoP connections are distinct in the sense that they are not related by any of the symmetries of (2). However, since $Γ ∗$ is $R ∗$-symmetric, the $R 1$-counterparts as well as the $R 2$-counterparts of these two EtoP connections exist as well. Note, in particular, that the $R 2$-counterparts provide return connections from $Γ ∗$ to $0$⁠. Hence, the EtoP connections in Figs. 4(a) and 4(b) together with their corresponding $R 1$- and $R 2$-counterparts generate $R 1$-symmetric and $R 1$-symmetry-broken, as well as $R 2$-symmetric and $R 2$-symmetry-broken homoclinic orbits, respectively.⁸ Some representative examples are shown in Figs. 4(c)–4(f) in terms of their $u 1$-traces, where we distinguished by color the two parts that follow a particular EtoP connection (or its symmetric counterpart). More specifically, panels (c1) and (c2) show solutions with $R 1$-symmetry that are generated from each of the two EtoP connections and their $R 1$-counterparts; here, each EtoP connection is effectively followed to its second maximum. Similarly, an EtoP connection and then the $R 1$-counterpart of the other combine to the $R 1$-symmetry-broken homoclinic solutions shown in panels (d1) and (d2). In complete analogy, panels (e1) and (e2) show homoclinic orbits with $R 2$-symmetry that are generated from each of the two EtoP connections and their $R 2$-counterparts. Combining an EtoP connection with the $R 2$-counterpart of the other gives the $R 2$-symmetry-broken homoclinic orbits shown in panels (f1) and (f2). We remark that the images of these solitons under the reversibility transformations given by $R 1$ and $R 2$ also exist; in the representation of the $u 1$-trace in Figs. 4(c)–4(f), these transformations are geometrically: reflection in the $u 1$ axis, and rotation by $π$ around the origin of the $(t, u 2)$-plane, respectively.

In the same style, Fig. 5 shows the EtoP connections to $Γ 1 +$ and associated homoclinic solutions that exist simultaneously. The difference is now that $Γ 1 +$ is not $R 2$-symmetric, and this means that the two distinct EtoP connections to $Γ 1 +$⁠, which are shown in panels (a1) and (a2) with a cutaway view of $S 1 +$⁠, can only be combined with their $R 1$-counterparts to generate heteroclinic cycles between $0$ and $Γ 1 +$⁠. (Of course, their images under the reversibility $R 2$ also exist, but they are heteroclinic cycles between $0$ and $Γ 1 −$⁠.) Panels (c1) and (c2) of Fig. 5 show examples of $R 1$-symmetric solitons that follow the respective EtoP connection to the first minimum near $Γ 1 +$⁠, and then follow the $R 1$-counterpart back to $0$⁠. In contrast, the solitons shown in panels (d1) and (d2) are composed by each of the EtoP connections in combination with the $R 1$-counterpart or the other EtoP connection, and this yields $R 1$-symmetry-broken solitons.

Because system (2) is reversible and Hamiltonian, all connecting orbits, the EtoP connections and associated homoclinic orbits, persist in open regions of parameter space.^13,31 Hence, once found, each such object can be continued as a solution branch in a chosen parameter, which we take to be the strength $β 2$ of the quadratic dispersion. In this way, one can compute a one-parameter bifurcation diagram that summarizes over what ranges of $β 2$ which types of connections and associated solitons can be found.^8, Figure 6 shows this kind of bifurcation diagram for the EtoP connections and different symmetric and symmetry-broken homoclinic orbits from Figs. 4 and 5. Here, the shown $β 2$-range includes the entire region II bounded by the bifurcations BD and HH, which are indicated by vertical lines and lie at $β 2≈±0.8164$ (for $μ=1$⁠). The different branches are represented by their $L 2$-norm $| |. | | 2$⁠, which is computed as an integral over the respective computed (part of the) connecting orbit.²⁶ The $L 2$-norm is a good choice to represent and distinguish branches of solitons because it increases with the number of maxima of the corresponding homoclinic orbits.

As Fig. 6 shows, the two EtoP connections to $Γ ∗$ from Figs. 4(a) and 4(b) lie on branches that emerge from BD and meet at a fold point at $β 2≈0.5753$ to form a single smooth curve. The location of this fold, which we refer to as EtoP $∗$, is highlighted in Fig. 6 by a vertical dashed line. The primary $R 1$-symmetric homoclinic orbit exists throughout regions I and II, forming a single branch that crosses BD and ends at HH. Note from Fig. 1(d) that, in region II, it can be interpreted as being part of the family of $R 1$-symmetric homoclinic orbits that are close to an EtoP cycle between $0$ and $Γ ∗$⁠. The other homoclinic orbits of this family, with more and more maxima as in Fig. 4(c), lie on branches that start at BD and meet in pairs at fold points to form smooth curves. In fact, these folds are points of symmetry breaking: they are also fold points where corresponding branches of $R 1$-symmetry-broken homoclinic orbits, as those in Fig. 4(d), meet to form single curves from and back to BD. Four curves of each family are shown in Fig. 6; note that two simultaneously existing symmetry-broken solutions have the same $L 2$-norm, which is why their two branches cannot be distinguished in Fig. 6. The situation for the families of $R 2$-symmetric and $R 2$-symmetry-broken homoclinic orbits is analogous. The primary $R 2$-symmetric homoclinic orbit from Fig. 1(e) is created at BD and ends at HH, and it can be interpreted as being part of the family of $R 2$-symmetric homoclinic orbits that follow an EtoP cycle between $0$ and $Γ ∗$⁠. As Fig. 6 shows, all other $R 2$-symmetric homoclinic orbits, such as those in Fig. 4(e), come in pairs that form single curves, which meet the curves of the corresponding pairs of $R 2$-symmetry-broken homoclinic orbits, as those in Fig. 4(e), at their respective fold points. Again, four curves of these families are shown in Fig. 6, and they illustrate that the folds of the families of $R 1$-symmetric and of the $R 2$-symmetric homoclinic orbits converge to the fold EtoP $∗$, from the left and from the right, respectively; compare with Fig. 1(b).

Similarly, the two EtoP connections to $Γ 1 +$ from Figs. 5(a) and 5(b) (and also those to $Γ 1 −$ by $R 2$-symmetry) lie on a single branch that emerges from and returns to BD, with a fold point that lies at $β 2≈0.6756$⁠. We refer to this fold as EtoP $1$, and it is also highlighted by a vertical dashed line in Fig. 6. For this type of $R 1$-symmetric EtoP connection, there only exist $R 1$-symmetric and $R 1$-symmetry-broken homoclinic orbits that come in pairs, as those in Figs. 5(c) and 5(d). Their branches, of which four are shown in Fig. 6, meet in the same way at joint fold points. These fold points are very close to and converge to the fold EtoP $1$ from the right; see Fig. 1(b).

We now show that one can find pairs of heteroclinic connections between two given saddle periodic orbits in the same energy level. More specifically, we compute and present such PtoP connections between the saddle periodic orbits $Γ ∗$ and $Γ 1 +$ and show that they emerge from the Belyakov–Devaney bifurcation BD and exist over a considerable $β 2$-range in region II.

Figure 7 shows in projection onto $( u 1, u 1, u 3)$-space a pair of PtoP connections between $Γ ∗$ and $Γ 1 +$ for $β 2=0.4$⁠. Each of these two PtoP connections is clearly seen to converge backward in time to $Γ ∗$ and forward in time to $Γ 1 +$⁠, and they are not related by symmetry. While the PtoP connection in Fig. 7(a) transitions swiftly from $Γ ∗$ to $Γ 1 +$⁠, the one in panel (b) has almost a complete loop near and around a third saddle periodic orbit $A Γ 1 +$ that also lies on the surface $S 1 +$⁠. All objects shown in Fig. 7 lie in the zero-energy level.

As we discuss next, such pairs of PtoP connections give rise to families of homoclinic orbits, to either $Γ ∗$ or $Γ 1 ±$ in this case. In the GNLSE, they correspond to generalized solitons²¹ with non-decaying oscillating tails, which are not confined to the zero-energy level. In Sec. V, we then show how PtoP connections in the zero-energy level can be combined with EtoP connections to form heteroclinic cycles from and back to $0$ that generate families of multi-oscillation solitons.

Figure 8 shows families of homoclinic orbits associated with the PtoP connections between $Γ ∗$ and $Γ 1 +$⁠, specifically examples of solution profiles and the associated bifurcation diagram. The $u 1$-traces of the two PtoP connections from Fig. 7 are shown in Fig. 8(a), which illustrates further how they connect $Γ ∗$ and $Γ 1 +$⁠. Notice the intermediate minimum in $u 1$ in panel (a2), which corresponds to the loop around the additional orbit $A Γ 1 +$ in Fig. 7(b). Since both $Γ ∗$ and $Γ 1 +$ are $R 1$-symmetric, the $R 1$-counterparts of the PtoP connections in Fig. 8(a), obtained by reflection in the $u 1$-axis, also exist: they correspond to return connections from $Γ 1 +$ to $Γ ∗$⁠. Hence, the two connections shown in panels (a) and their corresponding $R 1$-counterparts create different heteroclinic cycles between $Γ ∗$ and $Γ 1 +$⁠, and these organize different types of homoclinic solutions to $Γ ∗$ and to $Γ 1 +$⁠, respectively.

Panels (b) and (c) of Fig. 8 show examples of homoclinic orbits to $Γ ∗$ that are generated by PtoP cycles consisting of the two PtoP connections in Fig. 8(a) and their $R 1$-counterparts. More specifically, the homoclinic solutions in panels (b1) and (c1) follow the PtoP connection in panel (a1) from $Γ ∗$ to near $Γ 1 +$ and then follow its $R 1$-counterpart back to $Γ ∗$⁠. Those in panels (b2) and (c2) follow, similarly, the PtoP connection in panel (a2) and its $R 1$-counterpart. Finally, panels (b3) and (c3) show non-symmetric homoclinic solutions to $Γ ∗$⁠, which follow the PtoP connection in panel (a1) and the $R 1$-counterpart of that in panel (a2). The difference between panels (b) and (c) is that these homoclinic solutions make one and two loops around $Γ 1 +$⁠, respectively. Indeed, by following the respective PtoP connections longer, homoclinic solutions to $Γ ∗$ with any number of loops around $Γ 1 +$ can be found.

Figure 8(d) shows the one-parameter bifurcation diagram in $β 2$ of the PtoP connections, as well as the $R 1$-symmetric and non-symmetric homoclinic solutions to $Γ ∗$⁠, represented by the $L 2$-norm $| | u 1 | | 2$ of a suitably truncated portion of the solutions. As in Fig. 6, dashed vertical lines indicate the bifurcations BD and HH that bound the $β 2$-range in region II where the equilibrium $0$ is a saddle-focus. As is the case for the EtoP connections, the two PtoP connections in panels (a1) and (a2) lie on branches that emerge from BD and meet at a fold point to form a single curve. This fold, which we refer to as PtoP, occurs at $β 2≈0.5511$ and is also indicated by a vertical dashed line in Fig. 8(d); the PtoP connections in panels (a1) and (a2) are from the upper and the lower branch of this curve, respectively, and they no longer exist to the right of the line PtoP. We remark that the PtoP connections, which actually have an infinite $L 2$-norm, are represented in panel (d) by the $L 2$-norm of their truncation with ten loops near each periodic solution $Γ ∗$ and $Γ 1 +$⁠.

Also shown in Fig. 8(d) are the four branches of $R 1$-symmetric and non-symmetric homoclinic solutions to $Γ ∗$ with up to four oscillations near $Γ 1 +$⁠. Since homoclinic connections to periodic orbits also have infinite $L 2$-norm, they are also represented here by the $L 2$-norm of a truncation, with four oscillations near $Γ ∗$ in this case. As was the case for branches of homoclinic orbits to $0$ in Sec. III, the different branches of increasing norm connect in pairs at fold bifurcations. More specifically, the $R 1$-symmetric homoclinic solutions to $Γ ∗$ lie on branches that meet at fold points to form smooth curves. Here, the upper and lower branches represent homoclinic solutions generated by the PtoP connection in panel (a1) and its $R 1$-counterpart, and that in panel (a2) and its $R 1$-counterpart, respectively. The example solutions in panels (b1), (c1) and (b2), (c2) are indicated in Fig. 8(d) by dots and crosses on the correspondingly colored branches. Similarly, the non-symmetric homoclinic solutions to $Γ ∗$ lie on branches that meet at fold points, which are also the fold points of the corresponding $R 1$-symmetric solutions with the same number of loops around $Γ 1 +$⁠. Hence, these folds are symmetry-breaking bifurcations of homoclinic orbits to $Γ ∗$⁠, which is why we refer to this family as $R 1$-symmetry-broken. Note that the pairs of branches of symmetry-broken solutions are not distinguished by the norm and appear to be on top of each other in Fig. 8(d); the example solutions in panels (b3) and (c3) are indicated by diamonds on the correspondingly colored branches. As the number of loops around $Γ 1 +$ increases, the folds of the respective curves accumulate on the fold of PtoP connections from the left; notice that all these fold points have $β 2$-values that are very close to one another: in fact, they agree up to four decimal places with the value $β 2≈0.5511$ of PtoP.

PtoP cycles organize homoclinic solutions to any periodic solutions involved in the cycle. In particular, and as Fig. 9 shows, the pair of PtoP connections from Figs. 7 and 8(a) also generates families of homoclinic solutions to the periodic solution $Γ 1 +$⁠. Example solutions with one and two oscillations near $Γ ∗$ are shown in Figs. 9(a) and 9(b), and panel (c) is the associated bifurcation diagram that also includes the curve of PtoP connections. Comparison with Figs. 8(b)–8(d) shows that the different families of $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions are organized in the same way. Specifically, Figs. 9(a1) and 9(b1) show $R 1$-symmetric homoclinic solutions that follow the $R 1$-counterpart of the PtoP connection in Fig. 8(a1) from $Γ 1 +$ to near $Γ ∗$ and then follow this PtoP connection back to $Γ ∗$⁠. Figure 9(a2) and 9(b2) show the same type of homoclinic orbit, but now for the other PtoP connection, and panels (a3) and (b3) show the corresponding $R 1$-symmetry-broken homoclinic solutions obtained by combining the two PtoP connections. The respective pairs of branches of homoclinic orbits with any number of loops around $Γ ∗$⁠, of which the first four are shown in Fig. 9(c), meet at fold points to form smooth curves; the locations of the different example solutions are again marked. The folds of the respective pairs coincide and are symmetry-breaking bifurcations of the homoclinic solutions to $Γ 1 +$⁠. As the number of loops around $Γ ∗$ increases, these fold points converge to the fold of the curve of PtoP connections; their $β 2$-values are indistinguishable (up to the accuracy of the continuation) from the value $β 2≈0.5511$ of the point PtoP.

Figures 8(d) and 9(c) show how the PtoP connections from Fig. 7 organize families of $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions to either $Γ ∗$ and $Γ 1 +$⁠. For increasing $β 2$⁠, the curves of all these homoclinic solutions have folds, effectively all at the point PtoP, and disappear before reaching the bifurcation HH. Continuing these branches of homoclinic solutions to periodic orbits for decreasing $β 2$ turns out to be numerically more challenging. They have been computed to considerably lower values than shown in Figs. 8(d) and 9(c), where their $L 2$-norm increases significantly for decreasing $β 2$⁠. Nevertheless, our numerical computations strongly suggest that the homoclinic solutions to $Γ ∗$ and $Γ 1 +$ are all also created at the Belyakov–Devaney bifurcation BD; however, this limit is not quite reached in the computation of the branches due to increased sensitivity of the continuations, which is arguably due to the oscillating nature of their tails.

Note that the $R 2$-symmetric counterparts of the homoclinic solutions in Figs. 8(b1)–8(c3) are also homoclinic solutions to $Γ ∗$⁠; however, they now make multiple loops around $Γ 1 −$⁠. Since they have the same $L 2$-norm, their curves in the one-parameter bifurcation diagram are identical to those shown in panel (d). The $R 2$-symmetric counterparts of the homoclinic solutions to $Γ 1 +$ in Figs. 9(a1)–9(b3) are simply the corresponding homoclinic solutions to $Γ 1 −$⁠. Again, these have the same $L 2$-norm and so occur along the curves already shown in the bifurcation diagram in panel (c); hence, each curve represents homoclinic solutions both to $Γ 1 +$ and to $Γ 1 −$⁠.

All these families of homoclinic orbits to periodic orbits of the ODE system (2) are generalized solitons (with non-decaying oscillating tails) of the GNLSE (1). In other words, as for the regular solitons with decaying tails,⁸ we find that the generalized solitons are also organized by the phenomenon of BD-truncated homoclinic snaking, where the homoclinic orbits are now to periodic orbits rather than to the equilibrium $0$⁠; compare Figs. 8(d) and 9(c) with Fig. 6. An important difference is that BD-truncated homoclinic snaking to periodic orbits exists not only in the zero-energy level, as we will show next.

Since $Γ ∗$ and $Γ 1 +$ lie in the zero-energy level, so do the PtoP connections in Fig. 7. However, the periodic orbits can be continued to obtain the surfaces $S ∗$ and $S 1 +$ from Figs. 2(c) and 2(d). As Fig. 10 shows, each of the two PtoP connections can similarly be continued in the energy $H$ to PtoP connections of corresponding periodic orbits in the respective energy level of $S ∗$ and $S 1 +$⁠. Here, panels (a) and (b) show the surfaces $S ∗$ and $S 1 +$ in $( u 1, u 2,H)$-space, together with a number of representative PtoP connections of the corresponding one-parameter family; these are not rendered as a surface for clarity of presentation, and the PtoP connections from Fig. 7 are highlighted. Note that one could also continue each homoclinic solution in Figs. 8 and 9 to find their respective continuations in $H$⁠. We rather continue here the two PtoP connections themselves because they organize the many homoclinic solutions to the periodic orbits in the surfaces $S ∗$ and $S 1 +$⁠. In this way, we present a global picture that allows one to make connections to theoretical results on the existence of surfaces of homoclinic solutions to periodic solutions in reversible systems.²² 

Figure 10 clearly shows that the PtoP cycles and associated homoclinic orbits can be found over a range of the energy $H$⁠; specifically, over the range of the joint existence of the continuations of the periodic orbits $Γ ∗$ and $Γ 1 +$⁠. As we continue the PtoP connections in the zero-energy level (black curves) toward positive and negative values of $H$⁠, they reach a maximum and a minimum energy level, respectively, which correspond to folds with respect to $H$ of periodic orbits on the surfaces $S ∗$ and $S 1 +$⁠. We remark that these continuations are quite delicate near such local extrema of the energy, and it is difficult to continue a branch past such a fold. Nevertheless, we are able to identify the respective PtoP connection past the fold and resume the continuation in this way. More specifically, as we continue the PtoP connection in Fig. 7(a) toward the positive values of $H$⁠, it reaches the energy level with $H≈3.3$⁠, which is the global maximum of $S ∗$⁠; see Fig. 10(a). When we continue it toward negative values of $H$⁠, it reaches the energy level with $H≈−6.2$⁠, which is the global minimum of $S 1 +$⁠. Likewise, when the connection in Fig. 7(b) is continued toward negative values of $H$⁠, it also reaches the energy level with $H≈−6.2$⁠; see Fig. 10(b). However, as we continue it toward positive values of $H$⁠, it only reaches the energy level with $H≈0.29$⁠. This is because this family of PtoP connections has an extra loop near the continuation of the $R 1$-symmetric periodic orbit $A Γ 1 +$ on $S 1 +$⁠, which only exists up to the local maximum of $S 1 +$ at energy level $H≈0.29$⁠. The close passage of the PtoP connection in Fig. 7(b) near $A Γ 1 +$ already demonstrates that it is possible to find PtoP connections between $Γ ∗$ and $A Γ 1 +$ and between $Γ 1 +$ and $A Γ 1 +$⁠. Overall, past any fold on the surface $S ∗$ or $S 1 +$ the respective PtoP connection connects different periodic orbits, and they all exist over a certain range of $H$⁠.

Any PtoP connections in a given energy level can be combined to create heteroclinic PtoP cycles between different periodic orbits and associated homoclinic orbits, that is, additional families of generalized solitons that exist over certain ranges of $H$ near $H=0$⁠. In particular, there are infinitely many families of homoclinic orbits to the periodic orbits $Γ ∗$ and $Γ 1 ±$ in the zero-energy level. This geometric picture allows us to conclude the following. As $Γ ∗$ and/or $Γ 1 ±$ approach the primary homoclinic orbits in the zero-energy level, the families of the associated homoclinic orbits to the periodic orbits $Γ ∗$ and/or $Γ 1 ±$ approach connecting orbits from and back to the primary homoclinic orbits to $0$⁠. Taking the limit may, hence, be a mechanism for creating homoclinic orbits to homoclinic orbits, which are also referred to as super-homoclinic orbits.³² These objects are limits of generalized solitons in the GNLSE, and it is an interesting question beyond the scope of this paper whether and how they act as organizing centers for the different types of solitons.

In the zero-energy level, one finds the PtoP connections between $Γ ∗$ and $Γ 1 +$ from Fig. 7, as well as their $R 1$- and $R 2$-counterparts that involve $Γ 1 −$⁠. Together, they are involved in creating more complicated heteroclinic PtoP cycles between the three periodic orbits $Γ ∗$⁠, $Γ 1 +$⁠, and $Γ 1 −$⁠. As we show now, these PtoP cycles can be combined with the EtoP connections from Sec. III effectively in any combination to create more complicated heteroclinic cycles from and back to $0$⁠. In turn, these organize a plethora of homoclinic orbits to $0$⁠, which all correspond to solitons of the GNLSE with decaying tails that feature episodes of oscillations near any of these periodic orbits. For this reason, we refer to them as multi-oscillation solitons.

This whole mechanism is summarized in Fig. 11, which shows a schematic diagram in the form of a graph of possible connections between $0$⁠, $Γ ∗$⁠, and $Γ ±$ at $β 2$-values where all of them exist, such as for $β=0.4$ as we used throughout. As we discussed in Secs. III and IV, there exist two distinct EtoP and PtoP connections for fixed parameter values, and they are represented in Fig. 11 by single directed edges that are numbered $①$ to $⑧$⁠; here, the arrows indicate the direction of time and reference to previous figures indicates connections we already encountered. The edges labeled $⑦$ and $⑧$ represent pairs of PtoP connections between the $R 2$-counterparts $Γ 1 +$ and $Γ 1 −$⁠, which are jointly represented by a single bi-colored circle. Furthermore, in this labeling, the even connections $②$⁠, $④$⁠, $⑥$ and $⑧$ correspond to the $R 1$- or $R 2$-counterparts of the odd connections $①$⁠, $③$⁠, $⑤$ and $⑦$⁠, respectively.

Any path in the directed graph in Fig. 11 that starts and ends at $0$ represents certain heteroclinic cycles with associated infinite families of homoclinic orbits to $0$ of the different symmetry types, which are all solitons of the GNLSE. Notice that the families of homoclinic solutions that were previously found⁸ and are summarized in Sec. III are organized by the cycles that are formed by only the connections $①$ and $②$ and by $③$ and $④$⁠. That is, they reach $Γ ∗$ or $Γ 1 +$ along the EtoP connection $①$ and $③$⁠, loop multiple type times close to the respective periodic orbit, and then converge back to $0$ along the EtoP connection $②$ and $④$⁠. Importantly, the existence of the PtoP connections $⑤$ and $⑧$ that connect the periodic orbits $Γ ∗$ and $Γ 1 ±$ allows for the construction of more complicated homoclinic solutions to $0$⁠. Each arrow represents two distinct connections from one saddle object to another, which makes this mechanism more complex than Fig. 11 suggests.

For each of these two types of more complex homoclinic solutions to $0$⁠, we present again representative examples as well as the bifurcation diagram in $β 2$ with the corresponding branches.

Figure 12 shows families of homoclinic solutions to $0$ for $β 2=0.4$ that are organized by the heteroclinic cycles with the edge sequence (⁠ $①$⁠, $⑤$⁠, $⑥$⁠, $②$⁠). This type is determined by the fact that the connections $⑥$ and $②$ are the $R 1$-counterparts of the connections $⑤$ and $①$⁠, respectively. As panels (a),(b) and (d),(f) of Fig. 12 show, multi-oscillation solitons of this type correspond to homoclinic orbits that leave $0$ to oscillate around $Γ ∗$⁠, transition to oscillate around $Γ 1 +$⁠, transition back and oscillate around $Γ ∗$ again, and finally return to $0$⁠. The connections can be made by either of the pair of EtoP from Fig. 4(a) and PtoP connections from Fig. 8(a), and the number of oscillations of each of these episodes differs from branch to branch, as is illustrated in the bifurcation diagram in Figs. 12(c) and 12(f).

More specifically, Fig. 12(a1) shows $R 1$-symmetric homoclinic solutions of this type: from $0$ it follows the EtoP connection in Fig. 5(a2); makes two loops close to $Γ ∗$⁠; follows the PtoP connection in Fig. 8(a1); makes one loop close to $Γ 1 +$⁠; follows the $R 1$-counterpart of the PtoP connection in Fig. 8(a1); makes two loops close to $Γ ∗$⁠; and finally follows the $R 1$-counterpart of the EtoP connection in Fig. 5(a2) back to $0$⁠. The homoclinic solution in panel (a3) of Fig. 12 differs only in that it follows the PtoP connection in Fig. 8(a2) and its $R 1$-counterpart. The homoclinic solution in panel (a2), on the other hand, is $R 1$-symmetry-broken, because it first follows the EtoP and PtoP connections in Figs. 5(a2) and 8(a1), respectively, and then converges back to $0$ by following the $R 1$-counterparts of the PtoP and EtoP connections in Figs. 8(a2) and 5(a2), respectively. Fixing the number of loops near $Γ ∗$ and increasing the number of loops near $Γ 1 +$ is a way to generate families of $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions, and panels (b1), (b2), and (b3) show the ones with one additional loop.

Figure 12(c) shows the bifurcation diagram in $β 2$ of these families, where solutions are again represented by their $L 2$-norm $| |. | | 2$⁠. Here, we show the curves of the homoclinic solutions with two to five loops near $Γ 1 +$⁠; the bifurcations BD and HH that bound region II are shown by dashed vertical lines, as are the fold points EtoP $∗$, EtoP $1$, and PtoP. For reference, all previously found branches of homoclinic solutions⁸ from Sec. III are shown in light gray. The two $R 1$-symmetric homoclinic solutions in panels (a1) and (a3) of Fig. 12, with the same number of loops near $Γ ∗$ and $Γ 1 +$⁠, respectively, lie on branches that connect at a fold point to form a smooth curve. As is the case for the simpler homoclinic orbits, this fold is also a point of symmetry breaking from which bifurcates a pair of branches of $R 1$-symmetry-broken homoclinic solutions; the one in panel (a2) lies on one of these branches (which coincide because they have the same $L 2$-norm). This picture of branches connecting as pairs at joint fold points is repeated in panel (c) when the number of loops near $Γ 1 +$ is increased; the positions of the homoclinic solutions in panels (b1)–(b3) are indicated on the respective branches.

More families of homoclinic solutions associated with type (⁠ $①$⁠, $⑤$⁠, $⑥$⁠, $②$⁠) can be generated by fixing the oscillations near $Γ 1 +$ and increasing the number of loops near $Γ ∗$⁠. As before, there exist $R 1$-symmetric and $R 1$-symmetry-broken families of homoclinic solutions, and examples of them are shown in panels (d) and (e) of Fig. 12. Notice that those in panels (d1)–(d3) are like the ones in panels (a1)–(a3), but with one further loop near $Γ ∗$⁠; similarly, the homoclinic solutions in panels (e1)–(e3) have two extra loops near $Γ ∗$⁠. As the bifurcation diagram in panel (f) shows, all of these $R 1$-symmetric homoclinic solutions also lie in pairs on single curves with two branches that meet at fold points, from which branches of the $R 1$-symmetry-broken homoclinic solutions bifurcate.

Figure 13 shows families of homoclinic solutions to $0$ that are organized by the heteroclinic cycles with the edge sequence (⁠ $③$⁠, $⑥$⁠, $⑤$⁠, $④$⁠)). Shown are examples of $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions in panels (a) and (b), and of $R 2$-symmetric and $R 2$-symmetry-broken homoclinic solutions in panels (d) and (e), respectively; the corresponding bifurcation diagrams in $β 2$ are shown in panels (c) and (f). For the $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions, the connections $⑤$ and $④$ are the $R 1$-counterparts of the connections $⑥$ and $③$⁠, respectively. All these homoclinic solutions loop close to $Γ 1 +$ before converging back to $0$⁠. For the $R 2$-symmetric and $R 2$-symmetry-broken homoclinic solutions, the connections $⑤$ and $④$ are the $R 2$-counterparts of the connections $⑥$ and $③$⁠, respectively; so now all these homoclinic solutions loop close to $Γ 1 −$ before converging back to $0$⁠.

To be more specific, the homoclinic solution in Fig. 13(a1) follows the EtoP connection in Fig. 5(a2) and makes two loops near $Γ 1 +$⁠. After that, it follows the $R 1$-counterpart of the PtoP connection in Fig. 8(a1) and makes one loop near $Γ ∗$⁠; subsequently, it follows the PtoP connection in Fig. 8(a1) and the $R 1$-counterpart of the EtoP connection in Fig. 5(a1) to converge back to $0$⁠. Similarly, the homoclinic solution in panel (a3) also makes use of the EtoP connection in Fig. 5(a1); however, it is now associated with the PtoP connection in Fig. 8(a2). The homoclinic solutions in panels (b1) and (b3) are similar to the ones in panels (a1) and (a3), but make one additional loop near $Γ ∗$⁠. The corresponding $R 1$-symmetry-broken homoclinic solutions are shown in panels (a2) and (b2), respectively. Similar to what we found in Fig. 12, one can also generate other families of homoclinic solutions by fixing the number of loops they make near $Γ ∗$ and by increasing the number of loops near $Γ 1 +$⁠. As the bifurcation diagram in Fig. 13(c) shows, all these $R 1$-symmetric and $R 1$-symmetry-broken homoclinic solutions also lie in pairs on single curves with two branches that meet at fold points.

It is also possible to generate families of $R 2$-symmetric and $R 2$-symmetry-broken homoclinic solutions of type (⁠ $③$⁠, $⑥$⁠, $⑤$⁠, $④$⁠), and examples are shown in panels (d) and (e) of Fig. 13. As panels (d1),(e1) and (d3),(e3) show, such $R 2$-symmetric homoclinic solutions initially loop twice close to $Γ 1 +$⁠, then follow the respective PtoP connections to loop close to $Γ ∗$ and, finally, loop three times close to $Γ 1 −$⁠, instead of $Γ 1 +$⁠, before returning to $0$⁠. As before, the corresponding $R 2$-symmetry-broken homoclinic solutions can be constructed by combining the two different PtoP connections, and they are shown in panels (d2) and (e2). The difference between the homoclinic orbits in panels (d) and (e) is that they have one and two loops near $Γ ∗$⁠, respectively. Taking more loops near $Γ ∗$ creates additional homoclinic orbits of the respective family, of which the first four are shown in the bifurcation diagram in Fig. 13(f). Also the $R 2$-symmetric homoclinic solutions lie in pairs on curves with two branches that meet at fold points, from which branches of the corresponding $R 2$-symmetry-broken homoclinic solutions emerge.

Figures 12 and 13 illustrate the existence of an entire menagerie of additional homoclinic orbits to $0$ that visit any of the periodic orbits $Γ ∗$ and $Γ 1 ±$⁠. These additional homoclinic solutions are indeed all multi-oscillation solitons of the GNLSE. More specifically, for any chosen edge sequence, one can find the associated families of symmetric and symmetry-broken homoclinic solutions; these include $R 2$-symmetric and $R 2$-symmetry-broken ones when the respective heteroclinic cycles are not $R 2$-symmetric itself. Figures 12 and 13 show certain families given by increasing only the number of loops near one of the constituent periodic orbits. However, we observe and conjecture that homoclinic solutions exist for any prescribed number of loops near any of the constituent periodic orbits of a given type. Note also that the $R 2$-counterparts of the shown homoclinic solutions also exist. We do not show them here because, on the level of the $u 1$-traces shown in these figures, they simply correspond to reflections of $u 1$ in the $t$-axis, so that maxima become minima and vice versa, and $Γ 1 +$ and $Γ 1 −$ are exchanged.

As the bifurcation diagrams in panels (c) and (f) of Figs. 12 and 13 illustrate, all these multi-oscillation solitons lie on curves organized by BD-truncated homoclinic snaking. Indeed, in spite of increasing numerical sensitivity due to their more complicated nature, we are able to continue all branches of multi-oscillation solitons for decreasing $β 2$ to quite near the point BD to confirm this. All these solutions require the existence of the respective EtoP and PtoP connections. Since the fold PtoP of the latter occurs for smaller $β 2$ than the folds EtoP $∗$, EtoP $1$ of the different EtoP connections, all these multi-oscillation solitons involving the periodic orbits $Γ ∗$ and $Γ 1 ±$ lie on branches that have folds near PtoP. In fact, as was also the case for the folds of generalized solitons in Sec. IV, the fold points of multi-oscillation solitons are indistinguishable in their $β 2$-value (as found by continuation) from that of the fold PtoP.

We showed that there exist infinitely many infinite families of solitons of the GNLSE, beyond the “regular” solitons that were found before.⁸ They arise due to the presence of robuts transitions given by PtoP connecting orbits between different saddle periodic orbits, and they come in two flavors: (1) generalized solitons with non-decaying oscillating tails and (2) multi-oscillation solitons that feature episodes of different types of oscillations, yet have decaying tails. We presented a systematic picture of how different PtoP connections can be combined with each other and with associated EtoP connections from/to the zero solution to generate new families of solitons with given symmetry properties.

These results were obtained with a dynamical system approach by translating solitons of the GNLSE into homoclinic orbits of the ODE system (2), which is Hamiltonian and features two reversible symmetries, $R 1$ and $R 2$⁠. In this framework, the two new types of solitons involve PtoP connections between saddle-type periodic orbits of the ODE. For the specific example of the primary periodic orbits $Γ ∗$ and $Γ 1 ±$ in the zero-energy level, we showed how PtoP connections between them give rise to infinite families of homoclinic orbits to either of the constituent periodic orbits. In terms of the GNLSE, these solutions are indeed generalized solitons with non-decaying oscillating tails; notably, these PtoP connections and associated generalized solitons are not confined to the zero-energy level. PtoP connections in the zero-energy level, on the other hand, can be combined with EtoP connections from $0$ to the respective periodic orbits to obtain heteroclinic cycles from and back to $0$⁠. These are the multi-oscillation solitons of the GNLSE, which feature different episodes of oscillations near any of the periodic orbits they visit. We stress that these solitons are not bound states or “molecules” of two or more interacting primary solitons.

All these additional solitons can be found and then continued in parameters. By continuation in the quadratic dispersion parameter $β 2$⁠, we showed that they are also directly associated with and organized by BD-truncated homoclinic snaking—even the generalized solitons that correspond to homoclinic solutions to periodic orbits. Hence, pairs of branches of any given family emerge from the Belyakov–Daveney bifurcation BD and meet at fold bifurcations, which are also symmetry-breaking bifurcations. The folds of a given family accumulate on the fold of constituent EtoP or PtoP connections with the lowest value of $β 2$⁠. Note, in particular, that all these additional solitons of the GNLSE exist for $β 2=0$⁠, that is, for the case of pure quartic dispersion.

We presented a systematic picture in the form a directed graph with numbered edges that encodes all possible heteroclinic cycles involving the equilibrium $0$ and the primary periodic orbits $Γ ∗$ and $Γ 1 ±$⁠. Any edge sequence in this graph from and back to $0$ generates associated families of multi-oscillation solitons of the respective type. While we supported it only for certain families, it is a natural conjecture that this statement holds in this generality. What is more, the graph can also be used to construct more complicated generalized solitons by selecting edge sequences of this graph that avoid the equilibrium. However, the story is even more complicated: already the basic surfaces $S ∗$ and $S 1 ±$ intersect the zero-energy level infinitely many times as they converge onto the (union of the) primary homoclinic orbits. Hence, the graph we presented can be extended with infinitely more periodic orbits of $S ∗$ and $S 1 ±$ and associated edges that represent EtoP connections from $0$ to them, as well as PtoP connections between all these additonal periodic orbits. We conjecture that any edge sequence in this much extended graph also generates families of multi-oscillation solitons of the corresponding type, as was shown for the cases presented here. But there is even more: any of the thus created solitons is itself the limit of a family of periodic orbits! We suspect that each of them, when continued in the energy, creates infinitely many periodic orbits in the zero-energy level that also feature EtoP and PtoP connections, and so on. The picture that emerges is truly a never-ending and “self-similar” plethora of additional solitons, which we conjecture are all created at the point BD as part of the overall phenomenon of BD-truncated homoclinic snaking of the ODE system (2) and, hence, the GNLSE (1). More generally, we conjecture that a Belyakov–Devaney bifurcation in any fourth-order Hamiltonian system with two reversible symmetries is an organizing center near which all of the solutions reported here can be found in the same way.

There are certainly avenues for future research that arise from our study. First, we already alluded to the geometry of the basic surfaces $S ∗$ and $S 1 ±$⁠, which we only considered here for the fixed parameter value $β 2=0.4$⁠. As $β 2$ is increased, these surfaces interact and their different parts connect differently in a sequence of transitions through specific types of symmetry-breaking bifurcations, period- $k$ multiplying bifurcations, and saddle-node bifurcations; these results will be presented elsewhere.²⁸ Second, it remains a considerable challenge to formally prove the existence of the EtoP and PtoP connections and all the different families of homoclinic orbits they generate, presumably starting with the specific ones we already found by numerical continuation. This could possibly be done by setting up Lin’s method as a theoretical tool (rather than a numerical one).^10,33 Moreover, our computations based on boundary-value problem formulations come with error bounds; hence, the existence of any connecting orbit we showed here can be proved, in principle, in a computer-assisted way.^34,35 Another issue is the stability of the different solitons as solutions of the PDE. Only the primary single-hump soliton is stable,^3,5 while all the other solitons we discovered appear to be unstable (but some are only weakly unstable⁸). Establishing this observation rigorously remains a considerable challenge for future work and will require state-of-the art techniques, such as the computation of Evans functions.³⁶ 

Localized structures in other spatially extended systems may be studied and explored, in the same spirit, with the numerical continuation approach employed here. In particular, recent experiments³⁷ have shown the feasibilities of creating waveguides with higher even-order dispersions, which are described by further generalizations of the GNLSE with corresponding non-zero dispersion terms of order six, eight, ten, or even higher. Other examples of interesting PDEs in this context are a widely studied version of the NLSE describing solitary wave formation in inhomogeneous media,^38–41 and the Lugiato–Lefever equation,¹⁴ which both are known to feature homoclinic snaking. Finally, we mention that an interplay between nonlinearity and quartic dispersion may occur in a wider class of system—beyond single optical waveguides. For example, a three-dimensional Gross–Pitaevskii solitary wave in a Bose–Einstein condensate may become stable with purely quadratic dispersion, realized via a shaken optical lattice.⁴² Moreover, spatially extended optical and photonic systems may feature intrinsic localized modes (ILMs),⁴³ which are in a sense “spatial analogs” of solitons; a specific example is the emergence of ILM in an optical planar waveguide array with sufficiently strong Kerr nonlinearity and quadratic dispersion. A future extension might be to consider spatially coupled waveguides with a quartic component of the dispersion, or dominantly quartic waveguides.

A. Giraldo was supported by the Korea Institute for Advanced Study under KIAS Individual Grant No. CG086101.

The authors have no conflicts to disclose.

Ravindra Bandara: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrus Giraldo: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Neil G. R. Broderick: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Bernd Krauskopf: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.