Generalized and multi-oscillation solitons in the
Nonlinear Schrödinger Equation with quartic dispersion

Optical solitons exist in photonic waveguides, notably, in optical fibers, due to the balance between dispersion and the Kerr nonlinearity, and they have been the focus of many recent experimental [35, 9] and theoretical [37, 38, 39, 25, 1, 6] studies. In the first instance, the dispersion has been modeled with a quadratic term. However, it has recently been discovered that waveguides with only quartic dispersion may support solitons with oscillating tails, which are much narrower in a well defined way (in terms of how the amplitude relates to the width [9]). This discovery of quartic solitons has led to considerable interest in characterizing the interplay between quartic and quadratic dispersion and the types of solitons this may entail [2, 6, 37, 38, 39]. To set the stage, pulse propagation along a waveguide with quadratic and quartic dispersion is governed by the Generalized Nonlinear Schrödinger Equation (GNLSE)

$$\frac{\partial A}{\partial z}=i\gamma|A|^{2}A-i\frac{\beta_{2}}{2}\frac{\partial^{2}A}{\partial t^{2}}+i\frac{\beta_{4}}{24}\frac{\partial^{4}A}{\partial t^{4}}$$

(1)

for the complex pulse envelope $A(z,t)$; here $\beta_{2}$ and $\beta_{4}$ are the respective dispersion coefficients, and $\gamma$ represents the strength of the Kerr nonlinearity, which is cubic.

Theoretical and numerical work on the partial differential equation (PDE) (1) has shown the persistence of a single-hump (pulse) soliton [39], and the existence of infinitely many solitons with different symmetry properties and numbers of humps [6, 33] for $\beta_{4}<0$ and a range of positive and negative values of $\beta_{2}$. These results are obtained by considering a traveling-wave ansatz, which transforms the PDE into a Hamiltonian and reversible system of ordinary differential equations (ODEs) [16, 17, 11]. For the GNLSE  (1) with $A(z,t)=u(t)e^{i\mu z}$, one thus obtains the system of ODEs [6]

$$\cfrac{d\mathbf{u}}{dt}=\left({\begin{array}[]{cc}u_{2}\\ u_{3}\\ u_{4}\\ \cfrac{24}{\beta_{4}}\left(\cfrac{\beta_{2}}{2}u_{3}+\mu u_{1}-\gamma u_{1}^{3}\right)\\ \end{array}}\right),$$

(2)

where $\mathbf{u}=(u_{1},u_{2},u_{3},u_{4})=\left(u,\frac{du}{dt},\frac{d^{2}u}{dt^{2}},\frac{d^{3}u}{dt^{3}}\right)$. System (2) is Hamiltonian (thus, also volume-preserving) with the energy function [6, 33]

$$H(\mathbf{u})=u_{2}u_{4}-\frac{1}{2}u_{3}^{2}-\left(\frac{6\beta_{2}u_{2}^{2}-6\gamma u_{1}^{4}+12\mu u_{1}^{2}}{\beta_{4}}\right)$$

as conserved quantity. Moreover, the system of ODEs (2) has two reversibilities given by the transformations

	$\displaystyle R_{1}:$	$\displaystyle(u_{1},u_{2},u_{3},u_{4})$	$\displaystyle\rightarrow(u_{1},-u_{2},u_{3},-u_{4}),$		(3)
	$\displaystyle R_{2}:$	$\displaystyle(u_{1},u_{2},u_{3},u_{4})$	$\displaystyle\rightarrow(-u_{1},u_{2},-u_{3},u_{4})$		(4)

with corresponding reversibility sections

	$\displaystyle\Sigma_{1}$	$\displaystyle=$	$\displaystyle\{\mathbf{u}\in\mathbb{R}^{4}:u_{2}=u_{4}=0\}\quad\mathrm{and}$		(5)
	$\displaystyle\Sigma_{2}$	$\displaystyle=$	$\displaystyle\{\mathbf{u}\in\mathbb{R}^{4}:u_{1}=u_{3}=0\}$		(6)

that are (pointwise) invariant under $R_{1}$ and $R_{2}$, respectively [16, 11, 34, 14]. Hence, if $\mathbf{u}(t)$ is a trajectory or solution of system (2) then both $R_{1}(\mathbf{u}(-t))$ and $R_{2}(\mathbf{u}(-t))$ are also solutions [11, 34, 14]. We say that $\mathbf{u}(t)$ is symmetric if it is invariant as a set under $R_{1}$ and/or $R_{2}$. Furthermore, we make the following further distinction of its symmetry properties: we say that $\mathbf{u}(t)$ is

• •

  $R_{1}$-symmetric when it is invariant under only $R_{1}$; here, $\mathbf{u}$ intersects $\Sigma_{1}$ transversally;
  

• •

  $R_{2}$-symmetric when it is invariant under only $R_{2}$; here, $\mathbf{u}$ intersects $\Sigma_{2}$ transversally.
  

• •

  $R^{*}$-symmetric when it is invariant under $R_{1}$ as well as $R_{2}$; here, $\mathbf{u}$ intersects both $\Sigma_{1}$ and $\Sigma_{2}$ transversally.
  

• •

  non-symmetric when it is not invariant (under either $R_{1}$ or $R_{2}$).

Furthermore, system (2) is equivariant under the state-space symmetry transformation

$$S:(u_{1},u_{2},u_{3},u_{4})\rightarrow(-u_{1},-u_{2},-u_{3},-u_{4}),$$

which is a point reflection in the origin $\mathbf{0}=(0,0,0,0)$ and given by $S=R_{1}\circ R_{2}=R_{2}\circ R_{1}$. Notice that invariance of a solution $\mathbf{u}(t)$ under $S$ does not necessarily imply $R^{*}$-symmetry; however, any $R^{*}$-symmetric solution is necessarily invariant under $S$.

Solitons of the GNLSE (1) are given by the $u_{1}$-component of homoclinic orbits to the origin $\mathbf{0}=(0,0,0,0)$ of system (2), as a result of the traveling-wave ansatz [6, 33]. Note that $\mathbf{0}$ is an equilibrium in the zero-energy level $H\equiv 0$ and, moreover, $\mathbf{0}\in\Sigma_{1}\cap\Sigma_{2}$. Hence, when $\mathbf{0}$ is a saddle, homoclinic orbits to it (approaching $\mathbf{0}$ in both forward and backward time) may be $R_{1}$-symmetric, $R_{2}$-symmetric, $R^{*}$-symmetric or non-symmetric [6]. 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 [16, 17, 11, 23, 14, 4, 15] to understand the mechanism in which homoclinic solutions arise and disappear as parameters are varied.

System (2) is known [6, 37, 1] to have a pair of primary $R_{1}$-symmetric homoclinic orbits to $\mathbf{0}$, which are each other’s counterparts under the reversibility given by $R_{2}$. The eigenvalues of the equilibrium $\mathbf{0}$ become complex conjugate while the pair of homoclinic orbits exists, and this is known as a Belyakov-Devaney (BD) bifurcation [16, 17]. Moreover, the primary homoclinic orbit disappears in a Hamiltonian Hopf (HH) bifurcation when $\mathbf{0}$ ceases to be a saddle and becomes a center. Without loss of generality, we may consider the $(\beta_{2},\mu)$-plane for $\beta_{4}=-1$ and $\gamma=1$ to find the relevant bifurcations, which all occur along semi-parabolas [6]. Figure 1(a) shows how the curves BD and HH together form a single (analytically known) parabola that divides the $(\beta_{2},\mu)$-plane into region I to III [6, 39]. In regions I and II the equilibrium $\mathbf{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 [6, 37, 1]. In region III, the equilibrium $\mathbf{0}$ is a center with purely imaginary eigenvalues and, hence, can no longer support homoclinic or other connecting orbits.

The Belyakov-Devaney bifurcation [16, 17] implies that infinitely many homoclinic solutions of different symmetry types exist to the right of the curve BD in Fig. 1(a) [6, 15]. 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 $\mathbf{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 $\beta_{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. 3, 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 [6]; 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 [41, 10]) and note that this phenomenon is also observed in a certain parameter regime of the Lugiato-Lefever equation [34].

In this paper, we establish the existence in the GNLSE (1) of heteroclinic connections between saddle periodic orbits in the same energy surface, which we refer to as PtoP connections. This type of connecting orbit also persists with parameters and plays an important role in the GNLSE for two main reasons. First of all, PtoP connections are organizing centers that ‘generate’ associated families of homoclinic orbits to saddle periodic orbits. Each such homoclinic orbit corresponds to a soliton with non-decaying periodic tails. This type of solution to a ‘periodic background’ is referred to as a generalized soliton in the context of optics [3]. Much like EtoP connections organize solitons with decaying tails, PtoP connections are shown to give rise to infinite families of generalized solitons with different symmetry properties.

Secondly, PtoP connections can be found in the zero-energy level. This 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 $\mathbf{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 $\mathbf{0}$ via another EtoP connection. In contrast to the solitons considered before [6], which visit a single periodic orbit, we refer to these solitons organized by PtoP connections as multi-oscillation solitons. As we will show, they are also created at BD, may have different symmetry properties, and disappear in fold bifurcations before HH is reached. Hence, multi-oscillation solitons are an integral part of the phenomenon of BD-truncated homoclinic snaking in system (2).

Similarly, PtoP connections in the zero-energy surface arise in pairs from the Belyakov-Devaney bifurcation BD and disappear in fold bifurcations before HH is reached. 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 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 [6, 34]. On the other hand, PtoP connections are not confined to the zero-energy level, and 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 [22].

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 [16, 17]. 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 $\beta_{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 $\mathbf{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 $\mathbf{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) [30, 6]. In an approach called Lin’s method [31], 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 [21, 6] to system (2). All these computations are performed with the continuation package Auto-07p [18] and its extension HomCont [13].

The paper is structured as follows. In Sec. 2, we introduce the three basic surfaces $\mathcal{S}_{1}^{\pm}$ and $\mathcal{S}_{*}$ of periodic orbits that accumulate on the primary homoclinic orbit and its $R_{2}$-counterpart. In Sec. 3, we then briefly summarise previous results [6] regarding how EtoP connections to specific periodic orbits in these surfaces organize associated families of homoclinic orbits to $\mathbf{0}$; here, we also show a one-parameter bifurcation diagram that clarifies over which $\beta_{2}$-ranges they exist. This sets the stage for the discussion in Sec. 4 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. 5 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 $\mathbf{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 $\beta_{2}$-ranges of existence. We conclude in Sec. 6 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 $\beta_{2}=0.4$, this pair of primary homoclinic solutions gives rise to three basic surfaces; here and throughout, we fix, without loss of generality, $\beta_{4}=-1$ and $\gamma=1$, and set $\mu=1$ since all bifurcation curves in the $(\beta_{2},\mu)$-plane are parabolas; note here that Eq. (1) can be rescaled appropriately for different signs of $\beta_{2}$ and $\beta_{4}$, as was shown previously [6]. 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 $\mathcal{S}_{1}^{-}$ and $\mathcal{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 $\mathbf{0}$ and the pair of primary homoclinic orbits, which are also shown in Fig. 2(b) and (c). Note that $\mathcal{S}_{1}^{-}$ and $\mathcal{S}_{1}^{+}$ are both $R_{1}$-symmetric and $R_{2}$-counterparts of one another [6]. Moreover, they each extend over only a bounded range of $H$, where the largest value of $H$ occurs at the equilibria $\mathbf{E^{-}}$ and $\mathbf{E^{+}}$. These are centers at $\mathbf{E_{\pm}}=\left(\pm\sqrt{\frac{\mu}{\gamma}},0,0,0\right)$ (so only exists if $\mu\gamma>0$) and each others $R_{2}$-counterparts with $H(\mathbf{E_{\pm}})=-6\mu^{2}/(\gamma\beta_{4})$. The union of the pair of primary homoclinic orbits is the limit of a third, $R_{*}$-symmetric surface $\mathcal{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 [6].

We refer to the the surfaces $\mathcal{S}_{1}^{\pm}$ and $\mathcal{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}\leq 0$ up to the section $\Sigma=\{(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 $\Sigma$, is obtained from that shown by applying the transformation $R_{1}$. In Fig. 3, this corresponds to a reflection in the $\Sigma$-plane; compare with Fig. 2(b)–(d).

This cutaway representation allows us to compute and show the intersection curves of the surfaces $\mathcal{S}_{1}^{\pm}$ and $\mathcal{S}_{*}$ with the plane $\Sigma$, which is very useful for understanding the geometry of these surfaces. The gray line in $\Sigma$ indicates the zero-energy level with $u_{2}=0$, which contains the basic homoclinic orbits. As Fig. 3 illustrates, the basic surfaces $\mathcal{S}_{1}^{\pm}$ and $\mathcal{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 $\mathbf{E_{\pm}}$ and decreasing the energy $H$, one encounters in the zero-energy level the $R_{1}$-symmetric and $R_{2}$-related pair of periodic orbits $\Gamma_{1}^{-}$ on $\mathcal{S}_{1}^{-}$ and $\Gamma_{1}^{+}$ on $\mathcal{S}_{1}^{+}$. Similarly, starting from large negative $H$ and increasing the energy, one encounters in the zero-energy level the $R_{*}$-symmetric periodic orbit $\Gamma_{*}$ on $\mathcal{S}_{*}$. The three periodic orbits $\Gamma_{1}^{\pm}$ and $\Gamma_{*}$ 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 $\Gamma_{1}^{\pm}$ and $\Gamma_{*}$ 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 $\Sigma$.

In this work, we demonstrate with the specific examples of the primary periodic orbits $\Gamma_{1}^{\pm}$ and $\Gamma_{*}$ for $\beta_{2}=0.4$ how connecting orbits between them give rise to PtoP connection, which in turn generate generalized and multi-oscillation solitons. Indeed, the same is true for the further periodic orbits of $\mathcal{S}_{1}^{\pm}$ and $\mathcal{S}_{*}$ with $H=0$, of which there are countably infinitely many. In contrast to the primary periodic orbits, further periodic orbits of $\mathcal{S}_{1}^{\pm}$ and $\mathcal{S}_{*}$ in the zero-energy level are multi-loop orbits with an increasing number of additional intersection points with the (three-dimensional) section $\Sigma$. How they arise from the associated structure of $\mathcal{S}_{1}^{\pm}$ and $\mathcal{S}_{*}$, and how these surfaces depend on the parameter $\beta_{2}$ are interesting subjects beyond the scope of this paper that will be discussed elsewhere [7].

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 $\mathbf{0}$ to the primary periodic orbits $\Gamma_{1}^{\pm}$ and $\Gamma_{*}$ [6], again for $\beta_{2}=0.4$ (and $\mu=1$). The dependence of the EtoP connections and associated solitons on $\beta_{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 $\Gamma_{*}$ and to $\Gamma_{1}^{\pm}$, together with their corresponding symmetric counterparts, generate heteroclinic cycles that organize homoclinic orbits to $\mathbf{0}$ under mild conditions [32, 12, 6, 33]. These homoclinic orbits closely follow an EtoP connection and make a number of loops around the respective periodic orbit before converging back to $\mathbf{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 $\Gamma_{*}$ and associated homoclinic solutions. Panels (a) and (b) show two EtoP connections from $\mathbf{0}$ to $\Gamma_{*}$ that exist simultaneously for $\beta_{2}=0.4$. Here, panels (a1) and (b1) show the EtoP connections in $(u_{1},u_{2},H)$-space together with a cutaway view of the surface $\mathcal{S}_{*}$; and panels (a2) and (b2) show their $u_{1}$-traces, which illustrates that the EtoP connections converge backward in time to $\mathbf{0}$ and forward in time to $\Gamma_{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 $\Gamma_{*}$ 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 $\Gamma_{*}$ to $\mathbf{0}$. Hence, the EtoP connections in Fig. 4(a) and (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 [6]. Some representative examples are shown in Fig. 4(c)–(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 Fig. 4(c)–(f), these transformations are geometrically: reflection in the $u_{1}$ axis, and rotation by $\pi$ around the origin of the $(t,u_{2})$-plane, respectively.

In the same style, Fig. 5 shows the EtoP connections to $\Gamma_{1}^{+}$ and associated homoclinic solutions that exist simultaneously. The difference is now that $\Gamma_{1}^{+}$ is not $R_{2}$-symmetric, and this means that the two distinct EtoP connections to $\Gamma_{1}^{+}$, which are shown in panels (a) and (b) with a cutaway view of $\mathcal{S}_{1}^{+}$, can only be combined with their $R_{1}$-counterparts to generate heteroclinic cycles between $\mathbf{0}$ and $\Gamma_{1}^{+}$. (Of course, their images under the reversibility $R_{2}$ also exist, but they are heteroclinic cycles between $\mathbf{0}$ and $\Gamma_{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 $\Gamma_{1}^{+}$, and then follow the $R_{1}$-counterpart back to $\mathbf{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 [11, 20]. Hence, once found, each such object can be continued as a solution branch in a chosen parameter, which we take to be the strength $\beta_{2}$ of the quadratic dispersion. In this way, one can compute a one-parameter bifurcation diagram that summarizes over what ranges of $\beta_{2}$ which types of connections and associated solitons can be found [6]. 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 $\beta_{2}$-range includes the entire region II bounded by the bifurcations BD and HH, which are indicated by vertical lines and lie at $\beta_{2}\approx\pm 0.8164$ (for $\mu=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 [18]. 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 $\Gamma_{*}$ from Figs. 4(a) and (b) lie on branches that emerge from BD and meet at a fold point at $\beta_{2}\approx 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 $\mathbf{0}$ and $\Gamma_{*}$. 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 very similar. 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 $\mathbf{0}$ and $\Gamma_{*}$. 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 $\Gamma_{1}^{+}$ from Figs. 5(a) and (b) (and also those to $\Gamma_{1}^{-}$ by $R_{2}$-symmetry) lie on a single branch that emerges from and returns to BD, with a fold point that lies at $\beta_{2}\approx 0.6756$. We refer to this fold as EtoP₁, 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 Fig. 5(c) and (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₁ 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 $\Gamma_{*}$ and $\Gamma_{1}^{+}$, and show that they emerge from the Belyakov-Devaney bifurcation BD and exist over a considerable $\beta_{2}$-range in region II.

Figure 7 shows in projection onto $(u_{1},u_{1},u_{3})$-space a pair of PtoP connections between $\Gamma_{*}$ and $\Gamma_{1}^{+}$ for $\beta_{2}=0.4$. Each of these two PtoP connections is clearly seen to converge backward in time to $\Gamma_{*}$ and forward in time to $\Gamma^{+}_{1}$, and they are not related by symmetry. While the PtoP connection in Fig. 7(a) transitions swiftly from $\Gamma_{*}$ to $\Gamma^{+}_{1}$, the one in panel (b) has almost a complete loop near and around a third saddle periodic orbit ${}^{A}\Gamma^{+}_{1}$ that also lies on the surface $\mathcal{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 $\Gamma_{*}$ or $\Gamma_{1}^{\pm}$ in this case. In the GNLSE, they correspond to generalized solitons [3] with non-decaying periodic tails, which are not confined to the zero-energy level. In Sec. 5, 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 $\mathbf{0}$ that generate families of multi-oscillation solitons.

In the zero-energy level one finds the PtoP connections between $\Gamma_{*}$ and $\Gamma^{+}_{1}$ from Fig. 7, as well as their $R_{1}$- and $R_{2}$-counterparts that involve $\Gamma^{-}_{1}$. Together they are involved in creating more complicated heteroclinic PtoP cycles between the three periodic orbits $\Gamma_{*}$, $\Gamma^{+}_{1}$ and $\Gamma^{-}_{1}$. As we show now, these PtoP cycles can be combined with the EtoP connections from Sec. 3 effectively in any combination to create more complicated heteroclinic cycles from and back to $\mathbf{0}$. In turn, these organize a plethora of homoclinic orbits to $\mathbf{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 $\mathbf{0}$, $\Gamma_{*}$, and $\Gamma^{\pm}$ at $\beta_{2}$-values where all of them exist, such as for $\beta=0.4$ as we used throughout. As we discussed in Sec 3 and Sec. 4, 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 \raisebox{-0.9pt}{1}⃝ to \raisebox{-0.9pt}{8}⃝; here the arrows indicate the direction of time and reference to previous figures indicates already the encountered connections. The edges labeled \raisebox{-0.9pt}{7}⃝ and \raisebox{-0.9pt}{8}⃝ represent pairs of PtoP connections between the $R_{2}$-counterparts $\Gamma^{+}_{1}$ and $\Gamma^{-}_{1}$, which are jointly represented by a single bi-colored circle. Furthermore, in this labeling, the even connections \raisebox{-0.9pt}{2}⃝, \raisebox{-0.9pt}{4}⃝, \raisebox{-0.9pt}{6}⃝ and \raisebox{-0.9pt}{8}⃝ correspond to the $R_{1}$- or $R_{2}$-counterparts of the odd connections \raisebox{-0.9pt}{1}⃝, \raisebox{-0.9pt}{3}⃝, \raisebox{-0.9pt}{5}⃝ and \raisebox{-0.9pt}{7}⃝, respectively.

Any path in the directed graph in Fig. 11 that starts and ends at $\mathbf{0}$ represents certain heteroclinic cycles with associated infinite families of homoclinic orbits to $\mathbf{0}$ of the different symmetry types, which are all solitons of the GNLSE. Notice that the families of homoclinic solutions that were previously found [6] and are summarized in Sec. 3 are organized by the cycles that are formed by only the connections \raisebox{-0.9pt}{1}⃝ and \raisebox{-0.9pt}{2}⃝, and by \raisebox{-0.9pt}{3}⃝ and \raisebox{-0.9pt}{4}⃝. That is, they reach $\Gamma_{*}$ or $\Gamma^{+}_{1}$ along the EtoP connection \raisebox{-0.9pt}{1}⃝ or \raisebox{-0.9pt}{3}⃝, loop multiple type times close to the respective periodic orbit, and then converge back to $\mathbf{0}$ along the EtoP connection \raisebox{-0.9pt}{2}⃝ or \raisebox{-0.9pt}{4}⃝. Importantly, the existence of the PtoP connections \raisebox{-0.9pt}{5}⃝ to \raisebox{-0.9pt}{8}⃝ that connect the periodic orbits $\Gamma_{*}$ and $\Gamma^{\pm}_{1}$ allows for more complicated homoclinic solutions to $\mathbf{0}$ to be constructed. Each arrow represents two distinct connections from one saddle object to another, which makes this mechanism more complex than Fig 11 suggests.

Our numerical evidence supports the observation/conjecture that for any path from $\mathbf{0}$ back to itself in Fig 11, the corresponding EtoP cycles with associated families of homoclinic solutions to $\mathbf{0}$ indeed exist in system (2). To illustrate this statement, we consider two specific paths and compute and show the corresponding families of homoclinic solutions. We first consider the heteroclinic cycles given by the edge sequence

$$\mathbf{0}\xrightarrow{\textrm{\raisebox{-0.5pt}{1}⃝}}\leavevmode\nobreak\ \Gamma^{*}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{5}⃝}}\leavevmode\nobreak\ {\Gamma^{+}_{1}}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{6}⃝}}\leavevmode\nobreak\ \Gamma^{*}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{2}⃝}}\leavevmode\nobreak\ \mathbf{0},$$

which we say is of type (\raisebox{-0.9pt}{1}⃝, \raisebox{-0.9pt}{5}⃝, \raisebox{-0.9pt}{6}⃝, \raisebox{-0.9pt}{2}⃝). Since the periodic orbit $\Gamma^{+}_{1}$ in the center is only $R_{1}$-symmetric, this cycle only organizes families of $R_{1}$-symmetric and $R_{1}$-symmetry-broken homoclinic solutions.

Secondly, we consider the heteroclinic cycles given by the edge sequence

$$\mathbf{0}\xrightarrow{\textrm{\raisebox{-0.5pt}{3}⃝}}\leavevmode\nobreak\ \Gamma^{+}_{1}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{6}⃝}}\leavevmode\nobreak\ {\Gamma_{*}}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{5}⃝}}\leavevmode\nobreak\ \Gamma^{+}_{1}/\Gamma^{-}_{1}\leavevmode\nobreak\ \xrightarrow{\textrm{\raisebox{-0.5pt}{4}⃝}}\leavevmode\nobreak\ \mathbf{0},$$

which are of type (\raisebox{-0.9pt}{3}⃝, \raisebox{-0.9pt}{6}⃝, \raisebox{-0.9pt}{5}⃝, \raisebox{-0.9pt}{4}⃝). Since now the periodic orbit $\Gamma_{*}$ at the center is $R_{*}$-symmetric, this cycle organizes families of both $R_{1}$- and $R_{2}$-symmetric homoclinic solutions and the corresponding symmetry-broken homoclinic solutions.

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

We showed that there exist infinitely many infinite families of generalized and multi-oscillation solitons of the GNLSE, in addition to the families of ‘regular’ solitons that were found before [6]. To achieve this, we used a dynamical system approach to translate solitons of the GNLSE into homoclinic orbits to the equilibrium $\mathbf{0}$ in the zero-energy level of the ODE system (2), which is Hamiltonian and features two reversible symmetries, $R_{1}$ and $R_{2}$. The new kinds of solitons crucially involve transitions between periodic solutions, which are given by PtoP connections between saddle-type periodic orbits of the ODE. For the specific example of the primary periodic orbits $\Gamma_{*}$ and $\Gamma_{1}^{\pm}$ 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 generalized solitons with non-decaying oscillating tails; notably, these PtoP connections and associated generalized solitons are not confined to the zero-energy level. Furthermore, PtoP connections in the zero-energy level can be combined with EtoP connections from $\mathbf{0}$ to the respective periodic orbits into heteroclinic cycles. As we argued, these generate infinite families of homoclinic solutions to $\mathbf{0}$, which are multi-oscillation solitons of the GNLSE with 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 solitons can be continued in parameters. Specifically, we showed by changing the strength $\beta_{2}$ of the quadratic dispersion that all solitons are organized as BD-truncated homoclinic snaking, even the generalized solitons that correspond to homoclinic solutions to periodic orbits. This means that, for all families of solitons, pairs of branches start at the Belyakov-Daveney bifurcation BD and meet at fold bifurcations that are also symmetry-breaking bifurcation of the respective solutions. The folds of a given family accumulate on the folds of the respective EtoP or ProP connections that are organizing it. Note, in particular, that all these solutions of the GNLSE exist for $\beta_{2}=0$, that is, for the case of pure quartic dispersion.

For heteroclinic cycles involving only the equilibrium $\mathbf{0}$ and the primary periodic orbits $\Gamma_{*}$ and $\Gamma_{1}^{\pm}$, we presented a directed graph with numbered edges to encode the possibilities: any edge sequence in this graph from and back to $\mathbf{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 generally holds. 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 $\mathcal{S}_{*}$ and $\mathcal{S}_{1}^{\pm}$ 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 $\mathcal{S}_{*}$ and $\mathcal{S}_{1}^{\pm}$ and associated edges that represent EtoP connections from $\mathbf{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 of all, we already alluded to the geometry of the basic surfaces $\mathcal{S}_{*}$ and $\mathcal{S}_{1}^{\pm}$, which we only considered here for the fixed parameter value $\beta_{2}=0.4$. As $\beta_{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 [7]. Secondly, 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) [33, 26]. 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 [40, 24]. Another issue is the stability of the different solitons as solutions of the PDE. Only the primary single-hump soliton is stable [37, 39], while all the other solitons we discovered appear to be unstable (but some are only weakly unstable [6]). 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 [8].

Finally, 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 [36] have shown the feasibilities of creating waveguides with higher even-order dispersions, which are described by further generalizations of the NLSE with corresponding non-zero dispersion terms of order six, eight, ten, or even higher. Other examples of interesting PDEs in this context are the widely studied class of the NLSE describing solitary wave formation in inhomogeneous media [29, 27, 19, 28], and the Lugiato-Lefever equation [34], which both are known to feature homoclinic snaking.

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

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